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

Last change on this file since 3810 was 3791, checked in by cmv, 15 years ago

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