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

Last change on this file since 4018 was 4018, checked in by cmv, 14 years ago
  • introduction du tirage gaussien par la methode de la ziggurat
  • ajout argument slim par defaut dans GaussianTail pour proteger d'une reccursion infinie genere par le tirage gaussien ziggurat dans le cas ou GaussianTail applique la methode par rejection (sdev<slim)

cmv, 21/09/2011
File size: 38.5 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 <typeinfo>
10
11#include "pexceptions.h"
12
13#include "randinterf.h"
14
15namespace SOPHYA {
16
17//-------------------------------------------------------------------------------
18// ------ Definition d'interface des classes de generateurs de nombres aleatoires
19/*!
20 \class RandomGeneratorInterface
21 \ingroup BaseTools
22 \brief Base class for random number generators
23
24 This class defines the interface for random number generator classes and
25 implements the generation of some specific distributions (Gaussian, Poisson ...)
26 through generation of random number with a flat distribution in the range [0,1[.
27
28 The sub classes inheriting from this class should implement the Next() method.
29
30 This base class manages also a global instance of a default generator.
31
32 \sa frand01 drand01 frandpm1 drandpm1
33 \sa Gaussian Poisson
34
35*/
36
37
38RandomGeneratorInterface* RandomGeneratorInterface::gl_rndgen_p = NULL;
39
40/*!
41 \brief: static method to set or change the intance of the global Random Generator object
42
43 This method should be called during initialization, before any call to global
44 functions for random number generation. The rgp object should be created using new.
45*/
46void RandomGeneratorInterface::SetGlobalRandGenP(RandomGeneratorInterface* rgp)
47{
48 if (rgp == NULL) return;
49 if (gl_rndgen_p) delete gl_rndgen_p;
50 gl_rndgen_p = rgp;
51 return;
52}
53
54RandomGeneratorInterface::RandomGeneratorInterface()
55{
56 SelectGaussianAlgo();
57 SelectPoissonAlgo();
58 SelectExponentialAlgo();
59}
60
61
62RandomGeneratorInterface::~RandomGeneratorInterface(void)
63{
64 // rien a faire
65}
66
67void RandomGeneratorInterface::ShowRandom()
68{
69 cout << " RandomGeneratorInterface::ShowRandom() typeid(this)=" << typeid(*this).name() << " @ "
70 << hex << (unsigned long)(this) << dec << endl;
71 return;
72}
73
74/////////////////////////////////////////////////////////////////////////
75/////////////////////////////////////////////////////////////////////////
76/////////////////////////////////////////////////////////////////////////
77
78/*
79r_8 RandomGeneratorInterface::Next()
80{
81 printf("RandomGeneratorInterface::Next(): undefined code !!!\n");
82 throw MathExc("RandomGeneratorInterface::Next(): undefined code !!!");
83}
84*/
85
86/////////////////////////////////////////////////////////////////////////
87/////////////////////////////////////////////////////////////////////////
88/////////////////////////////////////////////////////////////////////////
89void RandomGeneratorInterface::GenerateSeedVector(int nseed,vector<uint_2>& seed,int lp)
90// renvoie un vecteur de nseed+2 entiers 32 bits
91// [0 - 2] = codage sur 48 bits du nombre (melange) de microsec depuis l'origine
92// [3 -> 3+ngene-1] = entiers aleatoires (poor man generator)
93//
94// L'initialiseur est donne par un codage du nombre de millisecondes
95// ecoulees depuis le 0 heure le 1er Janvier 1970 UTC (cf gettimeofday).
96// Seuls les 48 bits de poids faible sont retenus.
97// Un melange des bits est ensuite effectue pour que les 3 nombres
98// (unsigned short) d'initialisation ne soient pas trop semblables.
99// Le nombre le plus grand que l'on peut mettre
100// dans un entier unsigned de N bits est: 2^N-1
101// 48 bits -> 2^48-1 = 281474976710655 musec = 3257.8j = 8.9y
102// -> meme initialisation tous les 8.9 ans a 1 microsec pres !
103{
104 if(lp>0) cout<<"RandomGeneratorInterface::GenerateSeedVector: nseed="<<nseed<<endl;
105
106 // ---
107 // --- les deux premiers mots remplis avec le temps
108 // ---
109 // On recupere le temps ecoule depuis l'origine code en sec+musec
110 struct timeval now;
111 gettimeofday(&now,0);
112 // Calcul du temps ecoule depuis l'origine en microsecondes
113 uint_8 tmicro70 = (uint_8)now.tv_sec*(uint_8)1000000 + (uint_8)now.tv_usec;
114 if(lp>1) cout<<"."<<now.tv_sec<<" sec + "<<now.tv_usec<<" musec = "<<tmicro70<<" musec"<<endl;
115 // Remplissage du tableau de 48 bits
116 uint_2 b[48]; uint_8 tdum = tmicro70;
117 for(int ip=0;ip<48;ip++) {b[ip] = tdum&1; tdum = (tdum>>1);}
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 // Melange des bits qui varient vite (poids faible, microsec)
124 // avec ceux variant lentement (poids fort, sec)
125 for(int ip=0;ip<16;ip++) {
126 if(ip%3==1) swap(b[ip],b[32+ip]);
127 else if(ip%3==2) swap(b[ip],b[16-ip]);
128 }
129 if(lp>2) {
130 cout<<"..b= ";
131 for(int ip=47;ip>=0;ip--) {cout<<b[ip]; if(ip%32==0 || ip%16==0) cout<<" ";}
132 cout<<endl;
133 }
134 // Remplissage
135 seed.resize(0);
136 for(int i=0;i<3;i++) {
137 seed.push_back(0);
138 uint_2 w = 1;
139 for(int ip=0;ip<16;ip++) {seed[i] += w*b[i*16+ip]; w *= 2;}
140 }
141 if(lp>0) cout<<"seed(time): "<<seed[0]<<" "<<seed[1]<<" "<<seed[2]<<endl;
142
143 // ---
144 // --- generation des nombres aleatoires complementaires (poor man generator)
145 // ---
146 //----------------------------------------------------------------------------//
147 // Ran088: L'Ecuyer's 1996 three-component Tausworthe generator "taus88"
148 // Returns an integer random number uniformly distributed within [0,4294967295]
149 // The period length is approximately 2^88 (which is 3*10^26).
150 // This generator is very fast and passes all standard statistical tests.
151 // Reference:
152 // (1) P. L'Ecuyer, Maximally equidistributed combined Tausworthe generators,
153 // Mathematics of Computation, 65, 203-213 (1996), see Figure 4.
154 // (2) recommended in:
155 // P. L'Ecuyer, Random number generation, chapter 4 of the
156 // Handbook on Simulation, Ed. Jerry Banks, Wiley, 1997.
157 //----------------------------------------------------------------------------//
158 if(nseed<=0) return;
159 // initialize seeds using the given seed value taking care of
160 // the requirements. The constants below are arbitrary otherwise
161 uint_4 seed0 = uint_4(tmicro70&0xFFFFFFFFULL);
162 if(lp>2) cout<<"seed0(time): "<<seed0<<endl;
163 uint_4 state_s1, state_s2, state_s3;
164 state_s1 = 1243598713U ^ seed0; if (state_s1 < 2) state_s1 = 1243598713U;
165 state_s2 = 3093459404U ^ seed0; if (state_s2 < 8) state_s2 = 3093459404U;
166 state_s3 = 1821928721U ^ seed0; if (state_s3 < 16) state_s3 = 1821928721U;
167 int nfill = 0, ico=0;
168 while(nfill<nseed) {
169 uint_4 s1 = state_s1, s2 = state_s2, s3 = state_s3;
170 // generate a random 32 bit number
171 s1 = ((s1 & -2) << 12) ^ (((s1 << 13) ^ s1) >> 19);
172 s2 = ((s2 & -8) << 4) ^ (((s2 << 2) ^ s2) >> 25);
173 s3 = ((s3 & -16) << 17) ^ (((s3 << 3) ^ s3) >> 11);
174 state_s1 = s1; state_s2 = s2; state_s3 = s3;
175 // le nombre aleatoire sur 32 bits est: s1^s2^s3
176 if(ico<15) {ico++; continue;} // des tirages blancs
177 uint_2 s = uint_2( (s1^s2^s3)&0xFFFFU );
178 seed.push_back(s);
179 if(lp>0) cout<<"seed(t88): "<<seed[3+nfill]<<endl;
180 nfill++;
181 }
182
183}
184
185void RandomGeneratorInterface::AutoInit(int lp)
186{
187 printf("RandomGeneratorInterface::AutoInit(): undefined code !!!\n");
188 throw MathExc("RandomGeneratorInterface::AutoInit(): undefined code !!!");
189}
190
191/////////////////////////////////////////////////////////////////////////
192/////////////////////////////////////////////////////////////////////////
193/////////////////////////////////////////////////////////////////////////
194
195r_8 RandomGeneratorInterface::Gaussian()
196{
197 switch (usegaussian_) {
198 case C_Gaussian_BoxMuller :
199 return GaussianBoxMuller();
200 break;
201 case C_Gaussian_RandLibSNorm :
202 return GaussianSNorm();
203 break;
204 case C_Gaussian_PolarBoxMuller :
205 return GaussianPolarBoxMuller();
206 break;
207 case C_Gaussian_RatioUnif :
208 return GaussianRatioUnif();
209 break;
210 case C_Gaussian_LevaRatioUnif :
211 return GaussianLevaRatioUnif();
212 break;
213 case C_Gaussian_Ziggurat128 :
214 return GaussianZiggurat128();
215 break;
216 default:
217 return GaussianBoxMuller();
218 break;
219 }
220}
221
222//--- Generation de nombre aleatoires suivant une distribution gaussienne
223r_8 RandomGeneratorInterface::GaussianBoxMuller()
224{
225 r_8 A=Next();
226 while (A==0.) A=Next();
227 return sqrt(-2.*log(A))*cos(2.*M_PI*Next());
228}
229
230//-------------------------------------------
231// Adapte de ranlib float snorm()
232// http://orion.math.iastate.edu/burkardt/c_src/ranlib/ranlib.c
233/*
234**********************************************************************
235 (STANDARD-) N O R M A L DISTRIBUTION
236**********************************************************************
237
238 FOR DETAILS SEE:
239
240 AHRENS, J.H. AND DIETER, U.
241 EXTENSIONS OF FORSYTHE'S METHOD FOR RANDOM
242 SAMPLING FROM THE NORMAL DISTRIBUTION.
243 MATH. COMPUT., 27,124 (OCT. 1973), 927 - 937.
244
245 ALL STATEMENT NUMBERS CORRESPOND TO THE STEPS OF ALGORITHM 'FL'
246 (M=5) IN THE ABOVE PAPER (SLIGHTLY MODIFIED IMPLEMENTATION)
247
248 Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of
249 SUNIF. The argument IR thus goes away.
250
251**********************************************************************
252 THE DEFINITIONS OF THE CONSTANTS A(K), D(K), T(K) AND
253 H(K) ARE ACCORDING TO THE ABOVEMENTIONED ARTICLE
254*/
255static double a_snorm[32] = {
256 0.0,3.917609E-2,7.841241E-2,0.11777,0.1573107,0.1970991,0.2372021,0.2776904,
257 0.3186394,0.36013,0.4022501,0.4450965,0.4887764,0.5334097,0.5791322,
258 0.626099,0.6744898,0.7245144,0.7764218,0.8305109,0.8871466,0.9467818,
259 1.00999,1.077516,1.150349,1.229859,1.318011,1.417797,1.534121,1.67594,
260 1.862732,2.153875
261};
262static double d_snorm[31] = {
263 0.0,0.0,0.0,0.0,0.0,0.2636843,0.2425085,0.2255674,0.2116342,0.1999243,
264 0.1899108,0.1812252,0.1736014,0.1668419,0.1607967,0.1553497,0.1504094,
265 0.1459026,0.14177,0.1379632,0.1344418,0.1311722,0.128126,0.1252791,
266 0.1226109,0.1201036,0.1177417,0.1155119,0.1134023,0.1114027,0.1095039
267};
268static float t_snorm[31] = {
269 7.673828E-4,2.30687E-3,3.860618E-3,5.438454E-3,7.0507E-3,8.708396E-3,
270 1.042357E-2,1.220953E-2,1.408125E-2,1.605579E-2,1.81529E-2,2.039573E-2,
271 2.281177E-2,2.543407E-2,2.830296E-2,3.146822E-2,3.499233E-2,3.895483E-2,
272 4.345878E-2,4.864035E-2,5.468334E-2,6.184222E-2,7.047983E-2,8.113195E-2,
273 9.462444E-2,0.1123001,0.136498,0.1716886,0.2276241,0.330498,0.5847031
274};
275static float h_snorm[31] = {
276 3.920617E-2,3.932705E-2,3.951E-2,3.975703E-2,4.007093E-2,4.045533E-2,
277 4.091481E-2,4.145507E-2,4.208311E-2,4.280748E-2,4.363863E-2,4.458932E-2,
278 4.567523E-2,4.691571E-2,4.833487E-2,4.996298E-2,5.183859E-2,5.401138E-2,
279 5.654656E-2,5.95313E-2,6.308489E-2,6.737503E-2,7.264544E-2,7.926471E-2,
280 8.781922E-2,9.930398E-2,0.11556,0.1404344,0.1836142,0.2790016,0.7010474
281};
282r_8 RandomGeneratorInterface::GaussianSNorm()
283{
284long i;
285double snorm,u,s,ustar,aa,w,y,tt;
286 u = Next();
287 s = 0.0;
288 if(u > 0.5) s = 1.0;
289 u += (u-s);
290 u = 32.0*u;
291 i = (long) (u);
292 if(i == 32) i = 31;
293 if(i == 0) goto S100;
294/*
295 START CENTER
296*/
297 ustar = u-(double)i;
298 aa = *(a_snorm+i-1);
299S40:
300 if(ustar <= *(t_snorm+i-1)) goto S60;
301 w = (ustar-*(t_snorm+i-1))**(h_snorm+i-1);
302S50:
303/*
304 EXIT (BOTH CASES)
305*/
306 y = aa+w;
307 snorm = y;
308 if(s == 1.0) snorm = -y;
309 return snorm;
310S60:
311/*
312 CENTER CONTINUED
313*/
314 u = Next();
315 w = u*(*(a_snorm+i)-aa);
316 tt = (0.5*w+aa)*w;
317 goto S80;
318S70:
319 tt = u;
320 ustar = Next();
321S80:
322 if(ustar > tt) goto S50;
323 u = Next();
324 if(ustar >= u) goto S70;
325 ustar = Next();
326 goto S40;
327S100:
328/*
329 START TAIL
330*/
331 i = 6;
332 aa = *(a_snorm+31);
333 goto S120;
334S110:
335 aa += *(d_snorm+i-1);
336 i += 1;
337S120:
338 u += u;
339 if(u < 1.0) goto S110;
340 u -= 1.0;
341S140:
342 w = u**(d_snorm+i-1);
343 tt = (0.5*w+aa)*w;
344 goto S160;
345S150:
346 tt = u;
347S160:
348 ustar = Next();
349 if(ustar > tt) goto S50;
350 u = Next();
351 if(ustar >= u) goto S150;
352 u = Next();
353 goto S140;
354}
355
356r_8 RandomGeneratorInterface::GaussianPolarBoxMuller()
357{
358double x1,x2,w;
359do {
360 x1 = 2.0 * Next() - 1.0;
361 x2 = 2.0 * Next() - 1.0;
362 w = x1 * x1 + x2 * x2;
363 } while ( w >= 1.0 || w==0. );
364return x1 * sqrt(-2.0*log(w)/w);
365}
366
367static double s2se_RatioUnif=sqrt(2./M_E) , epm135_RatioUnif=exp(-1.35) , ep1q_RatioUnif=exp(1./4.);
368r_8 RandomGeneratorInterface::GaussianRatioUnif()
369{
370double u,v,x;
371while(true) {
372 do {u = Next();} while ( u == 0. );
373 v = (2.0*Next()-1.0)*s2se_RatioUnif;
374 x = v/u;
375 if(x*x <= 5.0-4.0*ep1q_RatioUnif*u) break;
376 if(x*x<4.0*epm135_RatioUnif/u+1.4)
377 if(v*v<-4.0*u*u*log(u)) break;
378}
379return x;
380}
381
382r_8 RandomGeneratorInterface::GaussianLevaRatioUnif()
383{
384double u,v,x,y,q;
385 do {
386 u = 1.-Next(); // in ]0,1]
387 v = Next()-0.5; // in [-0.5, 0.5[
388 v *= 1.7156;
389 x = u - 0.449871;
390 y = ((v<0)?-v:v) + 0.386595;
391 q = x*x + y*(0.19600*y - 0.25472*x);
392 } while( q>=0.27597 && (q>0.27846 || v*v>-4.0*u*u*log(u)) );
393 return v/u;
394}
395
396//-------------------------------------------------------------------------------
397// Definition des tableaux pour tirage Gaussien methode Ziggurat N=128 bandes
398// G.Marsaglia and W.W.Tsang "The ziggurat method for generating random variables
399// D.Thomas et al. ACM Computing Surveys, Vol 39, No 4, Article 11, October 2007
400// http://en.wikipedia.org/wiki/Ziggurat_algorithm (tres bien explique)
401// Calcul des tableaux avec programme "cmvziggurat.cc"
402//-------------------------------------------------------------------------------
403const int N_ZIGGURRAT_GAUSS128 = 128;
404static const double X_ZIGGURRAT_GAUSS128[N_ZIGGURRAT_GAUSS128+1] = {
4050.0000000000000000, 0.2723208647046734, 0.3628714310284243, 0.4265479863033096, 0.4774378372537916,
4060.5206560387251481, 0.5586921783755209, 0.5929629424419807, 0.6243585973090908, 0.6534786387150446,
4070.6807479186459064, 0.7064796113136101, 0.7309119106218833, 0.7542306644345121, 0.7765839878761502,
4080.7980920606262765, 0.8188539066833194, 0.8389522142812090, 0.8584568431780525, 0.8774274290977171,
4090.8959153525662399, 0.9139652510088031, 0.9316161966013551, 0.9489026254979132, 0.9658550793881319,
4100.9825008035027615, 0.9988642334806447, 1.0149673952393006, 1.0308302360564565, 1.0464709007525812,
4111.0619059636836206, 1.0771506248819389, 1.0922188768965548, 1.1071236475235364, 1.1218769225722551,
4121.1364898520030764, 1.1509728421389769, 1.1653356361550478, 1.1795873846544616, 1.1937367078237728,
4131.2077917504067583, 1.2217602305309634, 1.2356494832544818, 1.2494664995643345, 1.2632179614460288,
4141.2769102735517004, 1.2905495919178738, 1.3041418501204223, 1.3176927832013436, 1.3312079496576772,
4151.3446927517457137, 1.3581524543224235, 1.3715922024197329, 1.3850170377251492, 1.3984319141236070,
4161.4118417124397606, 1.4252512545068619, 1.4386653166774619, 1.4520886428822168, 1.4655259573357950,
4171.4789819769830983, 1.4924614237746157, 1.5059690368565506, 1.5195095847593711, 1.5330878776675565,
4181.5467087798535037, 1.5603772223598409, 1.5740982160167500, 1.5878768648844011, 1.6017183802152775,
4191.6156280950371333, 1.6296114794646788, 1.6436741568569830, 1.6578219209482079, 1.6720607540918526,
4201.6863968467734867, 1.7008366185643014, 1.7153867407081171, 1.7300541605582440, 1.7448461281083769,
4211.7597702248942324, 1.7748343955807697, 1.7900469825946195, 1.8054167642140493, 1.8209529965910054,
4221.8366654602533845, 1.8525645117230873, 1.8686611409895424, 1.8849670357028696, 1.9014946531003181,
4231.9182573008597323, 1.9352692282919006, 1.9525457295488893, 1.9701032608497135, 1.9879595741230611,
4242.0061338699589673, 2.0246469733729340, 2.0435215366506698, 2.0627822745039639, 2.0824562379877247,
4252.1025731351849992, 2.1231657086697902, 2.1442701823562618, 2.1659267937448412, 2.1881804320720208,
4262.2110814088747279, 2.2346863955870573, 2.2590595738653296, 2.2842740596736570, 2.3104136836950024,
4272.3375752413355309, 2.3658713701139877, 2.3954342780074676, 2.4264206455302118, 2.4590181774083506,
4282.4934545220919508, 2.5300096723854670, 2.5690336259216395, 2.6109722484286135, 2.6564064112581929,
4292.7061135731187225, 2.7611693723841539, 2.8231253505459666, 2.8943440070186708, 2.9786962526450171,
4303.0832288582142140, 3.2230849845786187, 3.4426198558966523, 3.7130862467403638
431};
432static const double Y_ZIGGURRAT_GAUSS128[N_ZIGGURRAT_GAUSS128+1] = {
4331.0000000000000000, 0.9635996931557651, 0.9362826817083690, 0.9130436479920363, 0.8922816508023012,
4340.8732430489268521, 0.8555006078850629, 0.8387836053106459, 0.8229072113952607, 0.8077382946961199,
4350.7931770117838580, 0.7791460859417020, 0.7655841739092348, 0.7524415591857027, 0.7396772436833371,
4360.7272569183545049, 0.7151515074204761, 0.7033360990258165, 0.6917891434460349, 0.6804918410064135,
4370.6694276673577053, 0.6585820000586529, 0.6479418211185500, 0.6374954773431442, 0.6272324852578138,
4380.6171433708265618, 0.6072195366326042, 0.5974531509518116, 0.5878370544418199, 0.5783646811267017,
4390.5690299910747210, 0.5598274127106941, 0.5507517931210546, 0.5417983550317235, 0.5329626593899870,
4400.5242405726789923, 0.5156282382498716, 0.5071220510813041, 0.4987186354765838, 0.4904148252893212,
4410.4822076463348383, 0.4740943006982492, 0.4660721526945706, 0.4581387162728716, 0.4502916436869266,
4420.4425287152802462, 0.4348478302546615, 0.4272469983095620, 0.4197243320540379, 0.4122780401070242,
4430.4049064208114880, 0.3976078564980422, 0.3903808082413892, 0.3832238110598833, 0.3761354695144541,
4440.3691144536682749, 0.3621594953730330, 0.3552693848515469, 0.3484429675498723, 0.3416791412350135,
4450.3349768533169710, 0.3283350983761522, 0.3217529158792085, 0.3152293880681574, 0.3087636380092518,
4460.3023548277894796, 0.2960021568498557, 0.2897048604458103, 0.2834622082260124, 0.2772735029218976,
4470.2711380791410251, 0.2650553022581618, 0.2590245673987105, 0.2530452985097656, 0.2471169475146965,
4480.2412389935477511, 0.2354109422657275, 0.2296323252343025, 0.2239026993871337, 0.2182216465563704,
4490.2125887730737359, 0.2070037094418736, 0.2014661100762031, 0.1959756531181102, 0.1905320403209136,
4500.1851349970107133, 0.1797842721249620, 0.1744796383324022, 0.1692208922389246, 0.1640078546849276,
4510.1588403711409350, 0.1537183122095865, 0.1486415742436969, 0.1436100800919329, 0.1386237799858510,
4520.1336826525846476, 0.1287867061971039, 0.1239359802039817, 0.1191305467087185, 0.1143705124498882,
4530.1096560210158177, 0.1049872554103545, 0.1003644410295455, 0.0957878491225781, 0.0912578008276347,
4540.0867746718955429, 0.0823388982429574, 0.0779509825146547, 0.0736115018847548, 0.0693211173941802,
4550.0650805852136318, 0.0608907703485663, 0.0567526634815385, 0.0526674019035031, 0.0486362958602840,
4560.0446608622008724, 0.0407428680747906, 0.0368843887869688, 0.0330878861465051, 0.0293563174402538,
4570.0256932919361496, 0.0221033046161116, 0.0185921027371658, 0.0151672980106720, 0.0118394786579823,
4580.0086244844129305, 0.0055489952208165, 0.0026696290839025, 0.0000000000000000
459};
460
461
462r_8 RandomGeneratorInterface::GaussianZiggurat128()
463//--------
464// On a "N = 128" bandes horizontales numerotees [0,N-1=127]
465// Les tableaux ont une taille "N + 1 = 129"
466// On tire un numero de bande dans [0,N-1=127]
467//--------
468// Pour choisir le signe sans avoir a retirer un aleatoire,
469// on utilise un digit du premier tirage qui n'est pas utilise
470// dans le choix du numero de bande "I":
471// U = [0,1[ , Ai=[0,9] (chiffres)
472// U = A1/10 + A2/100 + A3/1000 + A4/10000 + A5/100000 + ...
473// 128*U = A1*12.8 + A2*1.28 + A3*0.128 + A4*0.0128 + A5*0.00128
474// pour Ai le plus grand possible cad 9
475// 128*U = 115.2 + 11.52 + 1.152 + 0.1152 + 0.01152
476// On voit que pour Ai = 9
477// 1 terme, A1 : I = int(128*U) = int(115.2) = 115
478// 2 termes, A1,2 : I = int(128*U) = int(115.2+11.52) = int(126.72) = 126
479// 3 termes, A1,2,3 : I = int(128*U) = int(115.2+11.52+1.152) = int(127.872) = 127
480// 4 termes, A1,2,3,4 : I = int(128*U) = int(115.2+11.52+1.152+0.1152) = int(127.9872) = 127
481// ==> le digit A4 ne sert pas dans la determination de "I"
482// On prend un digit de + pour avoir de la marge -> A5 cad le dernier digit de int(U*10^5)
483//--------
484{
485 while(1) {
486 double U;
487 // -- choix de l'intervalle et de l'abscisse "x"
488 int I = N_ZIGGURRAT_GAUSS128;
489 while(I>=N_ZIGGURRAT_GAUSS128) {
490 U = Next();
491 I = int(N_ZIGGURRAT_GAUSS128*U);
492 }
493
494 // -- choix du signe (cf blabla ci-dessus)
495 double s = ( (int(U*100000)&1) == 0 ) ? -1.: 1.;
496
497 // -- choix de l'abscisse "x" dans l'intervalle
498 double x = Next() * X_ZIGGURRAT_GAUSS128[I+1];
499
500 // -- x est dans l'interieur de la bande
501 if(x<X_ZIGGURRAT_GAUSS128[I]) return s * x;
502
503 // -- x n'est pas a l'interieur de la bande mais dans la partie a 2 possibilites
504
505 // l'intervalle est celui qui contient la queue a l'infini
506 // On s'assure que la partie "rejection sur la gaussienne" ne sera pas appelle
507 // (cad que slim=1. < X[127]) pour eviter les recursions infinies (possibles?)
508 if(I==N_ZIGGURRAT_GAUSS128-1) // cas de la bande de la queue I=127
509 return s * GaussianTail(X_ZIGGURRAT_GAUSS128[N_ZIGGURRAT_GAUSS128-1],1.);
510
511 // on tire "y" uniforme a l'interieur des ordonnees de la bande choisie
512 // et on regarde si on est en-dessous ou au-dessus de la pdf
513 double y = Y_ZIGGURRAT_GAUSS128[I+1]
514 + Next()*(Y_ZIGGURRAT_GAUSS128[I]-Y_ZIGGURRAT_GAUSS128[I+1]);
515 double pdf = exp(-0.5*x*x);
516 if(pdf>y) return s * x;
517
518 // echec, on est au-dessus de la pdf -> on re-essaye
519 }
520}
521
522r_8 RandomGeneratorInterface::GaussianTail(double s,double slim)
523{
524 /* Returns a gaussian random variable larger than a
525 * This implementation does one-sided upper-tailed deviates.
526 */
527
528 if (s < slim)
529 {
530 /* For small s, use a direct rejection method. The limit s < 1
531 can be adjusted to optimise the overall efficiency */
532 double x;
533 do
534 {
535 x = Gaussian();
536 }
537 while (x < s);
538 return x;
539 }
540 else
541 {
542 /* Use the "supertail" deviates from the last two steps
543 * of Marsaglia's rectangle-wedge-tail method, as described
544 * in Knuth, v2, 3rd ed, pp 123-128. (See also exercise 11, p139,
545 * and the solution, p586.)
546 */
547 double u, v, x;
548 do
549 {
550 u = Next();
551 do
552 {
553 v = Next();
554 }
555 while (v == 0.0);
556 x = sqrt (s * s - 2 * log (v));
557 }
558 while (x * u > s);
559 return x;
560 }
561}
562
563/////////////////////////////////////////////////////////////////////////
564/////////////////////////////////////////////////////////////////////////
565/////////////////////////////////////////////////////////////////////////
566
567uint_8 RandomGeneratorInterface::Poisson(double mu, double mumax)
568{
569 switch (usepoisson_) {
570 case C_Poisson_Simple :
571 return PoissonSimple(mu,mumax);
572 break;
573 case C_Poisson_Ahrens :
574 return PoissonAhrens(mu);
575 break;
576 default:
577 return PoissonSimple(mu,mumax);
578 break;
579 }
580}
581
582
583//--- Generation de nombre aleatoires suivant une distribution de Poisson
584uint_8 RandomGeneratorInterface::PoissonSimple(double mu,double mumax)
585{
586 double pp,ppi,x;
587
588 if((mumax>0.)&&(mu>=mumax)) {
589 pp = sqrt(mu);
590 while( (x=pp*Gaussian()) < -mu );
591 return (uint_8)(mu+x+0.5);
592 }
593 else {
594 uint_8 n;
595 ppi = pp = exp(-mu);
596 x = Next();
597 n = 0;
598 while (x > ppi) {
599 n++;
600 pp = mu*pp/(double)n;
601 ppi += pp;
602 }
603 return n;
604 }
605 return 0; // pas necessaire ?
606}
607
608
609static double a0_poiahr = -0.5;
610static double a1_poiahr = 0.3333333;
611static double a2_poiahr = -0.2500068;
612static double a3_poiahr = 0.2000118;
613static double a4_poiahr = -0.1661269;
614static double a5_poiahr = 0.1421878;
615static double a6_poiahr = -0.1384794;
616static double a7_poiahr = 0.125006;
617static double fact_poiahr[10] = {
618 1.0,1.0,2.0,6.0,24.0,120.0,720.0,5040.0,40320.0,362880.0};
619uint_8 RandomGeneratorInterface::PoissonAhrens(double mu)
620/*
621**********************************************************************
622 long ignpoi(float mu)
623 GENerate POIsson random deviate
624 Function
625 Generates a single random deviate from a Poisson
626 distribution with mean AV.
627 Arguments
628 av --> The mean of the Poisson distribution from which
629 a random deviate is to be generated.
630 genexp <-- The random deviate.
631 Method
632 Renames KPOIS from TOMS as slightly modified by BWB to use RANF
633 instead of SUNIF.
634 For details see:
635 Ahrens, J.H. and Dieter, U.
636 Computer Generation of Poisson Deviates
637 From Modified Normal Distributions.
638 ACM Trans. Math. Software, 8, 2
639 (June 1982),163-179
640**********************************************************************
641**********************************************************************
642
643
644 P O I S S O N DISTRIBUTION
645
646
647**********************************************************************
648**********************************************************************
649
650 FOR DETAILS SEE:
651
652 AHRENS, J.H. AND DIETER, U.
653 COMPUTER GENERATION OF POISSON DEVIATES
654 FROM MODIFIED NORMAL DISTRIBUTIONS.
655 ACM TRANS. MATH. SOFTWARE, 8,2 (JUNE 1982), 163 - 179.
656
657 (SLIGHTLY MODIFIED VERSION OF THE PROGRAM IN THE ABOVE ARTICLE)
658
659**********************************************************************
660 INTEGER FUNCTION IGNPOI(IR,MU)
661 INPUT: IR=CURRENT STATE OF BASIC RANDOM NUMBER GENERATOR
662 MU=MEAN MU OF THE POISSON DISTRIBUTION
663 OUTPUT: IGNPOI=SAMPLE FROM THE POISSON-(MU)-DISTRIBUTION
664 MUPREV=PREVIOUS MU, MUOLD=MU AT LAST EXECUTION OF STEP P OR B.
665 TABLES: COEFFICIENTS A0-A7 FOR STEP F. FACTORIALS FACT
666 COEFFICIENTS A(K) - FOR PX = FK*V*V*SUM(A(K)*V**K)-DEL
667 SEPARATION OF CASES A AND B
668*/
669{
670uint_8 ignpoi,j,k,kflag,l,m;
671double b1,b2,c,c0,c1,c2,c3,d,del,difmuk,e,fk,fx,fy,g,omega,p,p0,px,py,q,s,
672 t,u,v,x,xx,pp[35];
673
674 if(mu < 10.0) goto S120;
675/*
676 C A S E A. (RECALCULATION OF S,D,L IF MU HAS CHANGED)
677*/
678 s = sqrt(mu);
679 d = 6.0*mu*mu;
680/*
681 THE POISSON PROBABILITIES PK EXCEED THE DISCRETE NORMAL
682 PROBABILITIES FK WHENEVER K >= M(MU). L=IFIX(MU-1.1484)
683 IS AN UPPER BOUND TO M(MU) FOR ALL MU >= 10 .
684*/
685 l = (uint_8) (mu-1.1484);
686/*
687 STEP N. NORMAL SAMPLE - SNORM(IR) FOR STANDARD NORMAL DEVIATE
688*/
689 g = mu+s*Gaussian();
690 if(g < 0.0) goto S20;
691 ignpoi = (uint_8) (g);
692/*
693 STEP I. IMMEDIATE ACCEPTANCE IF IGNPOI IS LARGE ENOUGH
694*/
695 if(ignpoi >= l) return ignpoi;
696/*
697 STEP S. SQUEEZE ACCEPTANCE - SUNIF(IR) FOR (0,1)-SAMPLE U
698*/
699 fk = (double)ignpoi;
700 difmuk = mu-fk;
701 u = Next();
702 if(d*u >= difmuk*difmuk*difmuk) return ignpoi;
703S20:
704/*
705 STEP P. PREPARATIONS FOR STEPS Q AND H.
706 (RECALCULATIONS OF PARAMETERS IF NECESSARY)
707 .3989423=(2*PI)**(-.5) .416667E-1=1./24. .1428571=1./7.
708 THE QUANTITIES B1, B2, C3, C2, C1, C0 ARE FOR THE HERMITE
709 APPROXIMATIONS TO THE DISCRETE NORMAL PROBABILITIES FK.
710 C=.1069/MU GUARANTEES MAJORIZATION BY THE 'HAT'-FUNCTION.
711*/
712 omega = 0.3989423/s;
713 b1 = 4.166667E-2/mu;
714 b2 = 0.3*b1*b1;
715 c3 = 0.1428571*b1*b2;
716 c2 = b2-15.0*c3;
717 c1 = b1-6.0*b2+45.0*c3;
718 c0 = 1.0-b1+3.0*b2-15.0*c3;
719 c = 0.1069/mu;
720 if(g < 0.0) goto S50;
721/*
722 'SUBROUTINE' F IS CALLED (KFLAG=0 FOR CORRECT RETURN)
723*/
724 kflag = 0;
725 goto S70;
726S40:
727/*
728 STEP Q. QUOTIENT ACCEPTANCE (RARE CASE)
729*/
730 if(fy-u*fy <= py*exp(px-fx)) return ignpoi;
731S50:
732/*
733 STEP E. EXPONENTIAL SAMPLE - SEXPO(IR) FOR STANDARD EXPONENTIAL
734 DEVIATE E AND SAMPLE T FROM THE LAPLACE 'HAT'
735 (IF T <= -.6744 THEN PK < FK FOR ALL MU >= 10.)
736*/
737 e = Exponential();
738 u = Next();
739 u += (u-1.0);
740 //t = 1.8+fsign(e,u);
741 t = 1.8 + (((u>0. && e<0.) || (u<0. && e>0.))?-e:e);
742 if(t <= -0.6744) goto S50;
743 ignpoi = (uint_8) (mu+s*t);
744 fk = (double)ignpoi;
745 difmuk = mu-fk;
746/*
747 'SUBROUTINE' F IS CALLED (KFLAG=1 FOR CORRECT RETURN)
748*/
749 kflag = 1;
750 goto S70;
751S60:
752/*
753 STEP H. HAT ACCEPTANCE (E IS REPEATED ON REJECTION)
754*/
755 if(c*fabs(u) > py*exp(px+e)-fy*exp(fx+e)) goto S50;
756 return ignpoi;
757S70:
758/*
759 STEP F. 'SUBROUTINE' F. CALCULATION OF PX,PY,FX,FY.
760 CASE IGNPOI .LT. 10 USES FACTORIALS FROM TABLE FACT
761*/
762 if(ignpoi >= 10) goto S80;
763 px = -mu;
764 py = pow(mu,(double)ignpoi)/ *(fact_poiahr+ignpoi);
765 goto S110;
766S80:
767/*
768 CASE IGNPOI .GE. 10 USES POLYNOMIAL APPROXIMATION
769 A0-A7 FOR ACCURACY WHEN ADVISABLE
770 .8333333E-1=1./12. .3989423=(2*PI)**(-.5)
771*/
772 del = 8.333333E-2/fk;
773 del -= (4.8*del*del*del);
774 v = difmuk/fk;
775 if(fabs(v) <= 0.25) goto S90;
776 px = fk*log(1.0+v)-difmuk-del;
777 goto S100;
778S90:
779 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;
780S100:
781 py = 0.3989423/sqrt(fk);
782S110:
783 x = (0.5-difmuk)/s;
784 xx = x*x;
785 fx = -0.5*xx;
786 fy = omega*(((c3*xx+c2)*xx+c1)*xx+c0);
787 if(kflag <= 0) goto S40;
788 goto S60;
789S120:
790/*
791 C A S E B. (START NEW TABLE AND CALCULATE P0 IF NECESSARY)
792*/
793// m = max(1L,(long) (mu));
794 m = (1ULL >= (uint_8)mu) ? 1ULL: (uint_8)mu;
795
796 l = 0;
797 p = exp(-mu);
798 q = p0 = p;
799S130:
800/*
801 STEP U. UNIFORM SAMPLE FOR INVERSION METHOD
802*/
803 u = Next();
804 ignpoi = 0;
805 if(u <= p0) return ignpoi;
806/*
807 STEP T. TABLE COMPARISON UNTIL THE END PP(L) OF THE
808 PP-TABLE OF CUMULATIVE POISSON PROBABILITIES
809 (0.458=PP(9) FOR MU=10)
810*/
811 if(l == 0) goto S150;
812 j = 1;
813//if(u > 0.458) j = min(l,m);
814 if(u > 0.458) j = ((l<=m)? l: m);
815 for(k=j; k<=l; k++) {
816 if(u <= *(pp+k-1)) goto S180;
817 }
818 if(l == 35) goto S130;
819S150:
820/*
821 STEP C. CREATION OF NEW POISSON PROBABILITIES P
822 AND THEIR CUMULATIVES Q=PP(K)
823*/
824 l += 1;
825 for(k=l; k<=35; k++) {
826 p = p*mu/(double)k;
827 q += p;
828 *(pp+k-1) = q;
829 if(u <= q) goto S170;
830 }
831 l = 35;
832 goto S130;
833S170:
834 l = k;
835S180:
836 ignpoi = k;
837 return ignpoi;
838}
839
840/////////////////////////////////////////////////////////////////////////
841/////////////////////////////////////////////////////////////////////////
842/////////////////////////////////////////////////////////////////////////
843
844r_8 RandomGeneratorInterface::Exponential()
845{
846 switch (useexpo_) {
847 case C_Exponential_Simple :
848 return ExpoSimple();
849 break;
850 case C_Exponential_Ahrens :
851 return ExpoAhrens();
852 break;
853 default:
854 return ExpoSimple();
855 break;
856 }
857}
858
859r_8 RandomGeneratorInterface::ExpoSimple(void)
860{
861 return -log(1.-Next());
862}
863
864
865static double q_expo[8] = {
866 0.6931472,0.9333737,0.9888778,0.9984959,0.9998293,0.9999833,0.9999986,1.0};
867r_8 RandomGeneratorInterface::ExpoAhrens(void)
868/*
869**********************************************************************
870**********************************************************************
871 (STANDARD-) E X P O N E N T I A L DISTRIBUTION
872**********************************************************************
873**********************************************************************
874
875 FOR DETAILS SEE:
876
877 AHRENS, J.H. AND DIETER, U.
878 COMPUTER METHODS FOR SAMPLING FROM THE
879 EXPONENTIAL AND NORMAL DISTRIBUTIONS.
880 COMM. ACM, 15,10 (OCT. 1972), 873 - 882.
881
882 ALL STATEMENT NUMBERS CORRESPOND TO THE STEPS OF ALGORITHM
883 'SA' IN THE ABOVE PAPER (SLIGHTLY MODIFIED IMPLEMENTATION)
884
885 Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of
886 SUNIF. The argument IR thus goes away.
887
888**********************************************************************
889 Q(N) = SUM(ALOG(2.0)**K/K!) K=1,..,N , THE HIGHEST N
890 (HERE 8) IS DETERMINED BY Q(N)=1.0 WITHIN STANDARD PRECISION
891*/
892{
893long i;
894double sexpo,a,u,ustar,umin;
895double *q1 = q_expo;
896 a = 0.0;
897 while((u=Next())==0.);
898 goto S30;
899S20:
900 a += *q1;
901S30:
902 u += u;
903 if(u <= 1.0) goto S20;
904 u -= 1.0;
905 if(u > *q1) goto S60;
906 sexpo = a+u;
907 return sexpo;
908S60:
909 i = 1;
910 ustar = Next();
911 umin = ustar;
912S70:
913 ustar = Next();
914 if(ustar < umin) umin = ustar;
915 i += 1;
916 if(u > *(q_expo+i-1)) goto S70;
917 sexpo = a+umin**q1;
918 return sexpo;
919}
920
921/////////////////////////////////////////////////////////////////////////
922/////////////////////////////////////////////////////////////////////////
923/////////////////////////////////////////////////////////////////////////
924
925int RandomGeneratorInterface::Gaussian2DRho(double &x,double &y,double mx,double my,double sx,double sy,double ro)
926/*
927++
928| Tirage de 2 nombres aleatoires x et y distribues sur une gaussienne 2D
929| de centre (mx,my), de coefficient de correlation rho (ro) et telle que
930| les sigmas finals des variables x et y soient sx,sy (ce sont
931| les valeurs des distributions marginales des variables aleatoires x et y
932| c'est a dire les sigmas des projections x et y de l'histogramme 2D
933| de la gaussienne). Retourne 0 si ok.
934|
935| - La densite de probabilite (normalisee a 1) sur laquelle on tire est:
936| N*exp[-0.5*{[(dx/sx)^2-2*ro/(sx*sy)*dx*dy+(dy/sy)^2]/(1-ro^2)}]
937| avec dx = x-mx, dy = y-my et N = 1/[2Pi*sx*sy*sqrt(1-ro^2)]
938| - Dans ce cas la distribution marginale est (ex en X):
939| 1/(sqrt(2Pi)*sx) * exp[-0.5*{dx^2/sx^2}]
940| - La matrice des covariances C des variables x,y est:
941| | sx^2 ro*sx*sy |
942| | | et det(C) = (1-ro^2)*sx^2*sy^2
943| | ro*sx*sy sy^2 |
944| - La matrice inverse C^(-1) est:
945| | 1/sx^2 -ro/(sx*sy) |
946| | | * 1/(1-ro^2)
947| | -ro/(sx*sy) 1/sy^2 |
948|
949| - Remarque:
950| le sigma que l'on obtient quand on fait une coupe de la gaussienne 2D
951| en y=0 (ou x=0) est: SX0(y=0) = sx*sqrt(1-ro^2) different de sx
952| SY0(x=0) = sy*sqrt(1-ro^2) different de sy
953| La distribution qui correspond a des sigmas SX0,SY0
954| pour les coupes en y=0,x=0 de la gaussienne 2D serait:
955| N*exp[-0.5*{ (dx/SX0)^2-2*ro/(SX0*SY0)*dx*dy+(dy/SY0)^2 }]
956| avec N = sqrt(1-ro^2)/(2Pi*SX0*SY0) et les variances
957| des variables x,y sont toujours
958| sx=SX0/sqrt(1-ro^2), sy=SY0/sqrt(1-ro^2)
959--
960*/
961{
962double a,b,sa;
963
964if( ro <= -1. || ro >= 1. ) return 1;
965
966while( (b=Flat01()) == 0. );
967b = sqrt(-2.*log(b));
968a = 2.*M_PI * Flat01();
969sa = sin(a);
970
971x = mx + sx*b*(sqrt(1.-ro*ro)*cos(a)+ro*sa);
972y = my + sy*b*sa;
973
974return 0;
975}
976
977void RandomGeneratorInterface::Gaussian2DAng(double &x,double &y,double mx,double my,double sa,double sb,double teta)
978/*
979++
980| Tirage de 2 nombres aleatoires x et y distribues sur une gaussienne 2D
981| de centre (x=mx,y=my), de sigmas grand axe et petit axe (sa,sb)
982| et dont le grand axe fait un angle teta (radian) avec l'axe des x.
983|
984| - La densite de probabilite (normalisee a 1) sur laquelle on tire est:
985| N*exp[-0.5*{ (A/sa)**2+(C/sc)**2 }], N=1/(2Pi*sa*sc)
986| ou A et B sont les coordonnees selon le grand axe et le petit axe
987| et teta = angle(x,A), le resultat subit ensuite une rotation d'angle teta.
988| - La matrice des covariances C des variables A,B est:
989| | sa^2 0 |
990| | | et det(C) = (1-ro^2)*sa^2*sb^2
991| | 0 sb^2 |
992| - La distribution x,y resultante est:
993| N*exp[-0.5*{[(dx/sx)^2-2*ro/(sx*sy)*dx*dy+(dy/sy)^2]/(1-ro^2)}]
994| ou N est donne dans NormCo et sx,sy,ro sont calcules a partir
995| de sa,sc,teta (voir fonctions paramga ou gaparam). La matrice des
996| covariances des variables x,y est donnee dans la fonction NormCo.
997--
998*/
999{
1000double c,s,X,Y;
1001
1002while( (s = Flat01()) == 0. );
1003s = sqrt(-2.*log(s));
1004c = 2.*M_PI * Flat01();
1005
1006X = sa*s*cos(c);
1007Y = sb*s*sin(c);
1008
1009c = cos(teta); s = sin(teta);
1010x = mx + c*X - s*Y;
1011y = my + s*X + c*Y;
1012}
1013
1014} /* namespace SOPHYA */
1015
1016
1017
1018/////////////////////////////////////////////////////////////////
1019/*
1020**** Remarques sur complex< r_8 > ComplexGaussian(double sig) ****
1021
1022--- variables gaussiennes x,y independantes
1023x gaussien: pdf f(x) = 1/(sqrt(2Pi) Sx) exp(-(x-Mx)^2/(2 Sx^2))
1024y gaussien: pdf f(y) = 1/(sqrt(2Pi) Sy) exp(-(y-My)^2/(2 Sy^2))
1025x,y independants --> pdf f(x,y) = f(x) f(y)
1026On a:
1027 <x> = Integrate[x*f(x)] = Mx
1028 <x^2> = Integrate[x^2*f(x)] = Mx^2 + Sx^2
1029
1030--- On cherche la pdf g(r,t) du module et de la phase
1031 x = r cos(t) , y = r sin(t)
1032 r=sqrt(x^2+y^2 , t=atan2(y,x)
1033 (r,t) --> (x,y): le Jacobien = r
1034
1035 g(r,t) = r f(x,y) = r f(x) f(y)
1036 = r/(2Pi Sx Sy) exp(-(x-Mx)^2/(2 Sx^2)) exp(-(y-My)^2/(2 Sy^2))
1037
1038- Le cas general est complique
1039 (cf D.Pelat cours DEA "bruits et signaux" section 4.5)
1040
1041- Cas ou "Mx = My = 0" et "Sx = Sy = S"
1042 c'est la pdf du module et de la phase d'un nombre complexe
1043 dont les parties reelles et imaginaires sont independantes
1044 et sont distribuees selon des gaussiennes de variance S^2
1045 g(r,t) = r/(2Pi S^2) exp(-r^2/(2 S^2))
1046 La distribution de "r" est donc:
1047 g(r) = Integrate[g(r,t),{t,0,2Pi}]
1048 = r/S^2 exp(-r^2/(2 S^2))
1049 La distribution de "t" est donc:
1050 g(t) = Integrate[g(r,t),{r,0,Infinity}]
1051 = 1 / 2Pi (distribution uniforme sur [0,2Pi[)
1052 Les variables aleatoires r,t sont independantes:
1053 g(r,t) = g(r) g(t)
1054On a:
1055 <r> = Integrate[r*g(r)] = sqrt(PI/2)*S
1056 <r^2> = Integrate[r^2*g(r)] = 2*S^2
1057 <r^3> = Integrate[r^3*g(r)] = 3*sqrt(PI/2)*S^3
1058 <r^4> = Integrate[r^4*g(r)] = 8*S^4
1059
1060- Attention:
1061La variable complexe "c = x+iy = r*exp(i*t)" ainsi definie verifie:
1062 <|c|^2> = <c c*> = <x^2+y^2> = <r^2> = 2 S^2
1063Si on veut generer une variable complexe gaussienne telle que
1064 <c c*> = s^2 alors il faut prendre S = s/sqrt(2) comme argument
1065
1066*/
Note: See TracBrowser for help on using the repository browser.