StRoot: StarClassLibrary/RanluxEngine.cc Source File

00001 /***************************************************************************
00002  *
00003  * $Id: RanluxEngine.cc,v 1.3 2003/09/02 17:59:34 perev Exp $
00004  *
00005  * Author: Original code from CLHEP by G. Cosmo
00006  *         modified for SCL bl
00007  ***************************************************************************
00008  *
00009  * Description:
00010  * -----------------------------------------------------------------------
00011  *                             HEP Random
00012  *                        --- RanluxEngine ---
00013  *                      class implementation file
00014  * -----------------------------------------------------------------------
00015  * This file is part of Geant4 (simulation toolkit for HEP).
00016  *
00017  * Ranlux random number generator originally implemented in FORTRAN77
00018  * by Fred James as part of the MATHLIB HEP library.
00019  * 'RanluxEngine' is designed to fit into the CLHEP random number
00020  * class structure.
00021  *
00022  ***************************************************************************
00023  *
00024  * $Log: RanluxEngine.cc,v $
00025  * Revision 1.3  2003/09/02 17:59:34  perev
00026  * gcc 3.2 updates + WarnOff
00027  *
00028  * Revision 1.2  1999/12/07 23:43:04  ullrich
00029  * Modified to get rid of warnings on Linux.
00030  *
00031  * Revision 1.1  1999/01/30 03:59:01  fisyak
00032  * Root Version of StarClassLibrary
00033  *
00034  * Revision 1.1  1999/01/23 00:29:12  ullrich
00035  * Initial Revision
00036  *
00037  **************************************************************************/
00038 #include "RanluxEngine.h"
00039 
00040 RanluxEngine::RanluxEngine(long seed, HepInt lux)
00041 : int_modulus(0x1000000),
00042   mantissa_bit_24((HepFloat) ::pow(0.5,24.)),
00043   mantissa_bit_12((HepFloat) ::pow(0.5,12.))
00044 {
00045    luxury = lux;
00046    setSeed(seed, luxury);
00047    setSeeds(&theSeed, luxury);
00048 }
00049 
00050 RanluxEngine::~RanluxEngine() {}
00051 
00052 RanluxEngine::RanluxEngine(const RanluxEngine &p)
00053 : int_modulus(0x1000000),
00054   mantissa_bit_24((HepFloat) ::pow(0.5,24.)),
00055   mantissa_bit_12((HepFloat) ::pow(0.5,12.))
00056 {
00057   if ((this != &p) && (&p)) {
00058     theSeed = p.getSeed();
00059     setSeeds(&theSeed, p.luxury);
00060     for (HepInt i=0; i<24; ++i)
00061       float_seed_table[i] = p.float_seed_table[i];
00062     nskip = p.nskip;
00063     luxury = p.luxury;
00064     i_lag = p.i_lag;  j_lag = p.j_lag;
00065     carry = p.carry;
00066     count24 = p.count24;
00067   }
00068 }
00069 
00070 RanluxEngine & RanluxEngine::operator = (const RanluxEngine &p)
00071 {
00072   if ((this != &p) && (&p)) {
00073     theSeed = p.getSeed();
00074     setSeeds(&theSeed, p.luxury);
00075     for (HepInt i=0; i<24; ++i)
00076       float_seed_table[i] = p.float_seed_table[i];
00077     nskip = p.nskip;
00078     luxury = p.luxury;
00079     i_lag = p.i_lag;  j_lag = p.j_lag;
00080     carry = p.carry;
00081     count24 = p.count24;
00082   }
00083   return *this;
00084 }
00085 
00086 void RanluxEngine::setSeed(long seed, HepInt lux) {
00087 
00088 // The initialisation is carried out using a Multiplicative
00089 // Congruential generator using formula constants of L'Ecuyer 
00090 // as described in "A review of pseudorandom number generators"
00091 // (Fred James) published in Computer Physics Communications 60 (1990)
00092 // pages 329-344
00093 
00094   const HepInt ecuyer_a = 53668;
00095   const HepInt ecuyer_b = 40014;
00096   const HepInt ecuyer_c = 12211;
00097   const HepInt ecuyer_d = 2147483563;
00098 
00099   const HepInt lux_levels[5] = {0,24,73,199,365};  
00100 
00101   long int_seed_table[24];
00102   long next_seed = seed;
00103   long k_multiple;
00104   register HepInt i;
00105   
00106 // number of additional random numbers that need to be 'thrown away'
00107 // every 24 numbers is set using luxury level variable.
00108 
00109   theSeed = seed;
00110   if( (lux > 4)||(lux < 0) ){
00111      if(lux >= 24){
00112         nskip = lux - 24;
00113      }else{
00114         nskip = lux_levels[3]; // corresponds to default luxury level
00115      }
00116   }else{
00117      luxury = lux;
00118      nskip = lux_levels[luxury];
00119   }
00120 
00121    
00122   for(i = 0;i != 24;i++){
00123      k_multiple = next_seed / ecuyer_a;
00124      next_seed = ecuyer_b * (next_seed - k_multiple * ecuyer_a) 
00125      - k_multiple * ecuyer_c ;
00126      if(next_seed < 0)next_seed += ecuyer_d;
00127      int_seed_table[i] = next_seed % int_modulus;
00128   }     
00129 
00130   for(i = 0;i != 24;i++)
00131     float_seed_table[i] = int_seed_table[i] * mantissa_bit_24;
00132 
00133   i_lag = 23;
00134   j_lag = 9;
00135   carry = 0. ;
00136 
00137   if( float_seed_table[23] == 0. ) carry = mantissa_bit_24;
00138   
00139   count24 = 0;
00140 }
00141 
00142 void RanluxEngine::setSeeds(const long *seeds, HepInt lux) {
00143 
00144    const HepInt ecuyer_a = 53668;
00145    const HepInt ecuyer_b = 40014;
00146    const HepInt ecuyer_c = 12211;
00147    const HepInt ecuyer_d = 2147483563;
00148 
00149    const HepInt lux_levels[5] = {0,24,73,199,365};
00150    HepInt i;
00151    long int_seed_table[24];
00152    long k_multiple,next_seed;
00153    const long *seedptr; 
00154 
00155    theSeeds = seeds;
00156    seedptr  = seeds;
00157  
00158    if(seeds == 0){
00159       setSeed(theSeed,lux);
00160       theSeeds = &theSeed;
00161       return;
00162    }
00163 
00164    theSeed = *seeds;
00165 
00166 // number of additional random numbers that need to be 'thrown away'
00167 // every 24 numbers is set using luxury level variable.
00168 
00169   if( (lux > 4)||(lux < 0) ){
00170      if(lux >= 24){
00171         nskip = lux - 24;
00172      }else{
00173         nskip = lux_levels[3]; // corresponds to default luxury level
00174      }
00175   }else{
00176      luxury = lux;
00177      nskip = lux_levels[luxury];
00178   }
00179       
00180   for( i = 0;(i != 24)&&(*seedptr != 0);i++){
00181       int_seed_table[i] = *seedptr % int_modulus;
00182       seedptr++;
00183   }                    
00184 
00185   if(i != 24){
00186      next_seed = int_seed_table[i-1];
00187      for(;i != 24;i++){
00188         k_multiple = next_seed / ecuyer_a;
00189         next_seed = ecuyer_b * (next_seed - k_multiple * ecuyer_a) 
00190         - k_multiple * ecuyer_c ;
00191         if(next_seed < 0)next_seed += ecuyer_d;
00192         int_seed_table[i] = next_seed % int_modulus;
00193      }          
00194   }
00195 
00196   for(i = 0;i != 24;i++)
00197     float_seed_table[i] = int_seed_table[i] * mantissa_bit_24;
00198 
00199   i_lag = 23;
00200   j_lag = 9;
00201   carry = 0. ;
00202 
00203   if( float_seed_table[23] == 0. ) carry = mantissa_bit_24;
00204   
00205   count24 = 0;
00206 }
00207 
00208 void RanluxEngine::saveStatus() const
00209 {
00210    ofstream outFile("Ranlux.conf", ios::out ) ;
00211 
00212    if (!outFile.bad()) {
00213      outFile << theSeed << endl;
00214      for (HepInt i=0; i<24; ++i)
00215        outFile << float_seed_table[i] << " ";
00216      outFile << endl;
00217      outFile << i_lag << " " << j_lag << endl;
00218      outFile << carry << " " << count24 << endl;
00219    }
00220 }
00221 
00222 void RanluxEngine::restoreStatus()
00223 {
00224    ifstream inFile("Ranlux.conf", ios::in);
00225 
00226    if (!inFile.bad() && !inFile.eof()) {
00227      inFile >> theSeed;
00228      for (HepInt i=0; i<24; ++i)
00229        inFile >> float_seed_table[i];
00230      inFile >> i_lag; inFile >> j_lag;
00231      inFile >> carry; inFile >> count24;
00232    }
00233 }
00234 
00235 void RanluxEngine::showStatus() const
00236 {
00237    cout << endl;
00238    cout << "--------- Ranlux engine status ---------" << endl;
00239    cout << " Initial seed = " << theSeed << endl;
00240    cout << " float_seed_table[] = ";
00241    for (HepInt i=0; i<24; ++i)
00242      cout << float_seed_table[i] << " ";
00243    cout << endl;
00244    cout << " i_lag = " << i_lag << ", j_lag = " << j_lag << endl;
00245    cout << " carry = " << carry << ", count24 = " << count24 << endl;
00246    cout << "----------------------------------------" << endl;
00247 }
00248 
00249 HepDouble RanluxEngine::flat() {
00250 
00251   HepFloat next_random;
00252   register HepFloat uni;
00253   register HepInt i;
00254 
00255   uni = float_seed_table[j_lag] - float_seed_table[i_lag] - carry;
00256   if(uni < 0. ){
00257      uni += 1.0;
00258      carry = mantissa_bit_24;
00259   }else{
00260      carry = 0.;
00261   }
00262 
00263   float_seed_table[i_lag] = uni;
00264   i_lag --;
00265   j_lag --;
00266   if(i_lag < 0) i_lag = 23;
00267   if(j_lag < 0) j_lag = 23;
00268 
00269   if( uni < mantissa_bit_12 ){
00270      uni += mantissa_bit_24 * float_seed_table[j_lag];
00271      if( uni == 0) uni = mantissa_bit_24 * mantissa_bit_24;
00272   }
00273   next_random = uni;
00274   count24 ++;
00275 
00276 // every 24th number generation, several random numbers are generated
00277 // and wasted depending upon the luxury level.
00278 
00279   if(count24 == 24 ){
00280      count24 = 0;
00281      for( i = 0; i != nskip ; i++){
00282          uni = float_seed_table[j_lag] - float_seed_table[i_lag] - carry;
00283          if(uni < 0. ){
00284             uni += 1.0;
00285             carry = mantissa_bit_24;
00286          }else{
00287             carry = 0.;
00288          }
00289          float_seed_table[i_lag] = uni;
00290          i_lag --;
00291          j_lag --;
00292          if(i_lag < 0)i_lag = 23;
00293          if(j_lag < 0) j_lag = 23;
00294       }
00295   } 
00296   return (HepDouble) next_random;
00297 }
00298 
00299 void RanluxEngine::flatArray(const HepInt size, HepDouble* vect)
00300 {
00301   HepFloat next_random;
00302   register HepFloat uni;
00303   register HepInt i;
00304   register HepInt index;
00305 
00306   for (index=0; index<size; ++index) {
00307     uni = float_seed_table[j_lag] - float_seed_table[i_lag] - carry;
00308     if(uni < 0. ){
00309        uni += 1.0;
00310        carry = mantissa_bit_24;
00311     }else{
00312        carry = 0.;
00313     }
00314 
00315     float_seed_table[i_lag] = uni;
00316     i_lag --;
00317     j_lag --;
00318     if(i_lag < 0) i_lag = 23;
00319     if(j_lag < 0) j_lag = 23;
00320 
00321     if( uni < mantissa_bit_12 ){
00322        uni += mantissa_bit_24 * float_seed_table[j_lag];
00323        if( uni == 0) uni = mantissa_bit_24 * mantissa_bit_24;
00324     }
00325     next_random = uni;
00326     vect[index] = (HepDouble)next_random;
00327     count24 ++;
00328 
00329 // every 24th number generation, several random numbers are generated
00330 // and wasted depending upon the luxury level.
00331 
00332     if(count24 == 24 ){
00333        count24 = 0;
00334        for( i = 0; i != nskip ; i++){
00335            uni = float_seed_table[j_lag] - float_seed_table[i_lag] - carry;
00336            if(uni < 0. ){
00337               uni += 1.0;
00338               carry = mantissa_bit_24;
00339            }else{
00340               carry = 0.;
00341            }
00342            float_seed_table[i_lag] = uni;
00343            i_lag --;
00344            j_lag --;
00345            if(i_lag < 0)i_lag = 23;
00346            if(j_lag < 0) j_lag = 23;
00347         }
00348     }
00349   }
00350 } 
00351 
00352 void
00353 #ifndef ST_NO_TEMPLATE_DEF_ARGS
00354 RanluxEngine::flatArray(vector<HepDouble>& vec)
00355 #else
00356 RanluxEngine::flatArray(vector<HepDouble,allocator<HepDouble> >& vec)
00357 #endif
00358 {
00359   HepFloat next_random;
00360   register HepFloat uni;
00361   register HepInt i;
00362 
00363   for (unsigned int index=0; index<vec.size(); ++index) {
00364     uni = float_seed_table[j_lag] - float_seed_table[i_lag] - carry;
00365     if(uni < 0. ){
00366        uni += 1.0;
00367        carry = mantissa_bit_24;
00368     }else{
00369        carry = 0.;
00370     }
00371 
00372     float_seed_table[i_lag] = uni;
00373     i_lag --;
00374     j_lag --;
00375     if(i_lag < 0) i_lag = 23;
00376     if(j_lag < 0) j_lag = 23;
00377 
00378     if( uni < mantissa_bit_12 ){
00379        uni += mantissa_bit_24 * float_seed_table[j_lag];
00380        if( uni == 0) uni = mantissa_bit_24 * mantissa_bit_24;
00381     }
00382     next_random = uni;
00383     vec[index] = (HepDouble)next_random;
00384     count24 ++;
00385 
00386 // every 24th number generation, several random numbers are generated
00387 // and wasted depending upon the luxury level.
00388 
00389     if(count24 == 24 ){
00390        count24 = 0;
00391        for( i = 0; i != nskip ; i++){
00392            uni = float_seed_table[j_lag] - float_seed_table[i_lag] - carry;
00393            if(uni < 0. ){
00394               uni += 1.0;
00395               carry = mantissa_bit_24;
00396            }else{
00397               carry = 0.;
00398            }
00399            float_seed_table[i_lag] = uni;
00400            i_lag --;
00401            j_lag --;
00402            if(i_lag < 0)i_lag = 23;
00403            if(j_lag < 0) j_lag = 23;
00404         }
00405     }
00406   }
00407 }