GSI Object Oriented Online Offline (Go4) GO4-6.4.0
Loading...
Searching...
No Matches
random-coll.c
Go to the documentation of this file.
1// $Id$
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
113unsigned p_iUniform(unsigned *seed)
114{
115unsigned lo, hi;
116
117lo= _Lo(_Lo(*seed) * _Lo(_Mult) + _Cons);
118hi= _Lo(_Hi(*seed) * _Lo(_Mult)) + _Lo(_Hi(_Mult) * _Lo(*seed)) +
119 _Hi((_Lo(*seed) * _Lo(_Mult)) + _Cons);
120
121*seed= (hi<<16 | lo);
122return( *seed );
123}
124
125static unsigned Scal=0xFFFFFFFF;
126double p_dUniform(unsigned* seed)
127{
128unsigned lo, hi;
129
130lo= _Lo(_Lo(*seed) * _Lo(_Mult) + _Cons);
131hi= _Lo(_Hi(*seed) * _Lo(_Mult)) + _Lo(_Hi(_Mult) * _Lo(*seed)) +
132 _Hi((_Lo(*seed) * _Lo(_Mult)) + _Cons);
133
134*seed= (hi<<16 | lo);
135return( ((double)*seed)/Scal );
136}
137
138
139
140
141double 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
194int 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 avoid 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
229double 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
271double p_dGammaInt (int a,unsigned int *seed)
272{
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
311double p_dGammaGen (double A, unsigned *seed)
312{
313if (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 }
350else /*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
378down 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
381fast 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
392int p_iGeometric (double p, unsigned int *seed)
393{
394 int i=1;
395
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
416mean
417 mean, standard deviation sigma.
418
419
420------------------------------------------------------------------------------*/
421
422double p_dNormal (double mean, double sigma,unsigned int *seed)
423{
424
425double u1, u2, v1, v2, s;
426static char flag = 0;
427static double x1, x2;
428
429if( flag )
430{
431flag = 0;
432return(x2*sigma + mean);
433}
434
435do {
436u1= p_dUniform(seed);
437u2= p_dUniform(seed);
438
439v1= 2* u1 -1;
440v2= 2* u2 -1;
441s= v1*v1 + v2*v2;
442} while( s>=1);
443
444s= sqrt( -2*log(s)/s);
445x1= v1* s;
446x2= v2* s;
447
448flag= 1;
449return( x1*sigma + mean);
450
451}
452
453int p_iPoisson(double mu, unsigned* seed)
454{
455double total=1, till;
456int count= -1 ;
457
458till= exp( -mu);
459while( total > till)
460{
461total *= p_dUniform(seed);
462count++;
463}
464
465return(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--------------------------------------------------------------------------
489Mark Feldman \\_==_o Skydivers do
490E-mail : u914097@student.canberra.edu.au \\ it at 120mph....
491--------------------------------------------------------------------------
492*/
493
double p_dUniform(unsigned *seed)
double p_dNormal(double mean, double sigma, unsigned int *seed)
int p_iBinomial(double p, int n, unsigned int *seed)
#define BETA_TURNING_POINT
Definition random-coll.c:88
unsigned p_iUniform(unsigned *seed)
double p_dGammaInt(int a, unsigned int *seed)
#define _Mult
Definition random-coll.c:90
#define _Lo(X)
Definition random-coll.c:93
int p_iGeometric(double p, unsigned int *seed)
#define _Hi(X)
Definition random-coll.c:94
#define _Cons
Definition random-coll.c:91
#define GAMMA_TURNING_POINT
double p_dExponential(double lambda, unsigned int *seed)
int p_iPoisson(double mu, unsigned *seed)
static unsigned Scal
double p_dGammaGen(double A, unsigned *seed)
#define GEOMETRIC_TURNING_POINT
double p_dBeta(double a, double b, unsigned int *seed)