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
8The Ranlux Random number generator class
9
10The algorithm for this random engine has been taken from the original
11implementation in FORTRAN by Fred James as part of CLHEP.
12
13The initialisation is carried out using a Multiplicative Congruential
14generator using formula constants of L'Ecuyer as described in "F.James,
15Comp. 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
31int TRandom1::fgMaxIndex = 215;
32#ifndef __CINT__
33const 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
314TRandom1::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
413void 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
421void 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////////////////////////////////////////////////////////////////////////////////
476///set seeds
477
478void 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) {
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
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 the luxury level provided in the constructor
622 SetSeed2(seed,fLuxury);
623}
int Int_t
Definition: RtypesCore.h:41
unsigned int UInt_t
Definition: RtypesCore.h:42
unsigned long ULong_t
Definition: RtypesCore.h:51
double Double_t
Definition: RtypesCore.h:55
long long Long64_t
Definition: RtypesCore.h:69
float Float_t
Definition: RtypesCore.h:53
#define ClassImp(name)
Definition: Rtypes.h:365
const UInt_t fgSeedTable[215][2]
Definition: TRandom1.cxx:33
The Ranlux Random number generator class.
Definition: TRandom1.h:27
const Double_t fMantissaBit12
Definition: TRandom1.h:42
Float_t fFloatSeedTable[24]
Definition: TRandom1.h:35
TRandom1()
default constructor
Definition: TRandom1.cxx:287
Int_t fIlag
Definition: TRandom1.h:32
Int_t fNskip
Definition: TRandom1.h:30
Int_t fCount24
Definition: TRandom1.h:34
const UInt_t * fTheSeeds
Definition: TRandom1.h:40
Int_t fJlag
Definition: TRandom1.h:33
virtual Double_t Rndm()
return a random number in ]0,1]
Definition: TRandom1.cxx:360
static Int_t fgMaxIndex
Definition: TRandom1.h:39
Float_t fCarry
Definition: TRandom1.h:36
const Double_t fMantissaBit24
Definition: TRandom1.h:41
virtual ~TRandom1()
destructor
Definition: TRandom1.cxx:341
const Int_t fIntModulus
Definition: TRandom1.h:37
virtual void SetSeeds(const UInt_t *seeds, Int_t lux=3)
set seeds
Definition: TRandom1.cxx:478
virtual void RndmArray(Int_t size, Float_t *vect)
return an array of random numbers in ]0,1]
Definition: TRandom1.cxx:413
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
static void GetTableSeeds(UInt_t *seeds, Int_t index)
static function returning the table of seeds
Definition: TRandom1.cxx:348
virtual void SetSeed(ULong_t seed)
Set the random generator seed.
Definition: TRandom1.cxx:619
static Int_t fgNumEngines
Definition: TRandom1.h:38
Int_t fLuxury
Definition: TRandom1.h:31
Random number generator class based on M.
Definition: TRandom3.h:27
virtual Double_t Rndm()
Machine independent random number generator.
Definition: TRandom3.cxx:100
UInt_t fSeed
Definition: TRandom.h:30
static constexpr double lux
TMath.
Definition: TMathBase.h:35
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:725