source: Sophya/trunk/SophyaLib/BaseTools/randinterf.cc@ 3612

Last change on this file since 3612 was 3611, checked in by ansari, 16 years ago

Petite erreur de frappe corrigee suite compil OSF - Reza 30/04/2009
File size: 26.4 KB
RevLine 
[3602]1#include "sopnamsp.h"
2#include "machdefs.h"
3#include <math.h>
4#include <stdlib.h>
5#include <sys/time.h>
6#include <time.h>
7#include <iostream>
8#include "pexceptions.h"
9
10#include "randinterf.h"
11
12namespace SOPHYA {
13
14//-------------------------------------------------------------------------------
15// ------ Definition d'interface des classes de generateurs de nombres aleatoires
[3604]16/*!
17 \class RandomGeneratorInterface
18 \ingroup BaseTools
19 \brief Base class for random number generators
[3602]20
[3604]21 This class defines the interface for random number generator classes and
22 implements the generation of some specific distributions (Gaussian, Poisson ...)
23 through generation of random number with a flat distribution in the range [0,1[.
24
25 The sub classes inheriting from this class should implement the Next() method.
26
27 This base class manages also a global instance of a default generator.
28
29 \sa frand01 drand01 frandpm1 drandpm1
30 \sa GauRnd PoissRand
31
32*/
33
34
35RandomGeneratorInterface* RandomGeneratorInterface::gl_rndgen_p = NULL;
36
37/*!
38 \brief: static method to set or change the intance of the global Random Generator object
39
40 This method should be called during initialization, before any call to global
41 functions for random number generation. The rgp object should be created using new.
42*/
43void RandomGeneratorInterface::SetGlobalRandGenP(RandomGeneratorInterface* rgp)
44{
45 if (rgp == NULL) return;
46 if (gl_rndgen_p) delete gl_rndgen_p;
47 gl_rndgen_p = rgp;
48 return;
49}
50
[3602]51RandomGeneratorInterface::RandomGeneratorInterface()
52{
[3604]53 SelectGaussianAlgo();
54 SelectPoissonAlgo();
55 SelectExponentialAlgo();
[3602]56}
57
58
59RandomGeneratorInterface::~RandomGeneratorInterface(void)
60{
61 // rien a faire
62}
63
64
65/////////////////////////////////////////////////////////////////////////
66/////////////////////////////////////////////////////////////////////////
67/////////////////////////////////////////////////////////////////////////
68
69r_8 RandomGeneratorInterface::Next()
70{
71 printf("RandomGeneratorInterface::Next(): undefined code !!!\n");
72 throw MathExc("RandomGeneratorInterface::Next(): undefined code !!!");
73}
74
75/////////////////////////////////////////////////////////////////////////
76/////////////////////////////////////////////////////////////////////////
77/////////////////////////////////////////////////////////////////////////
78void RandomGeneratorInterface::GenerateSeedVector(int nseed,vector<uint_2>& seed,int lp)
79// renvoie un vecteur de nseed+2 entiers 32 bits
80// [0 - 2] = codage sur 48 bits du nombre (melange) de microsec depuis l'origine
81// [3 -> 3+ngene-1] = entiers aleatoires (poor man generator)
82//
83// L'initialiseur est donne par un codage du nombre de millisecondes
84// ecoulees depuis le 0 heure le 1er Janvier 1970 UTC (cf gettimeofday).
85// Seuls les 48 bits de poids faible sont retenus.
86// Un melange des bits est ensuite effectue pour que les 3 nombres
87// (unsigned short) d'initialisation ne soient pas trop semblables.
88// Le nombre le plus grand que l'on peut mettre
89// dans un entier unsigned de N bits est: 2^N-1
90// 48 bits -> 2^48-1 = 281474976710655 musec = 3257.8j = 8.9y
91// -> meme initialisation tous les 8.9 ans a 1 microsec pres !
92{
93 if(lp>0) cout<<"RandomGeneratorInterface::GenerateSeedVector: nseed="<<nseed<<endl;
94
95 // ---
96 // --- les deux premiers mots remplis avec le temps
97 // ---
98 // On recupere le temps ecoule depuis l'origine code en sec+musec
99 struct timeval now;
100 gettimeofday(&now,0);
101 // Calcul du temps ecoule depuis l'origine en microsecondes
102 uint_8 tmicro70 = (uint_8)now.tv_sec*(uint_8)1000000 + (uint_8)now.tv_usec;
103 if(lp>1) cout<<"."<<now.tv_sec<<" sec + "<<now.tv_usec<<" musec = "<<tmicro70<<" musec"<<endl;
104 // Remplissage du tableau de 48 bits
105 uint_2 b[48]; uint_8 tdum = tmicro70;
106 for(int ip=0;ip<48;ip++) {b[ip] = tdum&1; tdum = (tdum>>1);}
107 if(lp>2) {
108 cout<<"..b= ";
109 for(int ip=47;ip>=0;ip--) {cout<<b[ip]; if(ip%32==0 || ip%16==0) cout<<" ";}
110 cout<<endl;
111 }
112 // Melange des bits qui varient vite (poids faible, microsec)
113 // avec ceux variant lentement (poids fort, sec)
114 for(int ip=0;ip<16;ip++) {
115 if(ip%3==1) swap(b[ip],b[32+ip]);
116 else if(ip%3==2) swap(b[ip],b[16-ip]);
117 }
118 if(lp>2) {
119 cout<<"..b= ";
120 for(int ip=47;ip>=0;ip--) {cout<<b[ip]; if(ip%32==0 || ip%16==0) cout<<" ";}
121 cout<<endl;
122 }
123 // Remplissage
124 seed.resize(0);
125 for(int i=0;i<3;i++) {
126 seed.push_back(0);
127 uint_2 w = 1;
128 for(int ip=0;ip<16;ip++) {seed[i] += w*b[i*16+ip]; w *= 2;}
129 }
130 if(lp>0) cout<<"seed(time): "<<seed[0]<<" "<<seed[1]<<" "<<seed[2]<<endl;
131
132 // ---
133 // --- generation des nombres aleatoires complementaires (poor man generator)
134 // ---
135 //----------------------------------------------------------------------------//
136 // Ran088: L'Ecuyer's 1996 three-component Tausworthe generator "taus88"
137 // Returns an integer random number uniformly distributed within [0,4294967295]
138 // The period length is approximately 2^88 (which is 3*10^26).
139 // This generator is very fast and passes all standard statistical tests.
140 // Reference:
141 // (1) P. L'Ecuyer, Maximally equidistributed combined Tausworthe generators,
142 // Mathematics of Computation, 65, 203-213 (1996), see Figure 4.
143 // (2) recommended in:
144 // P. L'Ecuyer, Random number generation, chapter 4 of the
145 // Handbook on Simulation, Ed. Jerry Banks, Wiley, 1997.
146 //----------------------------------------------------------------------------//
147 if(nseed<=0) return;
148 // initialize seeds using the given seed value taking care of
149 // the requirements. The constants below are arbitrary otherwise
150 uint_4 seed0 = uint_4(tmicro70&0xFFFFFFFFULL);
151 if(lp>2) cout<<"seed0(time): "<<seed0<<endl;
152 uint_4 state_s1, state_s2, state_s3;
153 state_s1 = 1243598713U ^ seed0; if (state_s1 < 2) state_s1 = 1243598713U;
154 state_s2 = 3093459404U ^ seed0; if (state_s2 < 8) state_s2 = 3093459404U;
155 state_s3 = 1821928721U ^ seed0; if (state_s3 < 16) state_s3 = 1821928721U;
156 int nfill = 0, ico=0;
157 while(nfill<nseed) {
158 uint_4 s1 = state_s1, s2 = state_s2, s3 = state_s3;
159 // generate a random 32 bit number
160 s1 = ((s1 & -2) << 12) ^ (((s1 << 13) ^ s1) >> 19);
161 s2 = ((s2 & -8) << 4) ^ (((s2 << 2) ^ s2) >> 25);
162 s3 = ((s3 & -16) << 17) ^ (((s3 << 3) ^ s3) >> 11);
163 state_s1 = s1; state_s2 = s2; state_s3 = s3;
164 // le nombre aleatoire sur 32 bits est: s1^s2^s3
165 if(ico<15) {ico++; continue;} // des tirages blancs
166 uint_2 s = uint_2( (s1^s2^s3)&0xFFFFU );
167 seed.push_back(s);
168 if(lp>0) cout<<"seed(t88): "<<seed[3+nfill]<<endl;
169 nfill++;
170 }
171
172}
173
174/////////////////////////////////////////////////////////////////////////
175/////////////////////////////////////////////////////////////////////////
176/////////////////////////////////////////////////////////////////////////
177
[3604]178r_8 RandomGeneratorInterface::Gaussian()
179{
180 switch (usegaussian_) {
181 case C_Gaussian_BoxMuller :
182 return GaussianBoxMuller();
183 break;
184 case C_Gaussian_RandLibSNorm :
185 return GaussianSNorm();
186 break;
187 case C_Gaussian_PolarBoxMuller :
188 return GaussianPolarBoxMuller();
189 break;
190 case C_Gaussian_RatioUnif :
191 return GaussianRatioUnif();
192 break;
193 case C_Gaussian_LevaRatioUnif :
194 return GaussianLevaRatioUnif();
195 break;
196 default:
197 return GaussianBoxMuller();
198 break;
199 }
200}
201
[3602]202//--- Generation de nombre aleatoires suivant une distribution gaussienne
203r_8 RandomGeneratorInterface::GaussianBoxMuller()
204{
205 r_8 A=Next();
206 while (A==0.) A=Next();
207 return sqrt(-2.*log(A))*cos(2.*M_PI*Next());
208}
209
210//-------------------------------------------
211// Adapte de ranlib float snorm()
212// http://orion.math.iastate.edu/burkardt/c_src/ranlib/ranlib.c
213/*
214**********************************************************************
215 (STANDARD-) N O R M A L DISTRIBUTION
216**********************************************************************
217
218 FOR DETAILS SEE:
219
220 AHRENS, J.H. AND DIETER, U.
221 EXTENSIONS OF FORSYTHE'S METHOD FOR RANDOM
222 SAMPLING FROM THE NORMAL DISTRIBUTION.
223 MATH. COMPUT., 27,124 (OCT. 1973), 927 - 937.
224
225 ALL STATEMENT NUMBERS CORRESPOND TO THE STEPS OF ALGORITHM 'FL'
226 (M=5) IN THE ABOVE PAPER (SLIGHTLY MODIFIED IMPLEMENTATION)
227
228 Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of
229 SUNIF. The argument IR thus goes away.
230
231**********************************************************************
232 THE DEFINITIONS OF THE CONSTANTS A(K), D(K), T(K) AND
233 H(K) ARE ACCORDING TO THE ABOVEMENTIONED ARTICLE
234*/
235static double a_snorm[32] = {
236 0.0,3.917609E-2,7.841241E-2,0.11777,0.1573107,0.1970991,0.2372021,0.2776904,
237 0.3186394,0.36013,0.4022501,0.4450965,0.4887764,0.5334097,0.5791322,
238 0.626099,0.6744898,0.7245144,0.7764218,0.8305109,0.8871466,0.9467818,
239 1.00999,1.077516,1.150349,1.229859,1.318011,1.417797,1.534121,1.67594,
240 1.862732,2.153875
241};
242static double d_snorm[31] = {
243 0.0,0.0,0.0,0.0,0.0,0.2636843,0.2425085,0.2255674,0.2116342,0.1999243,
244 0.1899108,0.1812252,0.1736014,0.1668419,0.1607967,0.1553497,0.1504094,
245 0.1459026,0.14177,0.1379632,0.1344418,0.1311722,0.128126,0.1252791,
246 0.1226109,0.1201036,0.1177417,0.1155119,0.1134023,0.1114027,0.1095039
247};
248static float t_snorm[31] = {
249 7.673828E-4,2.30687E-3,3.860618E-3,5.438454E-3,7.0507E-3,8.708396E-3,
250 1.042357E-2,1.220953E-2,1.408125E-2,1.605579E-2,1.81529E-2,2.039573E-2,
251 2.281177E-2,2.543407E-2,2.830296E-2,3.146822E-2,3.499233E-2,3.895483E-2,
252 4.345878E-2,4.864035E-2,5.468334E-2,6.184222E-2,7.047983E-2,8.113195E-2,
253 9.462444E-2,0.1123001,0.136498,0.1716886,0.2276241,0.330498,0.5847031
254};
255static float h_snorm[31] = {
256 3.920617E-2,3.932705E-2,3.951E-2,3.975703E-2,4.007093E-2,4.045533E-2,
257 4.091481E-2,4.145507E-2,4.208311E-2,4.280748E-2,4.363863E-2,4.458932E-2,
258 4.567523E-2,4.691571E-2,4.833487E-2,4.996298E-2,5.183859E-2,5.401138E-2,
259 5.654656E-2,5.95313E-2,6.308489E-2,6.737503E-2,7.264544E-2,7.926471E-2,
260 8.781922E-2,9.930398E-2,0.11556,0.1404344,0.1836142,0.2790016,0.7010474
261};
262r_8 RandomGeneratorInterface::GaussianSNorm()
263{
264long i;
265double snorm,u,s,ustar,aa,w,y,tt;
266 u = Next();
267 s = 0.0;
268 if(u > 0.5) s = 1.0;
269 u += (u-s);
270 u = 32.0*u;
271 i = (long) (u);
272 if(i == 32) i = 31;
273 if(i == 0) goto S100;
274/*
275 START CENTER
276*/
277 ustar = u-(double)i;
278 aa = *(a_snorm+i-1);
279S40:
280 if(ustar <= *(t_snorm+i-1)) goto S60;
281 w = (ustar-*(t_snorm+i-1))**(h_snorm+i-1);
282S50:
283/*
284 EXIT (BOTH CASES)
285*/
286 y = aa+w;
287 snorm = y;
288 if(s == 1.0) snorm = -y;
289 return snorm;
290S60:
291/*
292 CENTER CONTINUED
293*/
294 u = Next();
295 w = u*(*(a_snorm+i)-aa);
296 tt = (0.5*w+aa)*w;
297 goto S80;
298S70:
299 tt = u;
300 ustar = Next();
301S80:
302 if(ustar > tt) goto S50;
303 u = Next();
304 if(ustar >= u) goto S70;
305 ustar = Next();
306 goto S40;
307S100:
308/*
309 START TAIL
310*/
311 i = 6;
312 aa = *(a_snorm+31);
313 goto S120;
314S110:
315 aa += *(d_snorm+i-1);
316 i += 1;
317S120:
318 u += u;
319 if(u < 1.0) goto S110;
320 u -= 1.0;
321S140:
322 w = u**(d_snorm+i-1);
323 tt = (0.5*w+aa)*w;
324 goto S160;
325S150:
326 tt = u;
327S160:
328 ustar = Next();
329 if(ustar > tt) goto S50;
330 u = Next();
331 if(ustar >= u) goto S150;
332 u = Next();
333 goto S140;
334}
335
336r_8 RandomGeneratorInterface::GaussianPolarBoxMuller()
337{
338double x1,x2,w;
339do {
340 x1 = 2.0 * Next() - 1.0;
341 x2 = 2.0 * Next() - 1.0;
342 w = x1 * x1 + x2 * x2;
343 } while ( w >= 1.0 || w==0. );
344return x1 * sqrt(-2.0*log(w)/w);
345}
346
347static double s2se_RatioUnif=sqrt(2./M_E) , epm135_RatioUnif=exp(-1.35) , ep1q_RatioUnif=exp(1./4.);
348r_8 RandomGeneratorInterface::GaussianRatioUnif()
349{
350double u,v,x;
351while(true) {
352 do {u = Next();} while ( u == 0. );
353 v = (2.0*Next()-1.0)*s2se_RatioUnif;
354 x = v/u;
355 if(x*x <= 5.0-4.0*ep1q_RatioUnif*u) break;
356 if(x*x<4.0*epm135_RatioUnif/u+1.4)
357 if(v*v<-4.0*u*u*log(u)) break;
358}
359return x;
360}
361
362r_8 RandomGeneratorInterface::GaussianLevaRatioUnif()
363{
364double u,v,x,y,q;
365 do {
366 u = 1.-Next(); // in ]0,1]
367 v = Next()-0.5; // in [-0.5, 0.5[
368 v *= 1.7156;
369 x = u - 0.449871;
370 y = ((v<0)?-v:v) + 0.386595;
371 q = x*x + y*(0.19600*y - 0.25472*x);
372 } while( q>=0.27597 && (q>0.27846 || v*v>-4.0*u*u*log(u)) );
373 return v/u;
374}
375
376r_8 RandomGeneratorInterface::GaussianTail(double s)
377{
378 /* Returns a gaussian random variable larger than a
379 * This implementation does one-sided upper-tailed deviates.
380 */
381
382 if (s < 1)
383 {
384 /* For small s, use a direct rejection method. The limit s < 1
385 can be adjusted to optimise the overall efficiency */
386 double x;
387 do
388 {
389 x = Gaussian();
390 }
391 while (x < s);
392 return x;
393 }
394 else
395 {
396 /* Use the "supertail" deviates from the last two steps
397 * of Marsaglia's rectangle-wedge-tail method, as described
398 * in Knuth, v2, 3rd ed, pp 123-128. (See also exercise 11, p139,
399 * and the solution, p586.)
400 */
401 double u, v, x;
402 do
403 {
404 u = Next();
405 do
406 {
407 v = Next();
408 }
409 while (v == 0.0);
410 x = sqrt (s * s - 2 * log (v));
411 }
412 while (x * u > s);
413 return x;
414 }
415}
416
417/////////////////////////////////////////////////////////////////////////
418/////////////////////////////////////////////////////////////////////////
419/////////////////////////////////////////////////////////////////////////
420
[3604]421uint_8 RandomGeneratorInterface::Poisson(double mu, double mumax)
422{
423 switch (usepoisson_) {
424 case C_Poisson_Simple :
425 return PoissonSimple(mu,mumax);
426 break;
427 case C_Poisson_Ahrens :
428 return PoissonAhrens(mu);
429 break;
430 default:
431 return PoissonSimple(mu,mumax);
432 break;
433 }
434}
435
436
[3602]437//--- Generation de nombre aleatoires suivant une distribution de Poisson
438uint_8 RandomGeneratorInterface::PoissonSimple(double mu,double mumax)
439{
440 double pp,ppi,x;
441
442 if((mumax>0.)&&(mu>=mumax)) {
443 pp = sqrt(mu);
444 while( (x=pp*Gaussian()) < -mu );
445 return (uint_8)(mu+x+0.5);
446 }
447 else {
448 uint_8 n;
449 ppi = pp = exp(-mu);
450 x = Next();
451 n = 0;
452 while (x > ppi) {
453 n++;
454 pp = mu*pp/(double)n;
455 ppi += pp;
456 }
457 return n;
458 }
459 return 0; // pas necessaire ?
460}
461
462
463static double a0_poiahr = -0.5;
464static double a1_poiahr = 0.3333333;
465static double a2_poiahr = -0.2500068;
466static double a3_poiahr = 0.2000118;
467static double a4_poiahr = -0.1661269;
468static double a5_poiahr = 0.1421878;
469static double a6_poiahr = -0.1384794;
470static double a7_poiahr = 0.125006;
471static double fact_poiahr[10] = {
472 1.0,1.0,2.0,6.0,24.0,120.0,720.0,5040.0,40320.0,362880.0};
473uint_8 RandomGeneratorInterface::PoissonAhrens(double mu)
474/*
475**********************************************************************
476 long ignpoi(float mu)
477 GENerate POIsson random deviate
478 Function
479 Generates a single random deviate from a Poisson
480 distribution with mean AV.
481 Arguments
482 av --> The mean of the Poisson distribution from which
483 a random deviate is to be generated.
484 genexp <-- The random deviate.
485 Method
486 Renames KPOIS from TOMS as slightly modified by BWB to use RANF
487 instead of SUNIF.
488 For details see:
489 Ahrens, J.H. and Dieter, U.
490 Computer Generation of Poisson Deviates
491 From Modified Normal Distributions.
492 ACM Trans. Math. Software, 8, 2
493 (June 1982),163-179
494**********************************************************************
495**********************************************************************
496
497
498 P O I S S O N DISTRIBUTION
499
500
501**********************************************************************
502**********************************************************************
503
504 FOR DETAILS SEE:
505
506 AHRENS, J.H. AND DIETER, U.
507 COMPUTER GENERATION OF POISSON DEVIATES
508 FROM MODIFIED NORMAL DISTRIBUTIONS.
509 ACM TRANS. MATH. SOFTWARE, 8,2 (JUNE 1982), 163 - 179.
510
511 (SLIGHTLY MODIFIED VERSION OF THE PROGRAM IN THE ABOVE ARTICLE)
512
513**********************************************************************
514 INTEGER FUNCTION IGNPOI(IR,MU)
515 INPUT: IR=CURRENT STATE OF BASIC RANDOM NUMBER GENERATOR
516 MU=MEAN MU OF THE POISSON DISTRIBUTION
517 OUTPUT: IGNPOI=SAMPLE FROM THE POISSON-(MU)-DISTRIBUTION
518 MUPREV=PREVIOUS MU, MUOLD=MU AT LAST EXECUTION OF STEP P OR B.
519 TABLES: COEFFICIENTS A0-A7 FOR STEP F. FACTORIALS FACT
520 COEFFICIENTS A(K) - FOR PX = FK*V*V*SUM(A(K)*V**K)-DEL
521 SEPARATION OF CASES A AND B
522*/
523{
[3611]524uint_8 ignpoi,j,k,kflag,l,m;
[3602]525double b1,b2,c,c0,c1,c2,c3,d,del,difmuk,e,fk,fx,fy,g,omega,p,p0,px,py,q,s,
526 t,u,v,x,xx,pp[35];
527
528 if(mu < 10.0) goto S120;
529/*
530 C A S E A. (RECALCULATION OF S,D,L IF MU HAS CHANGED)
531*/
532 s = sqrt(mu);
533 d = 6.0*mu*mu;
534/*
535 THE POISSON PROBABILITIES PK EXCEED THE DISCRETE NORMAL
536 PROBABILITIES FK WHENEVER K >= M(MU). L=IFIX(MU-1.1484)
537 IS AN UPPER BOUND TO M(MU) FOR ALL MU >= 10 .
538*/
539 l = (uint_8) (mu-1.1484);
540/*
541 STEP N. NORMAL SAMPLE - SNORM(IR) FOR STANDARD NORMAL DEVIATE
542*/
543 g = mu+s*Gaussian();
544 if(g < 0.0) goto S20;
545 ignpoi = (uint_8) (g);
546/*
547 STEP I. IMMEDIATE ACCEPTANCE IF IGNPOI IS LARGE ENOUGH
548*/
549 if(ignpoi >= l) return ignpoi;
550/*
551 STEP S. SQUEEZE ACCEPTANCE - SUNIF(IR) FOR (0,1)-SAMPLE U
552*/
553 fk = (double)ignpoi;
554 difmuk = mu-fk;
555 u = Next();
556 if(d*u >= difmuk*difmuk*difmuk) return ignpoi;
557S20:
558/*
559 STEP P. PREPARATIONS FOR STEPS Q AND H.
560 (RECALCULATIONS OF PARAMETERS IF NECESSARY)
561 .3989423=(2*PI)**(-.5) .416667E-1=1./24. .1428571=1./7.
562 THE QUANTITIES B1, B2, C3, C2, C1, C0 ARE FOR THE HERMITE
563 APPROXIMATIONS TO THE DISCRETE NORMAL PROBABILITIES FK.
564 C=.1069/MU GUARANTEES MAJORIZATION BY THE 'HAT'-FUNCTION.
565*/
566 omega = 0.3989423/s;
567 b1 = 4.166667E-2/mu;
568 b2 = 0.3*b1*b1;
569 c3 = 0.1428571*b1*b2;
570 c2 = b2-15.0*c3;
571 c1 = b1-6.0*b2+45.0*c3;
572 c0 = 1.0-b1+3.0*b2-15.0*c3;
573 c = 0.1069/mu;
574 if(g < 0.0) goto S50;
575/*
576 'SUBROUTINE' F IS CALLED (KFLAG=0 FOR CORRECT RETURN)
577*/
578 kflag = 0;
579 goto S70;
580S40:
581/*
582 STEP Q. QUOTIENT ACCEPTANCE (RARE CASE)
583*/
584 if(fy-u*fy <= py*exp(px-fx)) return ignpoi;
585S50:
586/*
587 STEP E. EXPONENTIAL SAMPLE - SEXPO(IR) FOR STANDARD EXPONENTIAL
588 DEVIATE E AND SAMPLE T FROM THE LAPLACE 'HAT'
589 (IF T <= -.6744 THEN PK < FK FOR ALL MU >= 10.)
590*/
[3609]591 e = Exponential();
[3602]592 u = Next();
593 u += (u-1.0);
594 //t = 1.8+fsign(e,u);
595 t = 1.8 + (((u>0. && e<0.) || (u<0. && e>0.))?-e:e);
596 if(t <= -0.6744) goto S50;
597 ignpoi = (uint_8) (mu+s*t);
598 fk = (double)ignpoi;
599 difmuk = mu-fk;
600/*
601 'SUBROUTINE' F IS CALLED (KFLAG=1 FOR CORRECT RETURN)
602*/
603 kflag = 1;
604 goto S70;
605S60:
606/*
607 STEP H. HAT ACCEPTANCE (E IS REPEATED ON REJECTION)
608*/
609 if(c*fabs(u) > py*exp(px+e)-fy*exp(fx+e)) goto S50;
610 return ignpoi;
611S70:
612/*
613 STEP F. 'SUBROUTINE' F. CALCULATION OF PX,PY,FX,FY.
614 CASE IGNPOI .LT. 10 USES FACTORIALS FROM TABLE FACT
615*/
616 if(ignpoi >= 10) goto S80;
617 px = -mu;
618 py = pow(mu,(double)ignpoi)/ *(fact_poiahr+ignpoi);
619 goto S110;
620S80:
621/*
622 CASE IGNPOI .GE. 10 USES POLYNOMIAL APPROXIMATION
623 A0-A7 FOR ACCURACY WHEN ADVISABLE
624 .8333333E-1=1./12. .3989423=(2*PI)**(-.5)
625*/
626 del = 8.333333E-2/fk;
627 del -= (4.8*del*del*del);
628 v = difmuk/fk;
629 if(fabs(v) <= 0.25) goto S90;
630 px = fk*log(1.0+v)-difmuk-del;
631 goto S100;
632S90:
633 px = fk*v*v*(((((((a7_poiahr*v+a6_poiahr)*v+a5_poiahr)*v+a4_poiahr)*v+a3_poiahr)*v+a2_poiahr)*v+a1_poiahr)*v+a0_poiahr)-del;
634S100:
635 py = 0.3989423/sqrt(fk);
636S110:
637 x = (0.5-difmuk)/s;
638 xx = x*x;
639 fx = -0.5*xx;
640 fy = omega*(((c3*xx+c2)*xx+c1)*xx+c0);
641 if(kflag <= 0) goto S40;
642 goto S60;
643S120:
644/*
645 C A S E B. (START NEW TABLE AND CALCULATE P0 IF NECESSARY)
646*/
647// m = max(1L,(long) (mu));
648 m = (1ULL >= (uint_8)mu) ? 1ULL: (uint_8)mu;
649
650 l = 0;
651 p = exp(-mu);
652 q = p0 = p;
653S130:
654/*
655 STEP U. UNIFORM SAMPLE FOR INVERSION METHOD
656*/
657 u = Next();
658 ignpoi = 0;
659 if(u <= p0) return ignpoi;
660/*
661 STEP T. TABLE COMPARISON UNTIL THE END PP(L) OF THE
662 PP-TABLE OF CUMULATIVE POISSON PROBABILITIES
663 (0.458=PP(9) FOR MU=10)
664*/
665 if(l == 0) goto S150;
666 j = 1;
667//if(u > 0.458) j = min(l,m);
668 if(u > 0.458) j = ((l<=m)? l: m);
669 for(k=j; k<=l; k++) {
670 if(u <= *(pp+k-1)) goto S180;
671 }
672 if(l == 35) goto S130;
673S150:
674/*
675 STEP C. CREATION OF NEW POISSON PROBABILITIES P
676 AND THEIR CUMULATIVES Q=PP(K)
677*/
678 l += 1;
679 for(k=l; k<=35; k++) {
680 p = p*mu/(double)k;
681 q += p;
682 *(pp+k-1) = q;
683 if(u <= q) goto S170;
684 }
685 l = 35;
686 goto S130;
687S170:
688 l = k;
689S180:
690 ignpoi = k;
691 return ignpoi;
692}
693
694/////////////////////////////////////////////////////////////////////////
695/////////////////////////////////////////////////////////////////////////
696/////////////////////////////////////////////////////////////////////////
697
[3604]698r_8 RandomGeneratorInterface::Exponential()
699{
700 switch (useexpo_) {
701 case C_Exponential_Simple :
702 return ExpoSimple();
703 break;
704 case C_Exponential_Ahrens :
705 return ExpoAhrens();
706 break;
707 default:
708 return ExpoSimple();
709 break;
710 }
711}
712
[3602]713r_8 RandomGeneratorInterface::ExpoSimple(void)
714{
715 return -log(1.-Next());
716}
717
718
719static double q_expo[8] = {
720 0.6931472,0.9333737,0.9888778,0.9984959,0.9998293,0.9999833,0.9999986,1.0};
721r_8 RandomGeneratorInterface::ExpoAhrens(void)
722/*
723**********************************************************************
724**********************************************************************
725 (STANDARD-) E X P O N E N T I A L DISTRIBUTION
726**********************************************************************
727**********************************************************************
728
729 FOR DETAILS SEE:
730
731 AHRENS, J.H. AND DIETER, U.
732 COMPUTER METHODS FOR SAMPLING FROM THE
733 EXPONENTIAL AND NORMAL DISTRIBUTIONS.
734 COMM. ACM, 15,10 (OCT. 1972), 873 - 882.
735
736 ALL STATEMENT NUMBERS CORRESPOND TO THE STEPS OF ALGORITHM
737 'SA' IN THE ABOVE PAPER (SLIGHTLY MODIFIED IMPLEMENTATION)
738
739 Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of
740 SUNIF. The argument IR thus goes away.
741
742**********************************************************************
743 Q(N) = SUM(ALOG(2.0)**K/K!) K=1,..,N , THE HIGHEST N
744 (HERE 8) IS DETERMINED BY Q(N)=1.0 WITHIN STANDARD PRECISION
745*/
746{
747long i;
748double sexpo,a,u,ustar,umin;
749double *q1 = q_expo;
750 a = 0.0;
751 while((u=Next())==0.);
752 goto S30;
753S20:
754 a += *q1;
755S30:
756 u += u;
757 if(u <= 1.0) goto S20;
758 u -= 1.0;
759 if(u > *q1) goto S60;
760 sexpo = a+u;
761 return sexpo;
762S60:
763 i = 1;
764 ustar = Next();
765 umin = ustar;
766S70:
767 ustar = Next();
768 if(ustar < umin) umin = ustar;
769 i += 1;
770 if(u > *(q_expo+i-1)) goto S70;
771 sexpo = a+umin**q1;
772 return sexpo;
773}
774
775
776} /* namespace SOPHYA */
777
778
779
780/////////////////////////////////////////////////////////////////
781/*
782**** Remarques sur complex< r_8 > ComplexGaussRan(double sig) ****
783
784--- variables gaussiennes x,y independantes
785x gaussien: pdf f(x) = 1/(sqrt(2Pi) Sx) exp(-(x-Mx)^2/(2 Sx^2))
786y gaussien: pdf f(y) = 1/(sqrt(2Pi) Sy) exp(-(y-My)^2/(2 Sy^2))
787x,y independants --> pdf f(x,y) = f(x) f(y)
788On a:
789 <x> = Integrate[x*f(x)] = Mx
790 <x^2> = Integrate[x^2*f(x)] = Mx^2 + Sx^2
791
792--- On cherche la pdf g(r,t) du module et de la phase
793 x = r cos(t) , y = r sin(t)
794 r=sqrt(x^2+y^2 , t=atan2(y,x)
795 (r,t) --> (x,y): le Jacobien = r
796
797 g(r,t) = r f(x,y) = r f(x) f(y)
798 = r/(2Pi Sx Sy) exp(-(x-Mx)^2/(2 Sx^2)) exp(-(y-My)^2/(2 Sy^2))
799
800- Le cas general est complique
801 (cf D.Pelat cours DEA "bruits et signaux" section 4.5)
802
803- Cas ou "Mx = My = 0" et "Sx = Sy = S"
804 c'est la pdf du module et de la phase d'un nombre complexe
805 dont les parties reelles et imaginaires sont independantes
806 et sont distribuees selon des gaussiennes de variance S^2
807 g(r,t) = r/(2Pi S^2) exp(-r^2/(2 S^2))
808 La distribution de "r" est donc:
809 g(r) = Integrate[g(r,t),{t,0,2Pi}]
810 = r/S^2 exp(-r^2/(2 S^2))
811 La distribution de "t" est donc:
812 g(t) = Integrate[g(r,t),{r,0,Infinity}]
813 = 1 / 2Pi (distribution uniforme sur [0,2Pi[)
814 Les variables aleatoires r,t sont independantes:
815 g(r,t) = g(r) g(t)
816On a:
817 <r> = Integrate[r*g(r)] = sqrt(PI/2)*S
818 <r^2> = Integrate[r^2*g(r)] = 2*S^2
819 <r^3> = Integrate[r^3*g(r)] = 3*sqrt(PI/2)*S^3
820 <r^4> = Integrate[r^4*g(r)] = 8*S^4
821
822- Attention:
823La variable complexe "c = x+iy = r*exp(i*t)" ainsi definie verifie:
824 <|c|^2> = <c c*> = <x^2+y^2> = <r^2> = 2 S^2
825Si on veut generer une variable complexe gaussienne telle que
826 <c c*> = s^2 alors il faut prendre S = s/sqrt(2) comme argument
827
828*/
Note: See TracBrowser for help on using the repository browser.