GSI Object Oriented Online Offline (Go4)  GO4-6.1.4
 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 2627 2019-10-01 08:02:45Z 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 (a,b,seed)
142 double a,b;
143 unsigned int *seed;
144 {
145  if (a*b < BETA_TURNING_POINT)
146  {
147  double u1,u2,p1,p2,x;
148 
149  p1 = 1/a;
150  p2 = 1/b;
151  do
152  {
153  u1 = pow (p_dUniform(seed),p1);
154  u2 = pow (p_dUniform(seed),p2);
155  } while ((x=u1+u2) > 1);
156  return (u1/x);
157  }
158  else
159  {
160  double x,y;
161 
162  x = p_dGammaGen (a,seed);
163  y = p_dGammaGen (b,seed);
164  return (x/(x+y));
165  }
166 }
167 
168 /*------------------------------------------------------------------------------
169 
170  This file contains the binomial random number generator.
171  Unfortunately, here I actually have to simulate it for want of
172  a better method. Thus if a binomial is required with
173  mean .3 and n=50,fifty uniform variables are required!
174 
175  Reference: The Art of Computer Programming, D.E. Knuth 2nd ed., 3.4.1f
176 
177 ------------------------------------------------------------------------------*/
178 
179 /*------------------------------------------------------------------------------
180  p_iBionamial() returns bionomially distributed random numbers
181  with n trials of events with a probability p.
182 
183  Input:
184  n : the numer of trials
185  p : the probability of success in each trial
186  seed : a pointer to an unsigned integer. This is used to
187  store previous uniform variate calculated for
188  calculation of next uniform variate.
189  It must not be changed by calling program.
190 
191  Output: a bionomially distributed integer, with n trials of
192  events with probability p.
193 
194 ------------------------------------------------------------------------------*/
195 
196 int p_iBinomial (p,n,seed)
197 unsigned int *seed;
198 int n;
199 double p;
200 {
201  int i,ans=0;
202 
203  for (i=0;i<n;i++)
204  if (p_dUniform (seed) <= p)
205  ans++;
206 
207  return (ans);
208 }
209 /*------------------------------------------------------------------------------
210 
211  This file contains exponential random number function.The method
212  is to get a uniform random number U in (0,1] and return -ln(U).
213  Care is taken to aviod U=0;
214 
215  Reference: The Art of Computer Programming, D.E. Knuth 2nd ed., 3.4.1d
216 
217 ------------------------------------------------------------------------------*/
218 
219 /*------------------------------------------------------------------------------
220  p_dExponential() returns exponentially distributed random numbers
221  with mean lambda.
222 
223  Input:
224  lambda: the mean.
225  seed : a pointer to an unsigned integer. This is used to
226  store previous uniform variate calculated for
227  calculation of next uniform variate.
228  It must not be changed by calling program.
229 
230  Output: an exponentially distributed double, with mean lambda
231 
232 ------------------------------------------------------------------------------*/
233 
234 double p_dExponential (double lambda,unsigned int *seed)
235 {
236  double u= p_dUniform( seed);
237 
238  while (u == (double)0)
239  u= p_dUniform( seed);
240 
241  return (-lambda*log(u));
242 }
243 /*----------------------------------------------------------------------
244 
245  The standard Gamma distribution has been implemented using
246  two functions:
247  <1> For integer orders
248  <2> For order in range (0..1)
249 
250  A gamma distribution of order X+Y is the same as that of the
251  sum of two independent gamma distributions of orders X and Y
252  respectivley. Hence, for general orders, we split the order
253  into it's integer and fractional parts, calculate their gamma
254  variates, and sum the variates obtained.
255 
256  Reference: The Art of Computer Programming, D.E. Knuth 2nd ed., 3.4.1e
257 ----------------------------------------------------------------------*/
258 
259 /*----------------------------------------------------------------------
260 
261  There are 2 methods for Gamma with integer A. The fisrt is to
262  get the sum of A independent exponential variates, is quick,
263  but grows linearly with A. The second method gets a closed
264  form and is more or less independent of A in speed. Therefore
265  for A>GAMMA_TURNING_POINT we should switch methods.
266 
267  GAMMA_TURNING_POINT is 7 for a Sun SS-10.
268 
269  Reference: D.E.D.E. Knuth,The art of computer programming,Vol 2
270  Chapter 3.4.1E algo A (p.129)
271 
272 ----------------------------------------------------------------------*/
273 
274 #define GAMMA_TURNING_POINT 7
275 
276 double p_dGammaInt (int a,unsigned int *seed)
277 {
278  if (a>GAMMA_TURNING_POINT)
279  {
280  double x,y,v,sqrt2a1,asub1,y2;
281 
282  sqrt2a1 = sqrt ((double)(2*a - 1));
283  asub1 = a-1;
284  while (1)
285  {
286  y = tan (3.1415926536*p_dUniform (seed));
287  y2 = y*sqrt2a1;
288  if ((x = y2 + asub1) <= 0)
289  continue;
290  v = p_dUniform (seed);
291  if (v <= (1+y*y)*exp((asub1)*log(x/asub1)-y2))
292  return (x);
293  }
294  }
295  else
296  {
297  double ans=1;
298  while (a--)
299  ans *= p_dUniform (seed);
300  return (-log(ans));
301  }
302 }
303 
304 /*----------------------------------------------------------------------
305 
306  GammaGen() takes any order A>0, and generates a gamma variate
307  as described above.
308 
309  The function incorporated in GammaGen() for generating gamma
310  variates in range (0..1) was given by Fishman. It is taken
311  from:
312  W.J.Kennedy,Jr & James E. Gentle, Statistical Computing, pg214
313 
314 ----------------------------------------------------------------------*/
315 
316 double p_dGammaGen (A,seed)
317 double A;
318 unsigned *seed;
319 {
320 if (A <= 5)
321  {
322  double b,p,inva,res=0;
323  double EXP= 2.71828182;
324  double a,f;
325  int g;
326 
327  f= floor( A);
328  a= A-f;
329  g= (int) f;
330 
331  if( a!=0)
332  {
333  inva= 1/a;
334  b= 1.0 + a/EXP;
335 
336  while(1)
337  {
338  p= b*p_dUniform(seed);
339  if( p<=1.0 )
340  {
341  res= pow(p, inva);
342  if( p_dUniform(seed) <= exp( -res)) break;
343  }
344  else
345  {
346  res= - log( (b-p)/a);
347  if( p_dUniform(seed) <= 1.0/pow(res, 1-a)) break;
348  }
349  }
350  }
351 
352  if ( g>0)
353  res += p_dGammaInt(g,seed);
354 
355  return(res );
356  }
357 else /*Use the same function as integer*/
358 {
359  double x,y,v,sqrt2a1,asub1,y2;
360 
361  sqrt2a1 = sqrt (2*A - 1);
362  asub1 = A-1;
363  while (1)
364  {
365  y = tan (3.1415926536*p_dUniform (seed));
366  y2 = y*sqrt2a1;
367  if ((x = y2 + asub1) <= 0)
368  continue;
369  v = p_dUniform (seed);
370  if (v <= (1+y*y)*exp((asub1)*log(x/asub1)-y2))
371  return (x);
372  }
373 }
374 }
375 /*------------------------------------------------------------------------------
376 
377  This file contains functions to generate geometric random numbers.
378  P(1) = p
379  P(2) = p(1-p)
380  P(i) = p(1-p)^(i-1)
381  Keep generating a uniform random number till it is lesser than p.
382  OR
383  Take ceil (log (U)/log (p)).Though this is slower,we will select this
384  if p < .20 because as p falls,the iterative algo tends to slow
385 down and
386  they are equally fast at p = GEOMETRIC_TURNING_POINT.
387  Note that when p is large,the first method is nearly twice as
388 fast as
389  the second but when it can be made arbitarily slow by reducing p.
390 
391  GEOMETRIC_TURNING_POINT is 0.38 for a Sun SS-10.
392 
393  ERROR HANDLING: if p>1 or p<=0 then we have an error
394 
395 ------------------------------------------------------------------------------*/
396 
397 #define GEOMETRIC_TURNING_POINT 0.38
398 
399 int p_iGeometric (double p, unsigned int *seed)
400 {
401  int i=1;
402 
403  if (p < GEOMETRIC_TURNING_POINT)
404  return ((int) ceil(log (p_dUniform(seed))/log (1-p)));
405 
406  while (p_dUniform(seed) > p)
407  i++;
408  return (i);
409 }
410 
411 /*------------------------------------------------------------------------------
412 
413  This file contains functions to generate normally distributed
414  variables. It uses the polar method developed by G.E.P. Box, M.E.
415  Muller and G. Marsaglia.
416 
417  Reference: D.E.D.E. Knuth, The Art of Computer Programming Ed. 2, 3.4.1C
418 
419 ------------------------------------------------------------------------------*/
420 
421 /*------------------------------------------------------------------------------
422  p_dNormal() returns a normally distributed double variable with
423 mean
424  mean, standard deviation sigma.
425 
426 
427 ------------------------------------------------------------------------------*/
428 
429 double p_dNormal (double mean, double sigma,unsigned int *seed)
430 {
431 
432 double u1, u2, v1, v2, s;
433 static char flag= 0;
434 static double x1, x2;
435 
436 if( flag )
437 {
438 flag=0;
439 return(x2*sigma + mean);
440 }
441 
442 do {
443 u1= p_dUniform(seed);
444 u2= p_dUniform(seed);
445 
446 v1= 2* u1 -1;
447 v2= 2* u2 -1;
448 s= v1*v1 + v2*v2;
449 } while( s>=1);
450 
451 s= sqrt( -2*log(s)/s);
452 x1= v1* s;
453 x2= v2* s;
454 
455 flag= 1;
456 return( x1*sigma + mean);
457 
458 }
459 
460 int p_iPoisson(double mu, unsigned* seed)
461 {
462 double total=1, till;
463 int count= -1 ;
464 
465 till= exp( -mu);
466 while( total > till)
467 {
468 total *= p_dUniform(seed);
469 count++;
470 }
471 
472 return(count);
473 }
474 
475 /*------------------------------------------------------------------------------
476 
477  This file contains functions for the generation of uniformly
478  distributed random numbbers. The linear congruential method is
479  used. Function used is
480 
481  X[n+1]= ( a * X[n] + c ) mod m
482 
483  with
484  a= 1664525 ,
485  c= 1 ,
486  m= 2^32 .
487 
488  Reference : D.E. Knuth, The Art of Computer Programming, Chapter 3.
489 
490 ------------------------------------------------------------------------------*/
491 
492 /*
493 
494 --
495 --------------------------------------------------------------------------
496 Mark Feldman \\_==_o Skydivers do
497 E-mail : u914097@student.canberra.edu.au \\ it at 120mph....
498 --------------------------------------------------------------------------
499 */
500 
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:429
#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:460
#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:276
double p_dGammaGen(double A, unsigned *seed)
Definition: random-coll.c:316
double p_dExponential(double lambda, unsigned int *seed)
Definition: random-coll.c:234
int p_iGeometric(double p, unsigned int *seed)
Definition: random-coll.c:399
#define GAMMA_TURNING_POINT
Definition: random-coll.c:274
int p_iBinomial(double p, int n, unsigned int *seed)
Definition: random-coll.c:196
#define GEOMETRIC_TURNING_POINT
Definition: random-coll.c:397