Logo ROOT  
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 
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;
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;
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;
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;
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 /// Set the state of the generator providing an array of seeds
476 /// The array of seeds can be of size 24 or less. In case of an array of n seeds with n < 24
477 /// the n+1 element must be equal to zero.
478 /// The other elements are the initialized using a Multiplicative
479 /// Congruential generator using formula constants of L'Ecuyer
480 /// as described in "A review of pseudorandom number generators"
481 /// (Fred James) published in Computer Physics Communications 60 (1990)
482 /// pages 329-344
483 ///
484 /// \param[in] seeds array of seeds (size from 1 to 24)
485 /// \param[in] lux Luxury level
486 ///
487 void TRandom1::SetSeeds(const UInt_t *seeds, int lux)
488 {
489  const int ecuyer_a = 53668;
490  const int ecuyer_b = 40014;
491  const int ecuyer_c = 12211;
492  const int ecuyer_d = 2147483563;
493 
494  const int lux_levels[5] = {0,24,73,199,365};
495  int i;
496  UInt_t int_seed_table[24];
497  Long64_t k_multiple,next_seed;
498  const UInt_t *seedptr;
499 
500  fTheSeeds = seeds;
501  seedptr = seeds;
502 
503  if(seeds == 0) {
504  SetSeed2(fSeed,lux);
505  fTheSeeds = &fSeed;
506  return;
507  }
508 
509  fSeed = *seeds;
510 
511 // number of additional random numbers that need to be 'thrown away'
512 // every 24 numbers is set using fLuxury level variable.
513 
514  if( (lux > 4)||(lux < 0) ) {
515  if(lux >= 24) {
516  fNskip = lux - 24;
517  } else {
518  fNskip = lux_levels[3]; // corresponds to default fLuxury level
519  }
520  } else {
521  fLuxury = lux;
522  fNskip = lux_levels[fLuxury];
523  }
524 
525  for( i = 0;(i != 24)&&(*seedptr != 0);i++) {
526  int_seed_table[i] = *seedptr % fIntModulus;
527  seedptr++;
528  }
529 
530  if(i != 24){
531  next_seed = int_seed_table[i-1];
532  for(;i != 24;i++) {
533  k_multiple = next_seed / ecuyer_a;
534  next_seed = ecuyer_b * (next_seed - k_multiple * ecuyer_a)
535  - k_multiple * ecuyer_c ;
536  if(next_seed < 0)next_seed += ecuyer_d;
537  int_seed_table[i] = next_seed % fIntModulus;
538  }
539  }
540 
541  for(i = 0;i != 24;i++)
542  fFloatSeedTable[i] = int_seed_table[i] * fMantissaBit24;
543 
544  fIlag = 23;
545  fJlag = 9;
546  fCarry = 0. ;
547 
548  if( fFloatSeedTable[23] == 0. ) fCarry = fMantissaBit24;
549 
550  fCount24 = 0;
551 }
552 
553 ////////////////////////////////////////////////////////////////////////////////
554 /// Set the state of the generator providing a single seed value and a
555 /// luxury level.
556 /// The initialisation of the other state values is carried out using a Multiplicative
557 /// Congruential generator using formula constants of L'Ecuyer
558 /// as described in "A review of pseudorandom number generators"
559 /// (Fred James) published in Computer Physics Communications 60 (1990)
560 /// pages 329-344
561 ///
562 /// Note: When the provided seed = 0, a random and unique seed is generated
563 ///
564 /// \param[in] seed seed value (note special case if seed=0)
565 /// \param[in] lux Luxury level
566 ///
567 ///
569 {
570  // case of seed == 0
571  // use a random seed based on TRandom3(0) which is based on the UUID
572  if (seed == 0) {
573  UInt_t randSeeds[25];
574  TRandom3 r2(0);
575  for (int j = 0; j < 24; ++j)
576  randSeeds[j] = static_cast<UInt_t>(4294967296. * r2.Rndm());
577  randSeeds[24] = 0;
578  SetSeeds(randSeeds, lux);
579  return;
580  }
581 
582  // seed List requires a zero terminated array in case the size of the array is less that the
583  // state size (24)
584  UInt_t seedList[2] = {seed, 0};
585  SetSeeds(seedList, lux);
586 }
587 
588 ////////////////////////////////////////////////////////////////////////////////
589 /// Set the state of the generator providing a single seed value and using
590 /// the luxury level defined when constructing the class
591 /// The initialisation of the other state values is carried out using a Multiplicative
592 /// Congruential generator.
593 /// Note: When seed = 0, a random and unique seed is generated
594 ///
595 /// \param[in] seed seed value (note special case if seed=0)
596 ///
597 
599 {
600  // Set RanLux seed using the luxury level provided in the constructor
601  SetSeed2(seed,fLuxury);
602 }
TRandom1::fNskip
Int_t fNskip
Definition: TRandom1.h:36
TRandom1::RndmArray
virtual void RndmArray(Int_t size, Float_t *vect)
return an array of random numbers in ]0,1]
Definition: TRandom1.cxx:413
TRandom1::fCarry
Float_t fCarry
Definition: TRandom1.h:42
TRandom1::fLuxury
Int_t fLuxury
Definition: TRandom1.h:37
TRandom1::Rndm
virtual Double_t Rndm()
return a random number in ]0,1]
Definition: TRandom1.cxx:360
ClassImp
#define ClassImp(name)
Definition: Rtypes.h:364
TRandom1::SetSeeds
virtual void SetSeeds(const UInt_t *seeds, Int_t lux=3)
Set the state of the generator providing an array of seeds The array of seeds can be of size 24 or le...
Definition: TRandom1.cxx:487
Long64_t
long long Long64_t
Definition: RtypesCore.h:73
TRandom1::~TRandom1
virtual ~TRandom1()
destructor
Definition: TRandom1.cxx:341
TGeant4Unit::lux
static constexpr double lux
Definition: TGeant4SystemOfUnits.h:326
Float_t
float Float_t
Definition: RtypesCore.h:57
TRandom::fSeed
UInt_t fSeed
Definition: TRandom.h:36
TRandom1::fMantissaBit12
const Double_t fMantissaBit12
Definition: TRandom1.h:48
TRandom1::fMantissaBit24
const Double_t fMantissaBit24
Definition: TRandom1.h:47
TRandom1::fTheSeeds
const UInt_t * fTheSeeds
Definition: TRandom1.h:46
TRandom1::TRandom1
TRandom1()
default constructor
Definition: TRandom1.cxx:287
TRandom1.h
TRandom1::fCount24
Int_t fCount24
Definition: TRandom1.h:40
TRandom3
Definition: TRandom3.h:27
TMath::Power
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:737
TRandom3.h
UInt_t
unsigned int UInt_t
Definition: RtypesCore.h:46
TRandom1::fgNumEngines
static Int_t fgNumEngines
Definition: TRandom1.h:44
double
double
Definition: Converters.cxx:921
ULong_t
unsigned long ULong_t
Definition: RtypesCore.h:55
TRandom1::fFloatSeedTable
Float_t fFloatSeedTable[24]
Definition: TRandom1.h:41
unsigned int
Double_t
double Double_t
Definition: RtypesCore.h:59
TRandom1::fIlag
Int_t fIlag
Definition: TRandom1.h:38
TRandom3::Rndm
virtual Double_t Rndm()
Machine independent random number generator.
Definition: TRandom3.cxx:99
fgSeedTable
const UInt_t fgSeedTable[215][2]
Definition: TRandom1.cxx:33
TRandom1::fJlag
Int_t fJlag
Definition: TRandom1.h:39
TRandom1
Definition: TRandom1.h:27
TRandom1::GetTableSeeds
static void GetTableSeeds(UInt_t *seeds, Int_t index)
static function returning the table of seeds
Definition: TRandom1.cxx:348
TRandom1::SetSeed
virtual void SetSeed(ULong_t seed)
Set the state of the generator providing a single seed value and using the luxury level defined when ...
Definition: TRandom1.cxx:598
TRandom1::fgMaxIndex
static Int_t fgMaxIndex
Definition: TRandom1.h:45
TRandom1::fIntModulus
const Int_t fIntModulus
Definition: TRandom1.h:43
TRandom1::SetSeed2
virtual void SetSeed2(UInt_t seed, Int_t lux=3)
Set the state of the generator providing a single seed value and a luxury level.
Definition: TRandom1.cxx:568
TMath
TMath.
Definition: TMathBase.h:35
TMath.h
int