// @(#)root/mathcore:$Id$ // Author: Rene Brun from CLHEP & CERNLIB 04/05/2006 /** \class TRandom1 The Ranlux Random number generator class The algorithm for this random engine has been taken from the original implementation in FORTRAN by Fred James as part of CLHEP. The initialisation is carried out using a Multiplicative Congruential generator using formula constants of L'Ecuyer as described in "F.James, Comp. Phys. Comm. 60 (1990) 329-344". @ingroup Random */ #include "TRandom1.h" #include "TRandom3.h" #include "TMath.h" #include // Number of instances with automatic seed selection int TRandom1::fgNumEngines = 0; // Maximum index into the seed table int TRandom1::fgMaxIndex = 215; #ifndef __CINT__ const UInt_t fgSeedTable[215][2] = { { 9876, 54321 }, { 1299961164, 253987020 }, { 669708517, 2079157264 }, { 190904760, 417696270 }, { 1289741558, 1376336092 }, { 1803730167, 324952955 }, { 489854550, 582847132 }, { 1348037628, 1661577989 }, { 350557787, 1155446919 }, { 591502945, 634133404 }, { 1901084678, 862916278 }, { 1988640932, 1785523494 }, { 1873836227, 508007031 }, { 1146416592, 967585720 }, { 1837193353, 1522927634 }, { 38219936, 921609208 }, { 349152748, 112892610 }, { 744459040, 1735807920 }, { 1983990104, 728277902 }, { 309164507, 2126677523 }, { 362993787, 1897782044 }, { 556776976, 462072869 }, { 1584900822, 2019394912 }, { 1249892722, 791083656 }, { 1686600998, 1983731097 }, { 1127381380, 198976625 }, { 1999420861, 1810452455 }, { 1972906041, 664182577 }, { 84636481, 1291886301 }, { 1186362995, 954388413 }, { 2141621785, 61738584 }, { 1969581251, 1557880415 }, { 1150606439, 136325185 }, { 95187861, 1592224108 }, { 940517655, 1629971798 }, { 215350428, 922659102 }, { 786161212, 1121345074 }, { 1450830056, 1922787776 }, { 1696578057, 2025150487 }, { 1803414346, 1851324780 }, { 1017898585, 1452594263 }, { 1184497978, 82122239 }, { 633338765, 1829684974 }, { 430889421, 230039326 }, { 492544653, 76320266 }, { 389386975, 1314148944 }, { 1720322786, 709120323 }, { 1868768216, 1992898523 }, { 443210610, 811117710 }, { 1191938868, 1548484733 }, { 616890172, 159787986 }, { 935835339, 1231440405 }, { 1058009367, 1527613300 }, { 1463148129, 1970575097 }, { 1795336935, 434768675 }, { 274019517, 605098487 }, { 483689317, 217146977 }, { 2070804364, 340596558 }, { 930226308, 1602100969 }, { 989324440, 801809442 }, { 410606853, 1893139948 }, { 1583588576, 1219225407 }, { 2102034391, 1394921405 }, { 2005037790, 2031006861 }, { 1244218766, 923231061 }, { 49312790, 775496649 }, { 721012176, 321339902 }, { 1719909107, 1865748178 }, { 1156177430, 1257110891 }, { 307561322, 1918244397 }, { 906041433, 360476981 }, { 1591375755, 268492659 }, { 461522398, 227343256 }, { 2145930725, 2020665454 }, { 1938419274, 1331283701 }, { 174405412, 524140103 }, { 494343653, 18063908 }, { 1025534808, 181709577 }, { 2048959776, 1913665637 }, { 950636517, 794796256 }, { 1828843197, 1335757744 }, { 211109723, 983900607 }, { 825474095, 1046009991 }, { 374915657, 381856628 }, { 1241296328, 698149463 }, { 1260624655, 1024538273 }, { 900676210, 1628865823 }, { 697951025, 500570753 }, { 1007920268, 1708398558 }, { 264596520, 624727803 }, { 1977924811, 674673241 }, { 1440257718, 271184151 }, { 1928778847, 993535203 }, { 1307807366, 1801502463 }, { 1498732610, 300876954 }, { 1617712402, 1574250679 }, { 1261800762, 1556667280 }, { 949929273, 560721070 }, { 1766170474, 1953522912 }, { 1849939248, 19435166 }, { 887262858, 1219627824 }, { 483086133, 603728993 }, { 1330541052, 1582596025 }, { 1850591475, 723593133 }, { 1431775678, 1558439000 }, { 922493739, 1356554404 }, { 1058517206, 948567762 }, { 709067283, 1350890215 }, { 1044787723, 2144304941 }, { 999707003, 513837520 }, { 2140038663, 1850568788 }, { 1803100150, 127574047 }, { 867445693, 1149173981 }, { 408583729, 914837991 }, { 1166715497, 602315845 }, { 430738528, 1743308384 }, { 1388022681, 1760110496 }, { 1664028066, 654300326 }, { 1767741172, 1338181197 }, { 1625723550, 1742482745 }, { 464486085, 1507852127 }, { 754082421, 1187454014 }, { 1315342834, 425995190 }, { 960416608, 2004255418 }, { 1262630671, 671761697 }, { 59809238, 103525918 }, { 1205644919, 2107823293 }, { 1615183160, 1152411412 }, { 1024474681, 2118672937 }, { 1703877649, 1235091369 }, { 1821417852, 1098463802 }, { 1738806466, 1529062843 }, { 620780646, 1654833544 }, { 1070174101, 795158254 }, { 658537995, 1693620426 }, { 2055317555, 508053916 }, { 1647371686, 1282395762 }, { 29067379, 409683067 }, { 1763495989, 1917939635 }, { 1602690753, 810926582 }, { 885787576, 513818500 }, { 1853512561, 1195205756 }, { 1798585498, 1970460256 }, { 1819261032, 1306536501 }, { 1133245275, 37901 }, { 689459799, 1334389069 }, { 1730609912, 1854586207 }, { 1556832175, 1228729041 }, { 251375753, 683687209 }, { 2083946182, 1763106152 }, { 2142981854, 1365385561 }, { 763711891, 1735754548 }, { 1581256466, 173689858 }, { 2121337132, 1247108250 }, { 1004003636, 891894307 }, { 569816524, 358675254 }, { 626626425, 116062841 }, { 632086003, 861268491 }, { 1008211580, 779404957 }, { 1134217766, 1766838261 }, { 1423829292, 1706666192 }, { 942037869, 1549358884 }, { 1959429535, 480779114 }, { 778311037, 1940360875 }, { 1531372185, 2009078158 }, { 241935492, 1050047003 }, { 272453504, 1870883868 }, { 390441332, 1057903098 }, { 1230238834, 1548117688 }, { 1242956379, 1217296445 }, { 515648357, 1675011378 }, { 364477932, 355212934 }, { 2096008713, 1570161804 }, { 1409752526, 214033983 }, { 1288158292, 1760636178 }, { 407562666, 1265144848 }, { 1071056491, 1582316946 }, { 1014143949, 911406955 }, { 203080461, 809380052 }, { 125647866, 1705464126 }, { 2015685843, 599230667 }, { 1425476020, 668203729 }, { 1673735652, 567931803 }, { 1714199325, 181737617 }, { 1389137652, 678147926 }, { 288547803, 435433694 }, { 200159281, 654399753 }, { 1580828223, 1298308945 }, { 1832286107, 169991953 }, { 182557704, 1046541065 }, { 1688025575, 1248944426 }, { 1508287706, 1220577001 }, { 36721212, 1377275347 }, { 1968679856, 1675229747 }, { 279109231, 1835333261 }, { 1358617667, 1416978076 }, { 740626186, 2103913602 }, { 1882655908, 251341858 }, { 648016670, 1459615287 }, { 780255321, 154906988 }, { 857296483, 203375965 }, { 1631676846, 681204578 }, { 1906971307, 1623728832 }, { 1541899600, 1168449797 }, { 1267051693, 1020078717 }, { 1998673940, 1298394942 }, { 1914117058, 1381290704 }, { 426068513, 1381618498 }, { 139365577, 1598767734 }, { 2129910384, 952266588 }, { 661788054, 19661356 }, { 1104640222, 240506063 }, { 356133630, 1676634527 }, { 242242374, 1863206182 }, { 957935844, 1490681416 }}; #endif ClassImp(TRandom1) //////////////////////////////////////////////////////////////////////////////// /// Luxury level is set in the same way as the original FORTRAN routine. /// level 0 (p=24): equivalent to the original RCARRY of Marsaglia /// and Zaman, very long period, but fails many tests. /// level 1 (p=48): considerable improvement in quality over level 0, /// now passes the gap test, but still fails spectral test. /// level 2 (p=97): passes all known tests, but theoretically still /// defective. /// level 3 (p=223): DEFAULT VALUE. Any theoretically possible /// correlations have very small chance of being observed. /// level 4 (p=389): highest possible luxury, all 24 bits chaotic. TRandom1::TRandom1(UInt_t seed, Int_t lux) : fIntModulus(0x1000000), fMantissaBit24( TMath::Power(0.5,24.) ), fMantissaBit12( TMath::Power(0.5,12.) ) { UInt_t seedlist[2]={0,0}; fTheSeeds = &fSeed; fLuxury = lux; SetSeed2(seed, fLuxury); // in case seed = 0 SetSeed2 calls already SetSeeds if (seed != 0) { // setSeeds() wants a zero terminated array! seedlist[0]=fSeed; seedlist[1]=0; SetSeeds(seedlist, fLuxury); } } //////////////////////////////////////////////////////////////////////////////// ///default constructor TRandom1::TRandom1() : fIntModulus(0x1000000), fMantissaBit24( TMath::Power(0.5,24.) ), fMantissaBit12( TMath::Power(0.5,12.) ) { fTheSeeds = &fSeed; UInt_t seed; UInt_t seedlist[2]={0,0}; fLuxury = 3; int cycle = abs(int(fgNumEngines/fgMaxIndex)); int curIndex = abs(int(fgNumEngines%fgMaxIndex)); fgNumEngines +=1; UInt_t mask = ((cycle & 0x007fffff) << 8); GetTableSeeds( seedlist, curIndex ); seed = seedlist[0]^mask; SetSeed2(seed, fLuxury); // setSeeds() wants a zero terminated array! seedlist[0]=fSeed; //<============= seedlist[1]=0; SetSeeds(seedlist, fLuxury); } //////////////////////////////////////////////////////////////////////////////// ///constructor TRandom1::TRandom1(int rowIndex, int colIndex, int lux) : fIntModulus(0x1000000), fMantissaBit24( TMath::Power(0.5,24.) ), fMantissaBit12( TMath::Power(0.5,12.) ) { fTheSeeds = &fSeed; UInt_t seed; UInt_t seedlist[2]={0,0}; fLuxury = lux; int cycle = abs(int(rowIndex/fgMaxIndex)); int row = abs(int(rowIndex%fgMaxIndex)); int col = abs(int(colIndex%2)); UInt_t mask = (( cycle & 0x000007ff ) << 20 ); GetTableSeeds( seedlist, row ); seed = ( seedlist[col] )^mask; SetSeed2(seed, fLuxury); // setSeeds() wants a zero terminated array! seedlist[0]=fSeed; seedlist[1]=0; SetSeeds(seedlist, fLuxury); } //////////////////////////////////////////////////////////////////////////////// ///destructor TRandom1::~TRandom1() { } //////////////////////////////////////////////////////////////////////////////// ///static function returning the table of seeds void TRandom1::GetTableSeeds(UInt_t* seeds, Int_t index) { if ((index >= 0) && (index < 215)) { seeds[0] = fgSeedTable[index][0]; seeds[1] = fgSeedTable[index][1]; } else seeds = 0; } //////////////////////////////////////////////////////////////////////////////// ///return a random number in ]0,1] Double_t TRandom1::Rndm(Int_t) { float next_random; float uni; int i; uni = fFloatSeedTable[fJlag] - fFloatSeedTable[fIlag] - fCarry; if(uni < 0. ) { uni += 1.0; fCarry = fMantissaBit24; } else { fCarry = 0.; } fFloatSeedTable[fIlag] = uni; fIlag --; fJlag --; if(fIlag < 0) fIlag = 23; if(fJlag < 0) fJlag = 23; if( uni < fMantissaBit12 ){ uni += fMantissaBit24 * fFloatSeedTable[fJlag]; if( uni == 0) uni = fMantissaBit24 * fMantissaBit24; } next_random = uni; fCount24 ++; // every 24th number generation, several random numbers are generated // and wasted depending upon the fLuxury level. if(fCount24 == 24 ) { fCount24 = 0; for( i = 0; i != fNskip ; i++) { uni = fFloatSeedTable[fJlag] - fFloatSeedTable[fIlag] - fCarry; if(uni < 0. ) { uni += 1.0; fCarry = fMantissaBit24; } else { fCarry = 0.; } fFloatSeedTable[fIlag] = uni; fIlag --; fJlag --; if(fIlag < 0)fIlag = 23; if(fJlag < 0) fJlag = 23; } } return (double) next_random; } //////////////////////////////////////////////////////////////////////////////// ///return an array of random numbers in ]0,1] void TRandom1::RndmArray(const Int_t size, Float_t *vect) { for (Int_t i=0;i 4)||(lux < 0) ) { if(lux >= 24) { fNskip = lux - 24; } else { fNskip = lux_levels[3]; // corresponds to default fLuxury level } } else { fLuxury = lux; fNskip = lux_levels[fLuxury]; } for( i = 0;(i != 24)&&(*seedptr != 0);i++) { int_seed_table[i] = *seedptr % fIntModulus; seedptr++; } if(i != 24){ next_seed = int_seed_table[i-1]; for(;i != 24;i++) { k_multiple = next_seed / ecuyer_a; next_seed = ecuyer_b * (next_seed - k_multiple * ecuyer_a) - k_multiple * ecuyer_c ; if(next_seed < 0)next_seed += ecuyer_d; int_seed_table[i] = next_seed % fIntModulus; } } for(i = 0;i != 24;i++) fFloatSeedTable[i] = int_seed_table[i] * fMantissaBit24; fIlag = 23; fJlag = 9; fCarry = 0. ; if( fFloatSeedTable[23] == 0. ) fCarry = fMantissaBit24; fCount24 = 0; } //////////////////////////////////////////////////////////////////////////////// /// The initialisation is carried out using a Multiplicative /// Congruential generator using formula constants of L'Ecuyer /// as described in "A review of pseudorandom number generators" /// (Fred James) published in Computer Physics Communications 60 (1990) /// pages 329-344 /// /// modified for the case of seed = 0. In that case a random 64 bits seed based on /// TUUID (using TRandom3(0) ) is generated in order to have a unique seed /// void TRandom1::SetSeed2(UInt_t seed, int lux) { const int ecuyer_a = 53668; const int ecuyer_b = 40014; const int ecuyer_c = 12211; const int ecuyer_d = 2147483563; const int lux_levels[5] = {0,24,73,199,365}; UInt_t int_seed_table[24]; // case of seed == 0 // use a random seed based on TRandom2(0) which is based on the UUID if (seed == 0) { UInt_t randSeeds[25]; TRandom3 r2(0); for (int j = 0; j < 24; ++j) randSeeds[j] = static_cast (4294967296.*r2.Rndm()); randSeeds[24] = 0; SetSeeds(randSeeds, lux); return; } Long64_t next_seed = seed; Long64_t k_multiple; int i; // number of additional random numbers that need to be 'thrown away' // every 24 numbers is set using fLuxury level variable. fSeed = seed; if( (lux > 4)||(lux < 0) ) { if(lux >= 24) { fNskip = lux - 24; } else { fNskip = lux_levels[3]; // corresponds to default fLuxury level } } else { fLuxury = lux; fNskip = lux_levels[fLuxury]; } for(i = 0;i != 24;i++) { k_multiple = next_seed / ecuyer_a; next_seed = ecuyer_b * (next_seed - k_multiple * ecuyer_a) - k_multiple * ecuyer_c ; if(next_seed < 0)next_seed += ecuyer_d; int_seed_table[i] = next_seed % fIntModulus; } for(i = 0;i != 24;i++) fFloatSeedTable[i] = int_seed_table[i] * fMantissaBit24; fIlag = 23; fJlag = 9; fCarry = 0. ; if( fFloatSeedTable[23] == 0. ) fCarry = fMantissaBit24; fCount24 = 0; } void TRandom1::SetSeed(UInt_t seed) { // Set RanLux seed using default luxury level SetSeed2(seed); }