GSI Object Oriented Online Offline (Go4)  GO4-6.2.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
random-coll.c
Go to the documentation of this file.
1 // $Id: random-coll.c 3499 2022-01-21 16:42:51Z linev $
2 //-----------------------------------------------------------------------
3 // The GSI Online Offline Object Oriented (Go4) Project
4 // Experiment Data Processing at EE department, GSI
5 //-----------------------------------------------------------------------
6 // Copyright (C) 2000- GSI Helmholtzzentrum fuer Schwerionenforschung GmbH
7 // Planckstr. 1, 64291 Darmstadt, Germany
8 // Contact: http://go4.gsi.de
9 //-----------------------------------------------------------------------
10 // This software can be used under the license agreements as stated
11 // in Go4License.txt file which is part of the distribution.
12 //-----------------------------------------------------------------------
13 
14 
15 /*----------------------------------------------------------------------
16  DISCLAIMER
17 
18  This code is provided ``as is'' without express or implied
19  warranty. The authors do not and cannot warrant the performance
20  of the code provided, or the results that may be obtained by
21  its use or its fitness for any specific use. In no event shall
22  the authors be liable for any loss or damages that arise
23  from the use of this code.
24 
25  BUGS:
26  Report to: moudgill@cs.cornell.edu
27 
28 ----------------------------------------------------------------------*/
29 
30 /*----------------------------------------------------------------------
31  This file contains the C implementation of various random number
32  generators as described in
33  Seminumerical Algorithms, Vol. 2, The Art of Computer Programming,
34  by D.E. Knuth.
35 
36  Included are:
37  Beta
38  Bionomial
39  Exponential
40  Gamma
41  Geometric
42  Normal
43  Uniform
44 
45  This was implemented by Mayan Moudgill, Tushar Deepak Chandra
46  and Sandeep Jain in 1984 as part of the project requirement
47  for the degree of Bachelor of Technology in Computer Science
48  and Engineering at the Indian Institute of Technology, Kanpur.
49 
50 ----------------------------------------------------------------------*/
51 
52 /*----------------------------------------------------------------------
53 
54  This file ccntains the beta random number generator.
55  The algorithm used depends on the size of the inputs;
56  for small a and b, we use the method of M.D. Johnk
57  (Metrica 8, 1964, pg 5-15), in which we repeatedly
58  take two uniformly distributed numbers till the sum of
59  their 1/a and 1/b powers is less than 1. Then we return
60  one number divided by their sum.
61  Otherwise, we generate two gamma distributed numbers,
62  with orders a and b and return one number divided by their
63  sum.
64 
65  Reference: The Art of Computer Programming, D.E. Knuth 2nd ed., 3.4.1e
66 
67 ----------------------------------------------------------------------*/
68 
69 /*----------------------------------------------------------------------
70 
71  p_dBeta() returns beta distributed random numbers
72  with mean parameters a and b.
73 
74  Input:
75  a,b : positive doubles.
76  seed : a pointer to an unsigned integer. This is used to
77  store previous uniform variate calculated for
78  calculation of next uniform variate.
79  It must not be changed by calling program.
80 
81  Output: an beta distributed double, with parameters a and b.
82 
83 ----------------------------------------------------------------------*/
84 
85 #include "random-coll.h"
86 #include <math.h>
87 
88 #define BETA_TURNING_POINT 2
89 
90 #define _Mult 1664525
91 #define _Cons 1
92 #define _Mask 0xFFFF
93 #define _Lo(X) ((X)&_Mask) /* the 16 LSB of X */
94 #define _Hi(X) (((X)>>16)&_Mask) /* the 16 MSB of X (if 32 bit)*/
95 
96 
97 /*------------------------------------------------------------------------------
98 
99  p_iUniform() and p_dUniform generate uniformly distributed random
100  numbers, the first returning an unsigned integer in range [0,2^33-1],
101  and the other a double in range [0.0, 1.0)
102 
103  Input: the pointer to an unsigned number; this is X[n]. This
104  location is used to store the value of the old seed,
105  hence should not be changed by calling routine.
106 
107  Output: a random number, X[n+1] calculated from Input using given
108  formula, converted into correct type (unsigned or double).
109 
110 
111 ------------------------------------------------------------------------------*/
112 
113 unsigned p_iUniform(unsigned *seed)
114 {
115 unsigned lo, hi;
116 
117 lo= _Lo(_Lo(*seed) * _Lo(_Mult) + _Cons);
118 hi= _Lo(_Hi(*seed) * _Lo(_Mult)) + _Lo(_Hi(_Mult) * _Lo(*seed)) +
119  _Hi((_Lo(*seed) * _Lo(_Mult)) + _Cons);
120 
121 *seed= (hi<<16 | lo);
122 return( *seed );
123 }
124 
125 static unsigned Scal=0xFFFFFFFF;
126 double p_dUniform(unsigned* seed)
127 {
128 unsigned lo, hi;
129 
130 lo= _Lo(_Lo(*seed) * _Lo(_Mult) + _Cons);
131 hi= _Lo(_Hi(*seed) * _Lo(_Mult)) + _Lo(_Hi(_Mult) * _Lo(*seed)) +
132  _Hi((_Lo(*seed) * _Lo(_Mult)) + _Cons);
133 
134 *seed= (hi<<16 | lo);
135 return( ((double)*seed)/Scal );
136 }
137 
138 
139 
140 
141 double p_dBeta (double a, double b, unsigned int *seed)
142 {
143  if (a*b < BETA_TURNING_POINT)
144  {
145  double u1,u2,p1,p2,x;
146 
147  p1 = 1/a;
148  p2 = 1/b;
149  do
150  {
151  u1 = pow (p_dUniform(seed),p1);
152  u2 = pow (p_dUniform(seed),p2);
153  } while ((x=u1+u2) > 1);
154  return (u1/x);
155  }
156  else
157  {
158  double x,y;
159 
160  x = p_dGammaGen (a,seed);
161  y = p_dGammaGen (b,seed);
162  return (x/(x+y));
163  }
164 }
165 
166 /*------------------------------------------------------------------------------
167 
168  This file contains the binomial random number generator.
169  Unfortunately, here I actually have to simulate it for want of
170  a better method. Thus if a binomial is required with
171  mean .3 and n=50,fifty uniform variables are required!
172 
173  Reference: The Art of Computer Programming, D.E. Knuth 2nd ed., 3.4.1f
174 
175 ------------------------------------------------------------------------------*/
176 
177 /*------------------------------------------------------------------------------
178  p_iBionamial() returns bionomially distributed random numbers
179  with n trials of events with a probability p.
180 
181  Input:
182  n : the numer of trials
183  p : the probability of success in each trial
184  seed : a pointer to an unsigned integer. This is used to
185  store previous uniform variate calculated for
186  calculation of next uniform variate.
187  It must not be changed by calling program.
188 
189  Output: a bionomially distributed integer, with n trials of
190  events with probability p.
191 
192 ------------------------------------------------------------------------------*/
193 
194 int p_iBinomial (double p, int n, unsigned int *seed)
195 {
196  int i,ans=0;
197 
198  for (i=0;i<n;i++)
199  if (p_dUniform (seed) <= p)
200  ans++;
201 
202  return (ans);
203 }
204 /*------------------------------------------------------------------------------
205 
206  This file contains exponential random number function.The method
207  is to get a uniform random number U in (0,1] and return -ln(U).
208  Care is taken to aviod U=0;
209 
210  Reference: The Art of Computer Programming, D.E. Knuth 2nd ed., 3.4.1d
211 
212 ------------------------------------------------------------------------------*/
213 
214 /*------------------------------------------------------------------------------
215  p_dExponential() returns exponentially distributed random numbers
216  with mean lambda.
217 
218  Input:
219  lambda: the mean.
220  seed : a pointer to an unsigned integer. This is used to
221  store previous uniform variate calculated for
222  calculation of next uniform variate.
223  It must not be changed by calling program.
224 
225  Output: an exponentially distributed double, with mean lambda
226 
227 ------------------------------------------------------------------------------*/
228 
229 double p_dExponential (double lambda,unsigned int *seed)
230 {
231  double u= p_dUniform( seed);
232 
233  while (u == (double)0)
234  u= p_dUniform( seed);
235 
236  return (-lambda*log(u));
237 }
238 /*----------------------------------------------------------------------
239 
240  The standard Gamma distribution has been implemented using
241  two functions:
242  <1> For integer orders
243  <2> For order in range (0..1)
244 
245  A gamma distribution of order X+Y is the same as that of the
246  sum of two independent gamma distributions of orders X and Y
247  respectivley. Hence, for general orders, we split the order
248  into it's integer and fractional parts, calculate their gamma
249  variates, and sum the variates obtained.
250 
251  Reference: The Art of Computer Programming, D.E. Knuth 2nd ed., 3.4.1e
252 ----------------------------------------------------------------------*/
253 
254 /*----------------------------------------------------------------------
255 
256  There are 2 methods for Gamma with integer A. The fisrt is to
257  get the sum of A independent exponential variates, is quick,
258  but grows linearly with A. The second method gets a closed
259  form and is more or less independent of A in speed. Therefore
260  for A>GAMMA_TURNING_POINT we should switch methods.
261 
262  GAMMA_TURNING_POINT is 7 for a Sun SS-10.
263 
264  Reference: D.E.D.E. Knuth,The art of computer programming,Vol 2
265  Chapter 3.4.1E algo A (p.129)
266 
267 ----------------------------------------------------------------------*/
268 
269 #define GAMMA_TURNING_POINT 7
270 
271 double p_dGammaInt (int a,unsigned int *seed)
272 {
273  if (a>GAMMA_TURNING_POINT)
274  {
275  double x,y,v,sqrt2a1,asub1,y2;
276 
277  sqrt2a1 = sqrt ((double)(2*a - 1));
278  asub1 = a-1;
279  while (1)
280  {
281  y = tan (3.1415926536*p_dUniform (seed));
282  y2 = y*sqrt2a1;
283  if ((x = y2 + asub1) <= 0)
284  continue;
285  v = p_dUniform (seed);
286  if (v <= (1+y*y)*exp((asub1)*log(x/asub1)-y2))
287  return (x);
288  }
289  }
290  else
291  {
292  double ans=1;
293  while (a--)
294  ans *= p_dUniform (seed);
295  return (-log(ans));
296  }
297 }
298 
299 /*----------------------------------------------------------------------
300 
301  GammaGen() takes any order A>0, and generates a gamma variate
302  as described above.
303 
304  The function incorporated in GammaGen() for generating gamma
305  variates in range (0..1) was given by Fishman. It is taken
306  from:
307  W.J.Kennedy,Jr & James E. Gentle, Statistical Computing, pg214
308 
309 ----------------------------------------------------------------------*/
310 
311 double p_dGammaGen (double A, unsigned *seed)
312 {
313 if (A <= 5)
314  {
315  double b,p,inva,res=0;
316  double EXP= 2.71828182;
317  double a,f;
318  int g;
319 
320  f= floor( A);
321  a= A-f;
322  g= (int) f;
323 
324  if( a!=0)
325  {
326  inva= 1/a;
327  b= 1.0 + a/EXP;
328 
329  while(1)
330  {
331  p= b*p_dUniform(seed);
332  if( p<=1.0 )
333  {
334  res= pow(p, inva);
335  if( p_dUniform(seed) <= exp( -res)) break;
336  }
337  else
338  {
339  res= - log( (b-p)/a);
340  if( p_dUniform(seed) <= 1.0/pow(res, 1-a)) break;
341  }
342  }
343  }
344 
345  if ( g>0)
346  res += p_dGammaInt(g,seed);
347 
348  return(res );
349  }
350 else /*Use the same function as integer*/
351 {
352  double x,y,v,sqrt2a1,asub1,y2;
353 
354  sqrt2a1 = sqrt (2*A - 1);
355  asub1 = A-1;
356  while (1)
357  {
358  y = tan (3.1415926536*p_dUniform (seed));
359  y2 = y*sqrt2a1;
360  if ((x = y2 + asub1) <= 0)
361  continue;
362  v = p_dUniform (seed);
363  if (v <= (1+y*y)*exp((asub1)*log(x/asub1)-y2))
364  return (x);
365  }
366 }
367 }
368 /*------------------------------------------------------------------------------
369 
370  This file contains functions to generate geometric random numbers.
371  P(1) = p
372  P(2) = p(1-p)
373  P(i) = p(1-p)^(i-1)
374  Keep generating a uniform random number till it is lesser than p.
375  OR
376  Take ceil (log (U)/log (p)).Though this is slower,we will select this
377  if p < .20 because as p falls,the iterative algo tends to slow
378 down and
379  they are equally fast at p = GEOMETRIC_TURNING_POINT.
380  Note that when p is large,the first method is nearly twice as
381 fast as
382  the second but when it can be made arbitarily slow by reducing p.
383 
384  GEOMETRIC_TURNING_POINT is 0.38 for a Sun SS-10.
385 
386  ERROR HANDLING: if p>1 or p<=0 then we have an error
387 
388 ------------------------------------------------------------------------------*/
389 
390 #define GEOMETRIC_TURNING_POINT 0.38
391 
392 int p_iGeometric (double p, unsigned int *seed)
393 {
394  int i=1;
395 
396  if (p < GEOMETRIC_TURNING_POINT)
397  return ((int) ceil(log (p_dUniform(seed))/log (1-p)));
398 
399  while (p_dUniform(seed) > p)
400  i++;
401  return (i);
402 }
403 
404 /*------------------------------------------------------------------------------
405 
406  This file contains functions to generate normally distributed
407  variables. It uses the polar method developed by G.E.P. Box, M.E.
408  Muller and G. Marsaglia.
409 
410  Reference: D.E.D.E. Knuth, The Art of Computer Programming Ed. 2, 3.4.1C
411 
412 ------------------------------------------------------------------------------*/
413 
414 /*------------------------------------------------------------------------------
415  p_dNormal() returns a normally distributed double variable with
416 mean
417  mean, standard deviation sigma.
418 
419 
420 ------------------------------------------------------------------------------*/
421 
422 double p_dNormal (double mean, double sigma,unsigned int *seed)
423 {
424 
425 double u1, u2, v1, v2, s;
426 static char flag= 0;
427 static double x1, x2;
428 
429 if( flag )
430 {
431 flag=0;
432 return(x2*sigma + mean);
433 }
434 
435 do {
436 u1= p_dUniform(seed);
437 u2= p_dUniform(seed);
438 
439 v1= 2* u1 -1;
440 v2= 2* u2 -1;
441 s= v1*v1 + v2*v2;
442 } while( s>=1);
443 
444 s= sqrt( -2*log(s)/s);
445 x1= v1* s;
446 x2= v2* s;
447 
448 flag= 1;
449 return( x1*sigma + mean);
450 
451 }
452 
453 int p_iPoisson(double mu, unsigned* seed)
454 {
455 double total=1, till;
456 int count= -1 ;
457 
458 till= exp( -mu);
459 while( total > till)
460 {
461 total *= p_dUniform(seed);
462 count++;
463 }
464 
465 return(count);
466 }
467 
468 /*------------------------------------------------------------------------------
469 
470  This file contains functions for the generation of uniformly
471  distributed random numbbers. The linear congruential method is
472  used. Function used is
473 
474  X[n+1]= ( a * X[n] + c ) mod m
475 
476  with
477  a= 1664525 ,
478  c= 1 ,
479  m= 2^32 .
480 
481  Reference : D.E. Knuth, The Art of Computer Programming, Chapter 3.
482 
483 ------------------------------------------------------------------------------*/
484 
485 /*
486 
487 --
488 --------------------------------------------------------------------------
489 Mark Feldman \\_==_o Skydivers do
490 E-mail : u914097@student.canberra.edu.au \\ it at 120mph....
491 --------------------------------------------------------------------------
492 */
493 
double p_dUniform(unsigned *seed)
Definition: random-coll.c:126
unsigned p_iUniform(unsigned *seed)
Definition: random-coll.c:113
double p_dNormal(double mean, double sigma, unsigned int *seed)
Definition: random-coll.c:422
#define _Lo(X)
Definition: random-coll.c:93
double p_dBeta(double a, double b, unsigned int *seed)
Definition: random-coll.c:141
static unsigned Scal
Definition: random-coll.c:125
#define _Mult
Definition: random-coll.c:90
#define _Cons
Definition: random-coll.c:91
tuple a
Definition: go4init.py:12
int p_iPoisson(double mu, unsigned *seed)
Definition: random-coll.c:453
#define BETA_TURNING_POINT
Definition: random-coll.c:88
#define _Hi(X)
Definition: random-coll.c:94
double p_dGammaInt(int a, unsigned int *seed)
Definition: random-coll.c:271
double p_dGammaGen(double A, unsigned *seed)
Definition: random-coll.c:311
double p_dExponential(double lambda, unsigned int *seed)
Definition: random-coll.c:229
int p_iGeometric(double p, unsigned int *seed)
Definition: random-coll.c:392
#define GAMMA_TURNING_POINT
Definition: random-coll.c:269
int p_iBinomial(double p, int n, unsigned int *seed)
Definition: random-coll.c:194
#define GEOMETRIC_TURNING_POINT
Definition: random-coll.c:390