ROOT  6.06/09
Reference Guide
TRandom1.cxx
Go to the documentation of this file.
1 // @(#)root/mathcore:$Id$
2 // Author: Rene Brun from CLHEP & CERNLIB 04/05/2006
3 
4 /**
5 
6 \class TRandom1
7 
8 The Ranlux Random number generator class
9 
10 The algorithm for this random engine has been taken from the original
11 implementation in FORTRAN by Fred James as part of CLHEP.
12 
13 The initialisation is carried out using a Multiplicative Congruential
14 generator using formula constants of L'Ecuyer as described in "F.James,
15 Comp. Phys. Comm. 60 (1990) 329-344".
16 
17 @ingroup Random
18 
19 */
20 
21 
22 #include "TRandom1.h"
23 #include "TRandom3.h"
24 #include "TMath.h"
25 #include <stdlib.h>
26 
27 // Number of instances with automatic seed selection
29 
30 // Maximum index into the seed table
31 int TRandom1::fgMaxIndex = 215;
32 #ifndef __CINT__
33 const UInt_t fgSeedTable[215][2] = {
34  { 9876, 54321 },
35  { 1299961164, 253987020 },
36  { 669708517, 2079157264 },
37  { 190904760, 417696270 },
38  { 1289741558, 1376336092 },
39  { 1803730167, 324952955 },
40  { 489854550, 582847132 },
41  { 1348037628, 1661577989 },
42  { 350557787, 1155446919 },
43  { 591502945, 634133404 },
44  { 1901084678, 862916278 },
45  { 1988640932, 1785523494 },
46  { 1873836227, 508007031 },
47  { 1146416592, 967585720 },
48  { 1837193353, 1522927634 },
49  { 38219936, 921609208 },
50  { 349152748, 112892610 },
51  { 744459040, 1735807920 },
52  { 1983990104, 728277902 },
53  { 309164507, 2126677523 },
54  { 362993787, 1897782044 },
55  { 556776976, 462072869 },
56  { 1584900822, 2019394912 },
57  { 1249892722, 791083656 },
58  { 1686600998, 1983731097 },
59  { 1127381380, 198976625 },
60  { 1999420861, 1810452455 },
61  { 1972906041, 664182577 },
62  { 84636481, 1291886301 },
63  { 1186362995, 954388413 },
64  { 2141621785, 61738584 },
65  { 1969581251, 1557880415 },
66  { 1150606439, 136325185 },
67  { 95187861, 1592224108 },
68  { 940517655, 1629971798 },
69  { 215350428, 922659102 },
70  { 786161212, 1121345074 },
71  { 1450830056, 1922787776 },
72  { 1696578057, 2025150487 },
73  { 1803414346, 1851324780 },
74  { 1017898585, 1452594263 },
75  { 1184497978, 82122239 },
76  { 633338765, 1829684974 },
77  { 430889421, 230039326 },
78  { 492544653, 76320266 },
79  { 389386975, 1314148944 },
80  { 1720322786, 709120323 },
81  { 1868768216, 1992898523 },
82  { 443210610, 811117710 },
83  { 1191938868, 1548484733 },
84  { 616890172, 159787986 },
85  { 935835339, 1231440405 },
86  { 1058009367, 1527613300 },
87  { 1463148129, 1970575097 },
88  { 1795336935, 434768675 },
89  { 274019517, 605098487 },
90  { 483689317, 217146977 },
91  { 2070804364, 340596558 },
92  { 930226308, 1602100969 },
93  { 989324440, 801809442 },
94  { 410606853, 1893139948 },
95  { 1583588576, 1219225407 },
96  { 2102034391, 1394921405 },
97  { 2005037790, 2031006861 },
98  { 1244218766, 923231061 },
99  { 49312790, 775496649 },
100  { 721012176, 321339902 },
101  { 1719909107, 1865748178 },
102  { 1156177430, 1257110891 },
103  { 307561322, 1918244397 },
104  { 906041433, 360476981 },
105  { 1591375755, 268492659 },
106  { 461522398, 227343256 },
107  { 2145930725, 2020665454 },
108  { 1938419274, 1331283701 },
109  { 174405412, 524140103 },
110  { 494343653, 18063908 },
111  { 1025534808, 181709577 },
112  { 2048959776, 1913665637 },
113  { 950636517, 794796256 },
114  { 1828843197, 1335757744 },
115  { 211109723, 983900607 },
116  { 825474095, 1046009991 },
117  { 374915657, 381856628 },
118  { 1241296328, 698149463 },
119  { 1260624655, 1024538273 },
120  { 900676210, 1628865823 },
121  { 697951025, 500570753 },
122  { 1007920268, 1708398558 },
123  { 264596520, 624727803 },
124  { 1977924811, 674673241 },
125  { 1440257718, 271184151 },
126  { 1928778847, 993535203 },
127  { 1307807366, 1801502463 },
128  { 1498732610, 300876954 },
129  { 1617712402, 1574250679 },
130  { 1261800762, 1556667280 },
131  { 949929273, 560721070 },
132  { 1766170474, 1953522912 },
133  { 1849939248, 19435166 },
134  { 887262858, 1219627824 },
135  { 483086133, 603728993 },
136  { 1330541052, 1582596025 },
137  { 1850591475, 723593133 },
138  { 1431775678, 1558439000 },
139  { 922493739, 1356554404 },
140  { 1058517206, 948567762 },
141  { 709067283, 1350890215 },
142  { 1044787723, 2144304941 },
143  { 999707003, 513837520 },
144  { 2140038663, 1850568788 },
145  { 1803100150, 127574047 },
146  { 867445693, 1149173981 },
147  { 408583729, 914837991 },
148  { 1166715497, 602315845 },
149  { 430738528, 1743308384 },
150  { 1388022681, 1760110496 },
151  { 1664028066, 654300326 },
152  { 1767741172, 1338181197 },
153  { 1625723550, 1742482745 },
154  { 464486085, 1507852127 },
155  { 754082421, 1187454014 },
156  { 1315342834, 425995190 },
157  { 960416608, 2004255418 },
158  { 1262630671, 671761697 },
159  { 59809238, 103525918 },
160  { 1205644919, 2107823293 },
161  { 1615183160, 1152411412 },
162  { 1024474681, 2118672937 },
163  { 1703877649, 1235091369 },
164  { 1821417852, 1098463802 },
165  { 1738806466, 1529062843 },
166  { 620780646, 1654833544 },
167  { 1070174101, 795158254 },
168  { 658537995, 1693620426 },
169  { 2055317555, 508053916 },
170  { 1647371686, 1282395762 },
171  { 29067379, 409683067 },
172  { 1763495989, 1917939635 },
173  { 1602690753, 810926582 },
174  { 885787576, 513818500 },
175  { 1853512561, 1195205756 },
176  { 1798585498, 1970460256 },
177  { 1819261032, 1306536501 },
178  { 1133245275, 37901 },
179  { 689459799, 1334389069 },
180  { 1730609912, 1854586207 },
181  { 1556832175, 1228729041 },
182  { 251375753, 683687209 },
183  { 2083946182, 1763106152 },
184  { 2142981854, 1365385561 },
185  { 763711891, 1735754548 },
186  { 1581256466, 173689858 },
187  { 2121337132, 1247108250 },
188  { 1004003636, 891894307 },
189  { 569816524, 358675254 },
190  { 626626425, 116062841 },
191  { 632086003, 861268491 },
192  { 1008211580, 779404957 },
193  { 1134217766, 1766838261 },
194  { 1423829292, 1706666192 },
195  { 942037869, 1549358884 },
196  { 1959429535, 480779114 },
197  { 778311037, 1940360875 },
198  { 1531372185, 2009078158 },
199  { 241935492, 1050047003 },
200  { 272453504, 1870883868 },
201  { 390441332, 1057903098 },
202  { 1230238834, 1548117688 },
203  { 1242956379, 1217296445 },
204  { 515648357, 1675011378 },
205  { 364477932, 355212934 },
206  { 2096008713, 1570161804 },
207  { 1409752526, 214033983 },
208  { 1288158292, 1760636178 },
209  { 407562666, 1265144848 },
210  { 1071056491, 1582316946 },
211  { 1014143949, 911406955 },
212  { 203080461, 809380052 },
213  { 125647866, 1705464126 },
214  { 2015685843, 599230667 },
215  { 1425476020, 668203729 },
216  { 1673735652, 567931803 },
217  { 1714199325, 181737617 },
218  { 1389137652, 678147926 },
219  { 288547803, 435433694 },
220  { 200159281, 654399753 },
221  { 1580828223, 1298308945 },
222  { 1832286107, 169991953 },
223  { 182557704, 1046541065 },
224  { 1688025575, 1248944426 },
225  { 1508287706, 1220577001 },
226  { 36721212, 1377275347 },
227  { 1968679856, 1675229747 },
228  { 279109231, 1835333261 },
229  { 1358617667, 1416978076 },
230  { 740626186, 2103913602 },
231  { 1882655908, 251341858 },
232  { 648016670, 1459615287 },
233  { 780255321, 154906988 },
234  { 857296483, 203375965 },
235  { 1631676846, 681204578 },
236  { 1906971307, 1623728832 },
237  { 1541899600, 1168449797 },
238  { 1267051693, 1020078717 },
239  { 1998673940, 1298394942 },
240  { 1914117058, 1381290704 },
241  { 426068513, 1381618498 },
242  { 139365577, 1598767734 },
243  { 2129910384, 952266588 },
244  { 661788054, 19661356 },
245  { 1104640222, 240506063 },
246  { 356133630, 1676634527 },
247  { 242242374, 1863206182 },
248  { 957935844, 1490681416 }};
249 #endif
250 
252 
253 ////////////////////////////////////////////////////////////////////////////////
254 /// Luxury level is set in the same way as the original FORTRAN routine.
255 /// level 0 (p=24): equivalent to the original RCARRY of Marsaglia
256 /// and Zaman, very long period, but fails many tests.
257 /// level 1 (p=48): considerable improvement in quality over level 0,
258 /// now passes the gap test, but still fails spectral test.
259 /// level 2 (p=97): passes all known tests, but theoretically still
260 /// defective.
261 /// level 3 (p=223): DEFAULT VALUE. Any theoretically possible
262 /// correlations have very small chance of being observed.
263 /// level 4 (p=389): highest possible luxury, all 24 bits chaotic.
264 
265 TRandom1::TRandom1(UInt_t seed, Int_t lux)
266  : fIntModulus(0x1000000),
267  fMantissaBit24( TMath::Power(0.5,24.) ),
268  fMantissaBit12( TMath::Power(0.5,12.) )
269 {
270  UInt_t seedlist[2]={0,0};
271 
272  fTheSeeds = &fSeed;
273  fLuxury = lux;
274  SetSeed2(seed, fLuxury);
275  // in case seed = 0 SetSeed2 calls already SetSeeds
276  if (seed != 0) {
277  // setSeeds() wants a zero terminated array!
278  seedlist[0]=fSeed;
279  seedlist[1]=0;
280  SetSeeds(seedlist, fLuxury);
281  }
282 }
283 
284 ////////////////////////////////////////////////////////////////////////////////
285 ///default constructor
286 
288  : fIntModulus(0x1000000),
289  fMantissaBit24( TMath::Power(0.5,24.) ),
290  fMantissaBit12( TMath::Power(0.5,12.) )
291 {
292  fTheSeeds = &fSeed;
293  UInt_t seed;
294  UInt_t seedlist[2]={0,0};
295 
296  fLuxury = 3;
297  int cycle = abs(int(fgNumEngines/fgMaxIndex));
298  int curIndex = abs(int(fgNumEngines%fgMaxIndex));
299  fgNumEngines +=1;
300  UInt_t mask = ((cycle & 0x007fffff) << 8);
301  GetTableSeeds( seedlist, curIndex );
302  seed = seedlist[0]^mask;
303  SetSeed2(seed, fLuxury);
304 
305  // setSeeds() wants a zero terminated array!
306  seedlist[0]=fSeed; //<=============
307  seedlist[1]=0;
308  SetSeeds(seedlist, fLuxury);
309 }
310 
311 ////////////////////////////////////////////////////////////////////////////////
312 ///constructor
313 
314 TRandom1::TRandom1(int rowIndex, int colIndex, int lux)
315  : fIntModulus(0x1000000),
316  fMantissaBit24( TMath::Power(0.5,24.) ),
317  fMantissaBit12( TMath::Power(0.5,12.) )
318 {
319  fTheSeeds = &fSeed;
320  UInt_t seed;
321  UInt_t seedlist[2]={0,0};
322 
323  fLuxury = lux;
324  int cycle = abs(int(rowIndex/fgMaxIndex));
325  int row = abs(int(rowIndex%fgMaxIndex));
326  int col = abs(int(colIndex%2));
327  UInt_t mask = (( cycle & 0x000007ff ) << 20 );
328  GetTableSeeds( seedlist, row );
329  seed = ( seedlist[col] )^mask;
330  SetSeed2(seed, fLuxury);
331 
332  // setSeeds() wants a zero terminated array!
333  seedlist[0]=fSeed;
334  seedlist[1]=0;
335  SetSeeds(seedlist, fLuxury);
336 }
337 
338 ////////////////////////////////////////////////////////////////////////////////
339 ///destructor
340 
342 {
343 }
344 
345 ////////////////////////////////////////////////////////////////////////////////
346 ///static function returning the table of seeds
347 
349 {
350  if ((index >= 0) && (index < 215)) {
351  seeds[0] = fgSeedTable[index][0];
352  seeds[1] = fgSeedTable[index][1];
353  }
354  else seeds = 0;
355 }
356 
357 ////////////////////////////////////////////////////////////////////////////////
358 ///return a random number in ]0,1]
359 
361 {
362  float next_random;
363  float uni;
364  int i;
365 
367  if(uni < 0. ) {
368  uni += 1.0;
369  fCarry = fMantissaBit24;
370  } else {
371  fCarry = 0.;
372  }
373 
374  fFloatSeedTable[fIlag] = uni;
375  fIlag --;
376  fJlag --;
377  if(fIlag < 0) fIlag = 23;
378  if(fJlag < 0) fJlag = 23;
379 
380  if( uni < fMantissaBit12 ){
382  if( uni == 0) uni = fMantissaBit24 * fMantissaBit24;
383  }
384  next_random = uni;
385  fCount24 ++;
386 
387 // every 24th number generation, several random numbers are generated
388 // and wasted depending upon the fLuxury level.
389 
390  if(fCount24 == 24 ) {
391  fCount24 = 0;
392  for( i = 0; i != fNskip ; i++) {
394  if(uni < 0. ) {
395  uni += 1.0;
396  fCarry = fMantissaBit24;
397  } else {
398  fCarry = 0.;
399  }
400  fFloatSeedTable[fIlag] = uni;
401  fIlag --;
402  fJlag --;
403  if(fIlag < 0)fIlag = 23;
404  if(fJlag < 0) fJlag = 23;
405  }
406  }
407  return (double) next_random;
408 }
409 
410 ////////////////////////////////////////////////////////////////////////////////
411 ///return an array of random numbers in ]0,1]
412 
413 void TRandom1::RndmArray(const Int_t size, Float_t *vect)
414 {
415  for (Int_t i=0;i<size;i++) vect[i] = Rndm();
416 }
417 
418 ////////////////////////////////////////////////////////////////////////////////
419 ///return an array of random numbers in ]0,1[
420 
421 void TRandom1::RndmArray(const Int_t size, Double_t *vect)
422 {
423  float next_random;
424  float uni;
425  int i;
426  int index;
427 
428  for (index=0; index<size; ++index) {
430  if(uni < 0. ) {
431  uni += 1.0;
432  fCarry = fMantissaBit24;
433  } else {
434  fCarry = 0.;
435  }
436 
437  fFloatSeedTable[fIlag] = uni;
438  fIlag --;
439  fJlag --;
440  if(fIlag < 0) fIlag = 23;
441  if(fJlag < 0) fJlag = 23;
442 
443  if( uni < fMantissaBit12 ){
445  if( uni == 0) uni = fMantissaBit24 * fMantissaBit24;
446  }
447  next_random = uni;
448  vect[index] = (double)next_random;
449  fCount24 ++;
450 
451 // every 24th number generation, several random numbers are generated
452 // and wasted depending upon the fLuxury level.
453 
454  if(fCount24 == 24 ) {
455  fCount24 = 0;
456  for( i = 0; i != fNskip ; i++) {
458  if(uni < 0. ) {
459  uni += 1.0;
460  fCarry = fMantissaBit24;
461  } else {
462  fCarry = 0.;
463  }
464  fFloatSeedTable[fIlag] = uni;
465  fIlag --;
466  fJlag --;
467  if(fIlag < 0)fIlag = 23;
468  if(fJlag < 0) fJlag = 23;
469  }
470  }
471  }
472 }
473 
474 
475 ////////////////////////////////////////////////////////////////////////////////
476 ///set seeds
477 
478 void TRandom1::SetSeeds(const UInt_t *seeds, int lux)
479 {
480  const int ecuyer_a = 53668;
481  const int ecuyer_b = 40014;
482  const int ecuyer_c = 12211;
483  const int ecuyer_d = 2147483563;
484 
485  const int lux_levels[5] = {0,24,73,199,365};
486  int i;
487  UInt_t int_seed_table[24];
488  Long64_t k_multiple,next_seed;
489  const UInt_t *seedptr;
490 
491  fTheSeeds = seeds;
492  seedptr = seeds;
493 
494  if(seeds == 0) {
495  SetSeed2(fSeed,lux);
496  fTheSeeds = &fSeed;
497  return;
498  }
499 
500  fSeed = *seeds;
501 
502 // number of additional random numbers that need to be 'thrown away'
503 // every 24 numbers is set using fLuxury level variable.
504 
505  if( (lux > 4)||(lux < 0) ) {
506  if(lux >= 24) {
507  fNskip = lux - 24;
508  } else {
509  fNskip = lux_levels[3]; // corresponds to default fLuxury level
510  }
511  } else {
512  fLuxury = lux;
513  fNskip = lux_levels[fLuxury];
514  }
515 
516  for( i = 0;(i != 24)&&(*seedptr != 0);i++) {
517  int_seed_table[i] = *seedptr % fIntModulus;
518  seedptr++;
519  }
520 
521  if(i != 24){
522  next_seed = int_seed_table[i-1];
523  for(;i != 24;i++) {
524  k_multiple = next_seed / ecuyer_a;
525  next_seed = ecuyer_b * (next_seed - k_multiple * ecuyer_a)
526  - k_multiple * ecuyer_c ;
527  if(next_seed < 0)next_seed += ecuyer_d;
528  int_seed_table[i] = next_seed % fIntModulus;
529  }
530  }
531 
532  for(i = 0;i != 24;i++)
533  fFloatSeedTable[i] = int_seed_table[i] * fMantissaBit24;
534 
535  fIlag = 23;
536  fJlag = 9;
537  fCarry = 0. ;
538 
539  if( fFloatSeedTable[23] == 0. ) fCarry = fMantissaBit24;
540 
541  fCount24 = 0;
542 }
543 
544 ////////////////////////////////////////////////////////////////////////////////
545 /// The initialisation is carried out using a Multiplicative
546 /// Congruential generator using formula constants of L'Ecuyer
547 /// as described in "A review of pseudorandom number generators"
548 /// (Fred James) published in Computer Physics Communications 60 (1990)
549 /// pages 329-344
550 ///
551 /// modified for the case of seed = 0. In that case a random 64 bits seed based on
552 /// TUUID (using TRandom3(0) ) is generated in order to have a unique seed
553 ///
554 
555 void TRandom1::SetSeed2(UInt_t seed, int lux)
556 {
557  const int ecuyer_a = 53668;
558  const int ecuyer_b = 40014;
559  const int ecuyer_c = 12211;
560  const int ecuyer_d = 2147483563;
561 
562  const int lux_levels[5] = {0,24,73,199,365};
563 
564  UInt_t int_seed_table[24];
565 
566  // case of seed == 0
567  // use a random seed based on TRandom2(0) which is based on the UUID
568  if (seed == 0) {
569  UInt_t randSeeds[25];
570  TRandom3 r2(0);
571  for (int j = 0; j < 24; ++j)
572  randSeeds[j] = static_cast<UInt_t> (4294967296.*r2.Rndm());
573  randSeeds[24] = 0;
574  SetSeeds(randSeeds, lux);
575  return;
576  }
577 
578 
579  Long64_t next_seed = seed;
580  Long64_t k_multiple;
581  int i;
582 
583  // number of additional random numbers that need to be 'thrown away'
584  // every 24 numbers is set using fLuxury level variable.
585 
586  fSeed = seed;
587  if( (lux > 4)||(lux < 0) ) {
588  if(lux >= 24) {
589  fNskip = lux - 24;
590  } else {
591  fNskip = lux_levels[3]; // corresponds to default fLuxury level
592  }
593  } else {
594  fLuxury = lux;
595  fNskip = lux_levels[fLuxury];
596  }
597 
598 
599  for(i = 0;i != 24;i++) {
600  k_multiple = next_seed / ecuyer_a;
601  next_seed = ecuyer_b * (next_seed - k_multiple * ecuyer_a)
602  - k_multiple * ecuyer_c ;
603  if(next_seed < 0)next_seed += ecuyer_d;
604  int_seed_table[i] = next_seed % fIntModulus;
605  }
606 
607  for(i = 0;i != 24;i++)
608  fFloatSeedTable[i] = int_seed_table[i] * fMantissaBit24;
609 
610  fIlag = 23;
611  fJlag = 9;
612  fCarry = 0. ;
613 
614  if( fFloatSeedTable[23] == 0. ) fCarry = fMantissaBit24;
615 
616  fCount24 = 0;
617 }
618 
620 {
621  // Set RanLux seed using default luxury level
622  SetSeed2(seed);
623 }
const UInt_t * fTheSeeds
Definition: TRandom1.h:42
Random number generator class based on M.
Definition: TRandom3.h:29
long long Long64_t
Definition: RtypesCore.h:69
Int_t fIlag
Definition: TRandom1.h:34
virtual void SetSeed2(UInt_t seed, Int_t lux=3)
The initialisation is carried out using a Multiplicative Congruential generator using formula constan...
Definition: TRandom1.cxx:555
float Float_t
Definition: RtypesCore.h:53
virtual void RndmArray(Int_t size, Float_t *vect)
return an array of random numbers in ]0,1]
Definition: TRandom1.cxx:413
const UInt_t fgSeedTable[215][2]
Definition: TRandom1.cxx:33
Float_t fCarry
Definition: TRandom1.h:38
const Int_t fIntModulus
Definition: TRandom1.h:39
virtual void SetSeeds(const UInt_t *seeds, Int_t lux=3)
set seeds
Definition: TRandom1.cxx:478
int Int_t
Definition: RtypesCore.h:41
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:501
The Ranlux Random number generator class.
Definition: TRandom1.h:29
virtual void SetSeed(UInt_t seed)
Set the random generator seed.
Definition: TRandom1.cxx:619
TRandom1()
default constructor
Definition: TRandom1.cxx:287
virtual ~TRandom1()
destructor
Definition: TRandom1.cxx:341
static Vc_ALWAYS_INLINE Vector< T > abs(const Vector< T > &x)
Definition: vector.h:450
virtual Double_t Rndm(Int_t i=0)
Machine independent random number generator.
Definition: TRandom3.cxx:93
Int_t fCount24
Definition: TRandom1.h:36
Int_t fJlag
Definition: TRandom1.h:35
const Double_t fMantissaBit12
Definition: TRandom1.h:44
Int_t fLuxury
Definition: TRandom1.h:33
ClassImp(TRandom1) TRandom1
Luxury level is set in the same way as the original FORTRAN routine.
Definition: TRandom1.cxx:251
unsigned int UInt_t
Definition: RtypesCore.h:42
const Double_t fMantissaBit24
Definition: TRandom1.h:43
Float_t fFloatSeedTable[24]
Definition: TRandom1.h:37
UInt_t fSeed
Definition: TRandom.h:32
double Double_t
Definition: RtypesCore.h:55
static Int_t fgNumEngines
Definition: TRandom1.h:40
static void GetTableSeeds(UInt_t *seeds, Int_t index)
static function returning the table of seeds
Definition: TRandom1.cxx:348
static Int_t fgMaxIndex
Definition: TRandom1.h:41
unsigned int r2[N_CITIES]
Definition: simanTSP.cxx:322
Int_t fNskip
Definition: TRandom1.h:32
virtual Double_t Rndm(Int_t i=0)
return a random number in ]0,1]
Definition: TRandom1.cxx:360