Main Page   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Compound Members   File Members  

/Go4EventServer/random-coll.c

Go to the documentation of this file.
00001 //---------------------------------------------------------------
00002 //        Go4 Release Package v2.10-5 (build 21005) 
00003 //                      03-Nov-2005
00004 //---------------------------------------------------------------
00005 //       The GSI Online Offline Object Oriented (Go4) Project
00006 //       Experiment Data Processing at DVEE department, GSI
00007 //---------------------------------------------------------------
00008 //
00009 //Copyright (C) 2000- Gesellschaft f. Schwerionenforschung, GSI
00010 //                    Planckstr. 1, 64291 Darmstadt, Germany
00011 //Contact:            http://go4.gsi.de
00012 //----------------------------------------------------------------
00013 //This software can be used under the license agreements as stated
00014 //in Go4License.txt file which is part of the distribution.
00015 //----------------------------------------------------------------
00016 /*----------------------------------------------------------------------
00017    DISCLAIMER
00018 
00019    This code is provided ``as is'' without express or implied
00020    warranty. The authors do not and cannot warrant the performance
00021    of the code provided, or the results that may be obtained by
00022    its use or its fitness for any specific use. In no event shall
00023    the authors be liable for any loss or damages that arise
00024    from the use of this code.
00025 
00026    BUGS:
00027    Report to: moudgill@cs.cornell.edu
00028 
00029 ----------------------------------------------------------------------*/
00030 
00031 /*----------------------------------------------------------------------
00032    This file contains the C implementation of various random number
00033    generators as described in
00034    Seminumerical Algorithms, Vol. 2, The Art of Computer Programming,
00035    by D.E. Knuth.
00036 
00037    Included are:
00038       Beta
00039       Bionomial
00040       Exponential
00041       Gamma
00042       Geometric
00043       Normal
00044       Uniform
00045 
00046    This was implemented by Mayan Moudgill, Tushar Deepak Chandra
00047    and Sandeep Jain in 1984 as part of the project requirement
00048    for the degree of Bachelor of Technology in Computer Science
00049    and Engineering at the Indian Institute of Technology, Kanpur.
00050 
00051 ----------------------------------------------------------------------*/
00052 
00053 /*----------------------------------------------------------------------
00054 
00055    This file ccntains the beta random number generator.
00056    The algorithm used depends on the size of the inputs;
00057    for small a and b, we use the method of M.D. Johnk
00058    (Metrica 8, 1964, pg 5-15), in which we repeatedly
00059    take two uniformly distributed numbers till the sum of
00060    their 1/a and 1/b powers is less than 1. Then we return
00061    one number divided by their sum.
00062    Otherwise, we generate two gamma distributed numbers,
00063    with orders a and b and return one number divided by their
00064    sum.
00065 
00066    Reference: The Art of Computer Programming, D.E. Knuth 2nd ed., 3.4.1e
00067 
00068 ----------------------------------------------------------------------*/
00069 
00070 /*----------------------------------------------------------------------
00071 
00072    p_dBeta() returns beta distributed random numbers
00073    with mean parameters a and b.
00074 
00075    Input:
00076       a,b   : positive doubles.
00077       seed  : a pointer to an unsigned integer. This is used to
00078          store previous uniform variate calculated for
00079          calculation of next uniform variate.
00080          It must not be changed by calling program.
00081 
00082    Output: an beta distributed double, with parameters a and b.
00083 
00084 ----------------------------------------------------------------------*/
00085 
00086 extern double p_dUniform ();
00087 extern double pow ();
00088 extern double p_dGammaGen ();
00089 #define BETA_TURNING_POINT 2
00090 
00091 double   p_dBeta (a,b,seed)
00092 double   a,b;
00093 unsigned int   *seed;
00094 {
00095    if (a*b < BETA_TURNING_POINT)
00096    {
00097       double   u1,u2,p1,p2,x;
00098 
00099       p1 = 1/a;
00100       p2 = 1/b;
00101       do
00102       {
00103          u1 = pow (p_dUniform(seed),p1);
00104          u2 = pow (p_dUniform(seed),p2);
00105       } while ((x=u1+u2) > 1);
00106       return (u1/x);
00107    }
00108    else
00109    {
00110       double   x,y;
00111 
00112       x = p_dGammaGen (a,seed);
00113       y = p_dGammaGen (b,seed);
00114       return (x/(x+y));
00115    }
00116 }
00117 
00118 /*------------------------------------------------------------------------------
00119 
00120    This file contains the binomial random number generator.
00121    Unfortunately, here I actually have to simulate it for want of
00122    a better method. Thus if a binomial is required with
00123    mean .3 and n=50,fifty uniform variables are required!
00124 
00125    Reference: The Art of Computer Programming, D.E. Knuth 2nd ed., 3.4.1f
00126 
00127 ------------------------------------------------------------------------------*/
00128 
00129 /*------------------------------------------------------------------------------
00130    p_iBionamial() returns bionomially distributed random numbers
00131    with n trials of events with a probability p.
00132 
00133    Input:
00134       n     : the numer of trials
00135       p     : the probability of success in each trial
00136       seed  : a pointer to an unsigned integer. This is used to
00137          store previous uniform variate calculated for
00138          calculation of next uniform variate.
00139          It must not be changed by calling program.
00140 
00141    Output: a bionomially distributed integer, with n trials of
00142       events with probability p.
00143 
00144 ------------------------------------------------------------------------------*/
00145 
00146 extern double   p_dUniform ();
00147 
00148 int   p_iBinomial (p,n,seed)
00149 unsigned int   *seed;
00150 int   n;
00151 double   p;
00152 {
00153    int   i,ans=0;
00154 
00155    for (i=0;i<n;i++)
00156       if (p_dUniform (seed) <= p)
00157          ans++;
00158 
00159    return (ans);
00160 }
00161 /*------------------------------------------------------------------------------
00162 
00163    This file contains exponential random number function.The method
00164    is to get a uniform random number U in (0,1] and return -ln(U).
00165    Care is taken to aviod U=0;
00166 
00167    Reference: The Art of Computer Programming, D.E. Knuth 2nd ed., 3.4.1d
00168 
00169 ------------------------------------------------------------------------------*/
00170 
00171 /*------------------------------------------------------------------------------
00172    p_dExponential() returns exponentially distributed random numbers
00173    with mean lambda.
00174 
00175    Input:
00176       lambda: the mean.
00177       seed  : a pointer to an unsigned integer. This is used to
00178          store previous uniform variate calculated for
00179          calculation of next uniform variate.
00180          It must not be changed by calling program.
00181 
00182    Output: an exponentially distributed double, with mean lambda
00183 
00184 ------------------------------------------------------------------------------*/
00185 
00186 extern double p_dUniform();
00187 extern double log();
00188 
00189 double   p_dExponential (lambda,seed)
00190 double   lambda;
00191 unsigned int   *seed;
00192 {
00193    double   u= p_dUniform( seed);
00194 
00195    while (u == (double)0)
00196       u= p_dUniform( seed);
00197 
00198    return (-lambda*log(u));
00199 }
00200 /*----------------------------------------------------------------------
00201 
00202    The standard Gamma distribution has been implemented using
00203    two functions:
00204       <1> For integer orders
00205       <2> For order in range (0..1)
00206 
00207    A gamma distribution of order X+Y is the same as that of the
00208    sum of two independent gamma distributions of orders X and Y
00209    respectivley. Hence, for general orders, we split the order
00210    into it's integer and fractional parts, calculate their gamma
00211    variates, and sum the variates obtained.
00212 
00213    Reference: The Art of Computer Programming, D.E. Knuth 2nd ed., 3.4.1e
00214 ----------------------------------------------------------------------*/
00215 
00216 /*----------------------------------------------------------------------
00217 
00218    There are 2 methods for Gamma with integer A. The fisrt is to
00219    get the sum of A independent exponential variates, is quick,
00220    but grows linearly with A. The second method gets a closed
00221    form and is more or less independent of A in speed.  Therefore
00222    for A>GAMMA_TURNING_POINT we should switch methods.
00223 
00224    GAMMA_TURNING_POINT is 7 for a Sun SS-10.
00225 
00226    Reference: D.E.D.E. Knuth,The art of computer programming,Vol 2
00227       Chapter 3.4.1E algo A (p.129)
00228 
00229 ----------------------------------------------------------------------*/
00230 
00231 #define   GAMMA_TURNING_POINT   7
00232 extern double   p_dUniform ();
00233 extern double   log ();
00234 extern double   exp ();
00235 extern double   sqrt ();
00236 extern double   tan ();
00237 extern double   pow();
00238 extern double   floor();
00239 
00240 double   p_dGammaInt (a,seed)
00241 int   a;
00242 unsigned int   *seed;
00243 {
00244    if (a>GAMMA_TURNING_POINT)
00245    {
00246       double   x,y,v,sqrt2a1,asub1,y2;
00247 
00248       sqrt2a1 = sqrt ((double)(2*a - 1));
00249       asub1 = a-1;
00250       while (1)
00251       {
00252          y = tan (3.1415926536*p_dUniform (seed));
00253          y2 = y*sqrt2a1;
00254          if ((x = y2 + asub1) <= 0)
00255             continue;
00256          v = p_dUniform (seed);
00257          if (v <= (1+y*y)*exp((asub1)*log(x/asub1)-y2))
00258             return (x);
00259       }
00260    }
00261    else
00262    {
00263       double ans=1;
00264       while (a--)
00265          ans *= p_dUniform (seed);
00266       return (-log(ans));
00267    }
00268 }
00269 
00270 /*----------------------------------------------------------------------
00271 
00272    GammaGen() takes any order A>0, and generates a gamma variate
00273    as described above.
00274 
00275    The function incorporated in GammaGen() for generating gamma
00276    variates in range (0..1) was given by Fishman. It is taken
00277    from:
00278    W.J.Kennedy,Jr & James E. Gentle, Statistical Computing, pg214
00279 
00280 ----------------------------------------------------------------------*/
00281 
00282 double p_dGammaGen (A,seed)
00283 double A;
00284 unsigned *seed;
00285 {
00286 if (A <= 5)
00287    {
00288    double b,p,inva,res=0;
00289    double EXP= 2.71828182;
00290    double a,f;
00291    int g;
00292 
00293    f= floor( A);
00294    a= A-f;
00295    g= (int) f;
00296 
00297    if( a!=0)
00298    {
00299    inva= 1/a;
00300    b= 1.0 + a/EXP;
00301 
00302    while(1)
00303    {
00304    p= b*p_dUniform(seed);
00305    if( p<=1.0 )
00306    {
00307    res= pow(p, inva);
00308    if( p_dUniform(seed) <= exp( -res)) break;
00309    }
00310    else
00311    {
00312    res= - log( (b-p)/a);
00313    if( p_dUniform(seed) <= 1.0/pow(res, 1-a)) break;
00314    }
00315    }
00316    }
00317 
00318    if ( g>0)
00319    res += p_dGammaInt(g,seed);
00320 
00321    return(res );
00322    }
00323 else /*Use the same function as integer*/
00324 {
00325    double   x,y,v,sqrt2a1,asub1,y2;
00326 
00327    sqrt2a1 = sqrt (2*A - 1);
00328    asub1 = A-1;
00329    while (1)
00330    {
00331       y = tan (3.1415926536*p_dUniform (seed));
00332       y2 = y*sqrt2a1;
00333       if ((x = y2 + asub1) <= 0)
00334          continue;
00335       v = p_dUniform (seed);
00336       if (v <= (1+y*y)*exp((asub1)*log(x/asub1)-y2))
00337          return (x);
00338    }
00339 }
00340 }
00341 /*------------------------------------------------------------------------------
00342 
00343    This file contains functions to generate geometric random numbers.
00344       P(1) = p
00345       P(2) = p(1-p)
00346       P(i) = p(1-p)^(i-1)
00347    Keep generating a uniform random number till it is lesser than p.
00348       OR
00349    Take ceil (log (U)/log (p)).Though this is slower,we will select this
00350    if p < .20 because as p falls,the iterative algo tends to slow
00351 down and
00352    they are equally fast at p = GEOMETRIC_TURNING_POINT.
00353    Note that when p is large,the first method is nearly twice as
00354 fast as
00355    the second but when it can be made arbitarily slow by reducing p.
00356 
00357    GEOMETRIC_TURNING_POINT is 0.38 for a Sun SS-10.
00358 
00359    ERROR HANDLING: if p>1 or p<=0 then we have an error
00360 
00361 ------------------------------------------------------------------------------*/
00362 
00363 extern double p_dUniform();
00364 extern double log();
00365 extern double ceil();
00366 #define   GEOMETRIC_TURNING_POINT  0.38
00367 
00368 int   p_iGeometric (p,seed)
00369 double   p;
00370 unsigned int   *seed;
00371 {
00372    int   i=1;
00373 
00374    if (p < GEOMETRIC_TURNING_POINT)
00375       return ((int) ceil(log (p_dUniform(seed))/log (1-p)));
00376 
00377    while (p_dUniform(seed) > p)
00378       i++;
00379    return (i);
00380 }
00381 
00382 /*------------------------------------------------------------------------------
00383 
00384    This file contains functions to generate normally distributed
00385    variables. It uses the polar method developed by G.E.P. Box, M.E.
00386    Muller and G. Marsaglia.
00387 
00388    Reference: D.E.D.E. Knuth, The Art of Computer Programming Ed. 2, 3.4.1C
00389 
00390 ------------------------------------------------------------------------------*/
00391 
00392 /*------------------------------------------------------------------------------
00393    p_dNormal() returns a normally distributed double variable with
00394 mean
00395    mean, standard deviation  sigma.
00396 
00397 
00398 ------------------------------------------------------------------------------*/
00399 
00400 extern double p_dUniform();
00401 extern double log();
00402 extern double sqrt();
00403 
00404 double p_dNormal (mean,sigma,seed)
00405 double   mean,sigma;
00406 unsigned int   *seed;
00407 {
00408 
00409 double u1, u2, v1, v2, s;
00410 static char flag= 0;
00411 static double x1, x2;
00412 
00413 if( flag )
00414 {
00415 flag=0;
00416 return(x2*sigma + mean);
00417 }
00418 
00419 do {
00420 u1= p_dUniform(seed);
00421 u2= p_dUniform(seed);
00422 
00423 v1= 2* u1 -1;
00424 v2= 2* u2 -1;
00425 s= v1*v1 + v2*v2;
00426 } while( s>=1);
00427 
00428 s= sqrt( -2*log(s)/s);
00429 x1= v1* s;
00430 x2= v2* s;
00431 
00432 flag= 1;
00433 return( x1*sigma + mean);
00434 
00435 }
00436 extern double exp();
00437 extern double p_dUniform();
00438 
00439 int p_iPoisson(mu,seed)
00440 double mu;
00441 unsigned *seed;
00442 {
00443 double total=1, till;
00444 int count= -1   ;
00445 
00446 till= exp( -mu);
00447 while( total > till)
00448 {
00449 total *= p_dUniform(seed);
00450 count++;
00451 }
00452 
00453 return(count);
00454 }
00455 
00456 /*------------------------------------------------------------------------------
00457 
00458    This file contains functions for the generation of uniformly
00459    distributed random numbbers. The linear congruential method is
00460    used. Function used is
00461 
00462       X[n+1]= ( a * X[n] + c ) mod m
00463 
00464    with
00465       a= 1664525 ,
00466       c= 1 ,
00467       m= 2^32 .
00468 
00469    Reference : D.E. Knuth, The Art of Computer Programming, Chapter 3.
00470 
00471 ------------------------------------------------------------------------------*/
00472 
00473 
00474 #define      _Mult   1664525
00475 #define      _Cons   1
00476 #define    _Mask   0xFFFF
00477 #define      _Lo(X)   (X&_Mask)      /* the 16 LSB of X */
00478 #define      _Hi(X)   ((X>>16)&_Mask)      /* the 16 MSB of X (if 32 bit)*/
00479 
00480 
00481 /*------------------------------------------------------------------------------
00482 
00483    p_iUniform() and p_dUniform generate uniformly distributed random
00484    numbers, the first returning an unsigned integer in range [0,2^33-1],
00485    and the other a double in range [0.0, 1.0)
00486 
00487    Input: the pointer to an unsigned number; this is X[n]. This
00488       location is used to store the value of the old seed,
00489       hence should not be changed by calling routine.
00490 
00491    Output: a random number, X[n+1] calculated from Input using given
00492       formula, converted into correct type (unsigned or double).
00493 
00494 
00495 ------------------------------------------------------------------------------*/
00496 
00497 unsigned p_iUniform ( seed)
00498 unsigned *seed;
00499 {
00500 unsigned lo, hi;
00501 
00502 lo= _Lo(_Lo(*seed) * _Lo(_Mult) + _Cons);
00503 hi= _Lo(_Hi(*seed) * _Lo(_Mult))+_Lo(_Hi(_Mult) * _Lo(*seed))+
00504     _Hi(_Lo(*seed) * _Lo(_Mult) + _Cons);
00505 
00506 *seed= (hi<<16 | lo);
00507 return( *seed );
00508 }
00509 
00510 
00511 static unsigned   Scal=0xFFFFFFFF;
00512 double p_dUniform( seed)
00513 unsigned * seed;
00514 {
00515 unsigned lo, hi;
00516 
00517 lo= _Lo(_Lo(*seed) * _Lo(_Mult) + _Cons);
00518 hi= _Lo(_Hi(*seed) * _Lo(_Mult))+_Lo(_Hi(_Mult) * _Lo(*seed))+
00519     _Hi(_Lo(*seed) * _Lo(_Mult) + _Cons);
00520 
00521 *seed= (hi<<16 | lo);
00522 return( ((double)*seed)/Scal );
00523 }
00524 /*
00525 
00526 --
00527 --------------------------------------------------------------------------
00528 Mark Feldman                               \\_==_o    Skydivers do
00529 E-mail : u914097@student.canberra.edu.au        \\    it at 120mph....
00530 --------------------------------------------------------------------------
00531 */
00532 
00533 
00534 //----------------------------END OF GO4 SOURCE FILE ---------------------

Generated on Tue Nov 8 10:55:53 2005 for Go4-v2.10-5 by doxygen1.2.15