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

Last change on this file since 3603 was 3602, checked in by cmv, 16 years ago

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