Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/externals/clhep/src/RanluxEngine.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /externals/clhep/src/RanluxEngine.cc (Version 11.3.0) and /externals/clhep/src/RanluxEngine.cc (Version 10.0)


                                                   >>   1 // $Id:$
  1 // -*- C++ -*-                                      2 // -*- C++ -*-
  2 //                                                  3 //
  3 // -------------------------------------------      4 // -----------------------------------------------------------------------
  4 //                             HEP Random           5 //                             HEP Random
  5 //                        --- RanluxEngine ---      6 //                        --- RanluxEngine ---
  6 //                      class implementation f      7 //                      class implementation file
  7 // -------------------------------------------      8 // -----------------------------------------------------------------------
  8 // This file is part of Geant4 (simulation too      9 // This file is part of Geant4 (simulation toolkit for HEP).
  9 //                                                 10 //
 10 // Ranlux random number generator originally i     11 // Ranlux random number generator originally implemented in FORTRAN77
 11 // by Fred James as part of the MATHLIB HEP li     12 // by Fred James as part of the MATHLIB HEP library.
 12 // 'RanluxEngine' is designed to fit into the      13 // 'RanluxEngine' is designed to fit into the CLHEP random number
 13 // class structure.                                14 // class structure.
 14                                                    15 
 15 // ===========================================     16 // ===============================================================
 16 // Adeyemi Adesanya - Created: 6th November 19     17 // Adeyemi Adesanya - Created: 6th November 1995
 17 // Gabriele Cosmo - Adapted & Revised: 22nd No     18 // Gabriele Cosmo - Adapted & Revised: 22nd November 1995
 18 // Adeyemi Adesanya - Added setSeeds() method:     19 // Adeyemi Adesanya - Added setSeeds() method: 2nd February 1996
 19 // Gabriele Cosmo - Added flatArray() method:      20 // Gabriele Cosmo - Added flatArray() method: 8th February 1996
 20 //                - Minor corrections: 31st Oc     21 //                - Minor corrections: 31st October 1996
 21 //                - Added methods for engine s     22 //                - Added methods for engine status: 19th November 1996
 22 //                - Fixed bug in setSeeds(): 1     23 //                - Fixed bug in setSeeds(): 15th September 1997
 23 // J.Marraffino   - Added stream operators and     24 // J.Marraffino   - Added stream operators and related constructor.
 24 //                  Added automatic seed selec     25 //                  Added automatic seed selection from seed table and
 25 //                  engine counter: 14th Feb 1     26 //                  engine counter: 14th Feb 1998
 26 //                - Fixed bug: setSeeds() requ     27 //                - Fixed bug: setSeeds() requires a zero terminated
 27 //                  array of seeds: 19th Feb 1     28 //                  array of seeds: 19th Feb 1998
 28 // Ken Smith      - Added conversion operators     29 // Ken Smith      - Added conversion operators:  6th Aug 1998
 29 // J. Marraffino  - Remove dependence on hepSt     30 // J. Marraffino  - Remove dependence on hepString class  13 May 1999
 30 // M. Fischler    - In restore, checkFile for      31 // M. Fischler    - In restore, checkFile for file not found    03 Dec 2004
 31 // M. Fischler    - Methods put, getfor instan     32 // M. Fischler    - Methods put, getfor instance save/restore       12/8/04    
 32 // M. Fischler    - split get() into tag valid     33 // M. Fischler    - split get() into tag validation and 
 33 //                  getState() for anonymous r     34 //                  getState() for anonymous restores           12/27/04    
 34 // M. Fischler    - put/get for vectors of ulo     35 // M. Fischler    - put/get for vectors of ulongs   3/14/05
 35 // M. Fischler    - State-saving using only in     36 // M. Fischler    - State-saving using only ints, for portability 4/12/05
 36 //                                                 37 //
 37 // ===========================================     38 // ===============================================================
 38                                                    39 
 39 #include "CLHEP/Random/Random.h"                   40 #include "CLHEP/Random/Random.h"
 40 #include "CLHEP/Random/RanluxEngine.h"             41 #include "CLHEP/Random/RanluxEngine.h"
 41 #include "CLHEP/Random/engineIDulong.h"            42 #include "CLHEP/Random/engineIDulong.h"
 42 #include "CLHEP/Utility/atomic_int.h"          << 
 43                                                << 
 44 #include <atomic>                              << 
 45 #include <cstdlib>  // for std::abs(int)       << 
 46 #include <iostream>                            << 
 47 #include <string.h> // for strcmp                  43 #include <string.h> // for strcmp
 48 #include <string>                              <<  44 #include <cstdlib>  // for std::abs(int)
 49 #include <vector>                              << 
 50                                                    45 
 51 namespace CLHEP {                                  46 namespace CLHEP {
 52                                                    47 
 53 namespace {                                    << 
 54   // Number of instances with automatic seed s << 
 55   CLHEP_ATOMIC_INT_TYPE numberOfEngines(0);    << 
 56                                                << 
 57   // Maximum index into the seed table         << 
 58   const int maxIndex = 215;                    << 
 59 }                                              << 
 60                                                    48 
 61 static const int MarkerLen = 64; // Enough roo     49 static const int MarkerLen = 64; // Enough room to hold a begin or end marker. 
 62                                                    50 
 63 std::string RanluxEngine::name() const {return     51 std::string RanluxEngine::name() const {return "RanluxEngine";}
 64                                                    52 
                                                   >>  53 // Number of instances with automatic seed selection
                                                   >>  54 int RanluxEngine::numEngines = 0;
                                                   >>  55 
                                                   >>  56 // Maximum index into the seed table
                                                   >>  57 int RanluxEngine::maxIndex = 215;
                                                   >>  58 
 65 RanluxEngine::RanluxEngine(long seed, int lux)     59 RanluxEngine::RanluxEngine(long seed, int lux)
 66 : HepRandomEngine()                                60 : HepRandomEngine()
 67 {                                                  61 {
 68    long seedlist[2]={0,0};                         62    long seedlist[2]={0,0};
 69                                                    63 
 70    luxury = lux;                                   64    luxury = lux;
 71    setSeed(seed, luxury);                          65    setSeed(seed, luxury);
 72                                                    66    
 73    // setSeeds() wants a zero terminated array     67    // setSeeds() wants a zero terminated array!
 74    seedlist[0]=theSeed;                            68    seedlist[0]=theSeed;
 75    seedlist[1]=0;                                  69    seedlist[1]=0;
 76    setSeeds(seedlist, luxury);                     70    setSeeds(seedlist, luxury);
 77 }                                                  71 }
 78                                                    72 
 79 RanluxEngine::RanluxEngine()                       73 RanluxEngine::RanluxEngine()
 80 : HepRandomEngine()                                74 : HepRandomEngine()
 81 {                                                  75 {
 82    long seed;                                      76    long seed;
 83    long seedlist[2]={0,0};                         77    long seedlist[2]={0,0};
 84                                                    78 
 85    luxury = 3;                                     79    luxury = 3;
 86    int numEngines = numberOfEngines++;         << 
 87    int cycle = std::abs(int(numEngines/maxInde     80    int cycle = std::abs(int(numEngines/maxIndex));
 88    int curIndex = std::abs(int(numEngines%maxI     81    int curIndex = std::abs(int(numEngines%maxIndex));
 89                                                <<  82    numEngines +=1;
 90    long mask = ((cycle & 0x007fffff) << 8);        83    long mask = ((cycle & 0x007fffff) << 8);
 91    HepRandom::getTheTableSeeds( seedlist, curI     84    HepRandom::getTheTableSeeds( seedlist, curIndex );
 92    seed = seedlist[0]^mask;                        85    seed = seedlist[0]^mask;
 93    setSeed(seed, luxury);                          86    setSeed(seed, luxury);
 94                                                    87    
 95    // setSeeds() wants a zero terminated array     88    // setSeeds() wants a zero terminated array!
 96    seedlist[0]=theSeed;                            89    seedlist[0]=theSeed;
 97    seedlist[1]=0;                                  90    seedlist[1]=0;
 98    setSeeds(seedlist, luxury);                     91    setSeeds(seedlist, luxury);
 99 }                                                  92 }
100                                                    93 
101 RanluxEngine::RanluxEngine(int rowIndex, int c     94 RanluxEngine::RanluxEngine(int rowIndex, int colIndex, int lux)
102 : HepRandomEngine()                                95 : HepRandomEngine()
103 {                                                  96 {
104    long seed;                                      97    long seed;
105    long seedlist[2]={0,0};                         98    long seedlist[2]={0,0};
106                                                    99 
107    luxury = lux;                                  100    luxury = lux;
108    int cycle = std::abs(int(rowIndex/maxIndex)    101    int cycle = std::abs(int(rowIndex/maxIndex));
109    int row = std::abs(int(rowIndex%maxIndex));    102    int row = std::abs(int(rowIndex%maxIndex));
110    int col = std::abs(int(colIndex%2));           103    int col = std::abs(int(colIndex%2));
111    long mask = (( cycle & 0x000007ff ) << 20 )    104    long mask = (( cycle & 0x000007ff ) << 20 );
112    HepRandom::getTheTableSeeds( seedlist, row     105    HepRandom::getTheTableSeeds( seedlist, row );
113    seed = ( seedlist[col] )^mask;                 106    seed = ( seedlist[col] )^mask;
114    setSeed(seed, luxury);                         107    setSeed(seed, luxury);
115                                                   108    
116    // setSeeds() wants a zero terminated array    109    // setSeeds() wants a zero terminated array!
117    seedlist[0]=theSeed;                           110    seedlist[0]=theSeed;
118    seedlist[1]=0;                                 111    seedlist[1]=0;
119    setSeeds(seedlist, luxury);                    112    setSeeds(seedlist, luxury);
120 }                                                 113 }
121                                                   114 
122 RanluxEngine::RanluxEngine( std::istream& is )    115 RanluxEngine::RanluxEngine( std::istream& is )
123 : HepRandomEngine()                               116 : HepRandomEngine()
124 {                                                 117 {
125   is >> *this;                                    118   is >> *this;
126 }                                                 119 }
127                                                   120 
128 RanluxEngine::~RanluxEngine() {}                  121 RanluxEngine::~RanluxEngine() {}
129                                                   122 
130 void RanluxEngine::setSeed(long seed, int lux)    123 void RanluxEngine::setSeed(long seed, int lux) {
131                                                   124 
132 // The initialisation is carried out using a M    125 // The initialisation is carried out using a Multiplicative
133 // Congruential generator using formula consta    126 // Congruential generator using formula constants of L'Ecuyer 
134 // as described in "A review of pseudorandom n    127 // as described in "A review of pseudorandom number generators"
135 // (Fred James) published in Computer Physics     128 // (Fred James) published in Computer Physics Communications 60 (1990)
136 // pages 329-344                                  129 // pages 329-344
137                                                   130 
138   const int ecuyer_a = 53668;                     131   const int ecuyer_a = 53668;
139   const int ecuyer_b = 40014;                     132   const int ecuyer_b = 40014;
140   const int ecuyer_c = 12211;                     133   const int ecuyer_c = 12211;
141   const int ecuyer_d = 2147483563;                134   const int ecuyer_d = 2147483563;
142                                                   135 
143   const int lux_levels[5] = {0,24,73,199,365};    136   const int lux_levels[5] = {0,24,73,199,365};  
144                                                   137 
145   long int_seed_table[24];                        138   long int_seed_table[24];
146   long next_seed = seed;                          139   long next_seed = seed;
147   long k_multiple;                                140   long k_multiple;
148   int i;                                          141   int i;
149                                                   142   
150 // number of additional random numbers that ne    143 // number of additional random numbers that need to be 'thrown away'
151 // every 24 numbers is set using luxury level     144 // every 24 numbers is set using luxury level variable.
152                                                   145 
153   theSeed = seed;                                 146   theSeed = seed;
154   if( (lux > 4)||(lux < 0) ){                     147   if( (lux > 4)||(lux < 0) ){
155      if(lux >= 24){                               148      if(lux >= 24){
156         nskip = lux - 24;                         149         nskip = lux - 24;
157      }else{                                       150      }else{
158         nskip = lux_levels[3]; // corresponds     151         nskip = lux_levels[3]; // corresponds to default luxury level
159      }                                            152      }
160   }else{                                          153   }else{
161      luxury = lux;                                154      luxury = lux;
162      nskip = lux_levels[luxury];                  155      nskip = lux_levels[luxury];
163   }                                               156   }
164                                                   157 
165                                                   158    
166   for(i = 0;i != 24;i++){                         159   for(i = 0;i != 24;i++){
167      k_multiple = next_seed / ecuyer_a;           160      k_multiple = next_seed / ecuyer_a;
168      next_seed = ecuyer_b * (next_seed - k_mul    161      next_seed = ecuyer_b * (next_seed - k_multiple * ecuyer_a) 
169      - k_multiple * ecuyer_c ;                    162      - k_multiple * ecuyer_c ;
170      if(next_seed < 0)next_seed += ecuyer_d;      163      if(next_seed < 0)next_seed += ecuyer_d;
171      int_seed_table[i] = next_seed % int_modul    164      int_seed_table[i] = next_seed % int_modulus;
172   }                                               165   }     
173                                                   166 
174   for(i = 0;i != 24;i++)                          167   for(i = 0;i != 24;i++)
175     float_seed_table[i] = int_seed_table[i] *     168     float_seed_table[i] = int_seed_table[i] * mantissa_bit_24();
176                                                   169 
177   i_lag = 23;                                     170   i_lag = 23;
178   j_lag = 9;                                      171   j_lag = 9;
179   carry = 0. ;                                    172   carry = 0. ;
180                                                   173 
181   if( float_seed_table[23] == 0. ) carry = man    174   if( float_seed_table[23] == 0. ) carry = mantissa_bit_24();
182                                                   175   
183   count24 = 0;                                    176   count24 = 0;
184 }                                                 177 }
185                                                   178 
186 void RanluxEngine::setSeeds(const long *seeds,    179 void RanluxEngine::setSeeds(const long *seeds, int lux) {
187                                                   180 
188    const int ecuyer_a = 53668;                    181    const int ecuyer_a = 53668;
189    const int ecuyer_b = 40014;                    182    const int ecuyer_b = 40014;
190    const int ecuyer_c = 12211;                    183    const int ecuyer_c = 12211;
191    const int ecuyer_d = 2147483563;               184    const int ecuyer_d = 2147483563;
192                                                   185 
193    const int lux_levels[5] = {0,24,73,199,365}    186    const int lux_levels[5] = {0,24,73,199,365};
194    int i;                                         187    int i;
195    long int_seed_table[24];                       188    long int_seed_table[24];
196    long k_multiple,next_seed;                     189    long k_multiple,next_seed;
197    const long *seedptr;                           190    const long *seedptr; 
198                                                   191 
199    theSeeds = seeds;                              192    theSeeds = seeds;
200    seedptr  = seeds;                              193    seedptr  = seeds;
201                                                   194  
202    if(seeds == 0){                                195    if(seeds == 0){
203       setSeed(theSeed,lux);                       196       setSeed(theSeed,lux);
204       theSeeds = &theSeed;                        197       theSeeds = &theSeed;
205       return;                                     198       return;
206    }                                              199    }
207                                                   200 
208    theSeed = *seeds;                              201    theSeed = *seeds;
209                                                   202 
210 // number of additional random numbers that ne    203 // number of additional random numbers that need to be 'thrown away'
211 // every 24 numbers is set using luxury level     204 // every 24 numbers is set using luxury level variable.
212                                                   205 
213   if( (lux > 4)||(lux < 0) ){                     206   if( (lux > 4)||(lux < 0) ){
214      if(lux >= 24){                               207      if(lux >= 24){
215         nskip = lux - 24;                         208         nskip = lux - 24;
216      }else{                                       209      }else{
217         nskip = lux_levels[3]; // corresponds     210         nskip = lux_levels[3]; // corresponds to default luxury level
218      }                                            211      }
219   }else{                                          212   }else{
220      luxury = lux;                                213      luxury = lux;
221      nskip = lux_levels[luxury];                  214      nskip = lux_levels[luxury];
222   }                                               215   }
223                                                   216       
224   for( i = 0;(i != 24)&&(*seedptr != 0);i++){     217   for( i = 0;(i != 24)&&(*seedptr != 0);i++){
225       int_seed_table[i] = *seedptr % int_modul    218       int_seed_table[i] = *seedptr % int_modulus;
226       seedptr++;                                  219       seedptr++;
227   }                                               220   }          
228                                                   221 
229   if(i != 24){                                    222   if(i != 24){
230      next_seed = int_seed_table[i-1];             223      next_seed = int_seed_table[i-1];
231      for(;i != 24;i++){                           224      for(;i != 24;i++){
232         k_multiple = next_seed / ecuyer_a;        225         k_multiple = next_seed / ecuyer_a;
233         next_seed = ecuyer_b * (next_seed - k_    226         next_seed = ecuyer_b * (next_seed - k_multiple * ecuyer_a) 
234         - k_multiple * ecuyer_c ;                 227         - k_multiple * ecuyer_c ;
235         if(next_seed < 0)next_seed += ecuyer_d    228         if(next_seed < 0)next_seed += ecuyer_d;
236         int_seed_table[i] = next_seed % int_mo    229         int_seed_table[i] = next_seed % int_modulus;
237      }                                            230      }          
238   }                                               231   }
239                                                   232 
240   for(i = 0;i != 24;i++)                          233   for(i = 0;i != 24;i++)
241     float_seed_table[i] = int_seed_table[i] *     234     float_seed_table[i] = int_seed_table[i] * mantissa_bit_24();
242                                                   235 
243   i_lag = 23;                                     236   i_lag = 23;
244   j_lag = 9;                                      237   j_lag = 9;
245   carry = 0. ;                                    238   carry = 0. ;
246                                                   239 
247   if( float_seed_table[23] == 0. ) carry = man    240   if( float_seed_table[23] == 0. ) carry = mantissa_bit_24();
248                                                   241   
249   count24 = 0;                                    242   count24 = 0;
250 }                                                 243 }
251                                                   244 
252 void RanluxEngine::saveStatus( const char file    245 void RanluxEngine::saveStatus( const char filename[] ) const
253 {                                                 246 {
254    std::ofstream outFile( filename, std::ios::    247    std::ofstream outFile( filename, std::ios::out ) ;
255   if (!outFile.bad()) {                           248   if (!outFile.bad()) {
256     outFile << "Uvec\n";                          249     outFile << "Uvec\n";
257     std::vector<unsigned long> v = put();         250     std::vector<unsigned long> v = put();
258     for (unsigned int i=0; i<v.size(); ++i) {     251     for (unsigned int i=0; i<v.size(); ++i) {
259       outFile << v[i] << "\n";                    252       outFile << v[i] << "\n";
260     }                                             253     }
261   }                                               254   }
262 }                                                 255 }
263                                                   256 
264 void RanluxEngine::restoreStatus( const char f    257 void RanluxEngine::restoreStatus( const char filename[] )
265 {                                                 258 {
266    std::ifstream inFile( filename, std::ios::i    259    std::ifstream inFile( filename, std::ios::in);
267    if (!checkFile ( inFile, filename, engineNa    260    if (!checkFile ( inFile, filename, engineName(), "restoreStatus" )) {
268      std::cerr << "  -- Engine state remains u    261      std::cerr << "  -- Engine state remains unchanged\n";
269      return;                                      262      return;
270    }                                              263    }
271   if ( possibleKeywordInput ( inFile, "Uvec",     264   if ( possibleKeywordInput ( inFile, "Uvec", theSeed ) ) {
272     std::vector<unsigned long> v;                 265     std::vector<unsigned long> v;
273     unsigned long xin;                            266     unsigned long xin;
274     for (unsigned int ivec=0; ivec < VECTOR_ST    267     for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) {
275       inFile >> xin;                              268       inFile >> xin;
276       if (!inFile) {                              269       if (!inFile) {
277         inFile.clear(std::ios::badbit | inFile    270         inFile.clear(std::ios::badbit | inFile.rdstate());
278         std::cerr << "\nRanluxEngine state (ve    271         std::cerr << "\nRanluxEngine state (vector) description improper."
279          << "\nrestoreStatus has failed."         272          << "\nrestoreStatus has failed."
280          << "\nInput stream is probably mispos    273          << "\nInput stream is probably mispositioned now." << std::endl;
281         return;                                   274         return;
282       }                                           275       }
283       v.push_back(xin);                           276       v.push_back(xin);
284     }                                             277     }
285     getState(v);                                  278     getState(v);
286     return;                                       279     return;
287   }                                               280   }
288                                                   281 
289    if (!inFile.bad() && !inFile.eof()) {          282    if (!inFile.bad() && !inFile.eof()) {
290 //     inFile >> theSeed;  removed -- encompas    283 //     inFile >> theSeed;  removed -- encompased by possibleKeywordInput
291      for (int i=0; i<24; ++i)                     284      for (int i=0; i<24; ++i)
292        inFile >> float_seed_table[i];             285        inFile >> float_seed_table[i];
293      inFile >> i_lag; inFile >> j_lag;            286      inFile >> i_lag; inFile >> j_lag;
294      inFile >> carry; inFile >> count24;          287      inFile >> carry; inFile >> count24;
295      inFile >> luxury; inFile >> nskip;           288      inFile >> luxury; inFile >> nskip;
296    }                                              289    }
297 }                                                 290 }
298                                                   291 
299 void RanluxEngine::showStatus() const             292 void RanluxEngine::showStatus() const
300 {                                                 293 {
301    std::cout << std::endl;                        294    std::cout << std::endl;
302    std::cout << "--------- Ranlux engine statu    295    std::cout << "--------- Ranlux engine status ---------" << std::endl;
303    std::cout << " Initial seed = " << theSeed     296    std::cout << " Initial seed = " << theSeed << std::endl;
304    std::cout << " float_seed_table[] = ";         297    std::cout << " float_seed_table[] = ";
305    for (int i=0; i<24; ++i)                       298    for (int i=0; i<24; ++i)
306      std::cout << float_seed_table[i] << " ";     299      std::cout << float_seed_table[i] << " ";
307    std::cout << std::endl;                        300    std::cout << std::endl;
308    std::cout << " i_lag = " << i_lag << ", j_l    301    std::cout << " i_lag = " << i_lag << ", j_lag = " << j_lag << std::endl;
309    std::cout << " carry = " << carry << ", cou    302    std::cout << " carry = " << carry << ", count24 = " << count24 << std::endl;
310    std::cout << " luxury = " << luxury << " ns    303    std::cout << " luxury = " << luxury << " nskip = " << nskip << std::endl;
311    std::cout << "-----------------------------    304    std::cout << "----------------------------------------" << std::endl;
312 }                                                 305 }
313                                                   306 
314 double RanluxEngine::flat() {                     307 double RanluxEngine::flat() {
315                                                   308 
316   float next_random;                              309   float next_random;
317   float uni;                                      310   float uni;
318   int i;                                          311   int i;
319                                                   312 
320   uni = float_seed_table[j_lag] - float_seed_t    313   uni = float_seed_table[j_lag] - float_seed_table[i_lag] - carry;
321   if(uni < 0. ){                                  314   if(uni < 0. ){
322      uni += 1.0;                                  315      uni += 1.0;
323      carry = mantissa_bit_24();                   316      carry = mantissa_bit_24();
324   }else{                                          317   }else{
325      carry = 0.;                                  318      carry = 0.;
326   }                                               319   }
327                                                   320 
328   float_seed_table[i_lag] = uni;                  321   float_seed_table[i_lag] = uni;
329   i_lag --;                                       322   i_lag --;
330   j_lag --;                                       323   j_lag --;
331   if(i_lag < 0) i_lag = 23;                       324   if(i_lag < 0) i_lag = 23;
332   if(j_lag < 0) j_lag = 23;                       325   if(j_lag < 0) j_lag = 23;
333                                                   326 
334   if( uni < mantissa_bit_12() ){                  327   if( uni < mantissa_bit_12() ){
335      uni += mantissa_bit_24() * float_seed_tab    328      uni += mantissa_bit_24() * float_seed_table[j_lag];
336      if( uni == 0) uni = mantissa_bit_24() * m    329      if( uni == 0) uni = mantissa_bit_24() * mantissa_bit_24();
337   }                                               330   }
338   next_random = uni;                              331   next_random = uni;
339   count24 ++;                                     332   count24 ++;
340                                                   333 
341 // every 24th number generation, several rando    334 // every 24th number generation, several random numbers are generated
342 // and wasted depending upon the luxury level.    335 // and wasted depending upon the luxury level.
343                                                   336 
344   if(count24 == 24 ){                             337   if(count24 == 24 ){
345      count24 = 0;                                 338      count24 = 0;
346      for( i = 0; i != nskip ; i++){               339      for( i = 0; i != nskip ; i++){
347          uni = float_seed_table[j_lag] - float    340          uni = float_seed_table[j_lag] - float_seed_table[i_lag] - carry;
348          if(uni < 0. ){                           341          if(uni < 0. ){
349             uni += 1.0;                           342             uni += 1.0;
350             carry = mantissa_bit_24();            343             carry = mantissa_bit_24();
351          }else{                                   344          }else{
352             carry = 0.;                           345             carry = 0.;
353          }                                        346          }
354          float_seed_table[i_lag] = uni;           347          float_seed_table[i_lag] = uni;
355    i_lag --;                                      348    i_lag --;
356          j_lag --;                                349          j_lag --;
357          if(i_lag < 0)i_lag = 23;                 350          if(i_lag < 0)i_lag = 23;
358          if(j_lag < 0) j_lag = 23;                351          if(j_lag < 0) j_lag = 23;
359       }                                           352       }
360   }                                               353   } 
361   return (double) next_random;                    354   return (double) next_random;
362 }                                                 355 }
363                                                   356 
364 void RanluxEngine::flatArray(const int size, d    357 void RanluxEngine::flatArray(const int size, double* vect)
365 {                                                 358 {
366   float next_random;                              359   float next_random;
367   float uni;                                      360   float uni;
368   int i;                                          361   int i;
369   int index;                                      362   int index;
370                                                   363 
371   for (index=0; index<size; ++index) {            364   for (index=0; index<size; ++index) {
372     uni = float_seed_table[j_lag] - float_seed    365     uni = float_seed_table[j_lag] - float_seed_table[i_lag] - carry;
373     if(uni < 0. ){                                366     if(uni < 0. ){
374        uni += 1.0;                                367        uni += 1.0;
375        carry = mantissa_bit_24();                 368        carry = mantissa_bit_24();
376     }else{                                        369     }else{
377        carry = 0.;                                370        carry = 0.;
378     }                                             371     }
379                                                   372 
380     float_seed_table[i_lag] = uni;                373     float_seed_table[i_lag] = uni;
381     i_lag --;                                     374     i_lag --;
382     j_lag --;                                     375     j_lag --;
383     if(i_lag < 0) i_lag = 23;                     376     if(i_lag < 0) i_lag = 23;
384     if(j_lag < 0) j_lag = 23;                     377     if(j_lag < 0) j_lag = 23;
385                                                   378 
386     if( uni < mantissa_bit_12() ){                379     if( uni < mantissa_bit_12() ){
387        uni += mantissa_bit_24() * float_seed_t    380        uni += mantissa_bit_24() * float_seed_table[j_lag];
388        if( uni == 0) uni = mantissa_bit_24() *    381        if( uni == 0) uni = mantissa_bit_24() * mantissa_bit_24();
389     }                                             382     }
390     next_random = uni;                            383     next_random = uni;
391     vect[index] = (double)next_random;            384     vect[index] = (double)next_random;
392     count24 ++;                                   385     count24 ++;
393                                                   386 
394 // every 24th number generation, several rando    387 // every 24th number generation, several random numbers are generated
395 // and wasted depending upon the luxury level.    388 // and wasted depending upon the luxury level.
396                                                   389 
397     if(count24 == 24 ){                           390     if(count24 == 24 ){
398        count24 = 0;                               391        count24 = 0;
399        for( i = 0; i != nskip ; i++){             392        for( i = 0; i != nskip ; i++){
400            uni = float_seed_table[j_lag] - flo    393            uni = float_seed_table[j_lag] - float_seed_table[i_lag] - carry;
401            if(uni < 0. ){                         394            if(uni < 0. ){
402               uni += 1.0;                         395               uni += 1.0;
403               carry = mantissa_bit_24();          396               carry = mantissa_bit_24();
404            }else{                                 397            }else{
405               carry = 0.;                         398               carry = 0.;
406            }                                      399            }
407            float_seed_table[i_lag] = uni;         400            float_seed_table[i_lag] = uni;
408            i_lag --;                              401            i_lag --;
409            j_lag --;                              402            j_lag --;
410            if(i_lag < 0)i_lag = 23;               403            if(i_lag < 0)i_lag = 23;
411            if(j_lag < 0) j_lag = 23;              404            if(j_lag < 0) j_lag = 23;
412         }                                         405         }
413     }                                             406     }
414   }                                               407   }
415 }                                                 408 } 
416                                                   409 
417 RanluxEngine::operator double() {              << 
418   return flat();                               << 
419 }                                              << 
420                                                << 
421 RanluxEngine::operator float() {               << 
422   return float( flat() );                      << 
423 }                                              << 
424                                                << 
425 RanluxEngine::operator unsigned int() {           410 RanluxEngine::operator unsigned int() {
426    return ((unsigned int)(flat() * exponent_bi    411    return ((unsigned int)(flat() * exponent_bit_32()) & 0xffffffff) |
427          (((unsigned int)(float_seed_table[i_l    412          (((unsigned int)(float_seed_table[i_lag]*exponent_bit_32())>>16) & 0xff);
428    // needed because Ranlux doesn't fill all b    413    // needed because Ranlux doesn't fill all bits of the double
429    // which therefore doesn't fill all bits of    414    // which therefore doesn't fill all bits of the integer.
430 }                                                 415 }
431                                                   416 
432 std::ostream & RanluxEngine::put ( std::ostrea    417 std::ostream & RanluxEngine::put ( std::ostream& os ) const
433 {                                                 418 {
434    char beginMarker[] = "RanluxEngine-begin";     419    char beginMarker[] = "RanluxEngine-begin";
435   os << beginMarker << "\nUvec\n";                420   os << beginMarker << "\nUvec\n";
436   std::vector<unsigned long> v = put();           421   std::vector<unsigned long> v = put();
437   for (unsigned int i=0; i<v.size(); ++i) {       422   for (unsigned int i=0; i<v.size(); ++i) {
438      os <<  v[i] <<  "\n";                        423      os <<  v[i] <<  "\n";
439   }                                               424   }
440   return os;                                      425   return os;  
441 }                                                 426 }
442                                                   427 
443 std::vector<unsigned long> RanluxEngine::put (    428 std::vector<unsigned long> RanluxEngine::put () const {
444   std::vector<unsigned long> v;                   429   std::vector<unsigned long> v;
445   v.push_back (engineIDulong<RanluxEngine>());    430   v.push_back (engineIDulong<RanluxEngine>());
446   for (int i=0; i<24; ++i) {                      431   for (int i=0; i<24; ++i) {
447     v.push_back                                   432     v.push_back
448       (static_cast<unsigned long>(float_seed_t    433       (static_cast<unsigned long>(float_seed_table[i]/mantissa_bit_24()));
449   }                                               434   }
450   v.push_back(static_cast<unsigned long>(i_lag    435   v.push_back(static_cast<unsigned long>(i_lag));
451   v.push_back(static_cast<unsigned long>(j_lag    436   v.push_back(static_cast<unsigned long>(j_lag));
452   v.push_back(static_cast<unsigned long>(carry    437   v.push_back(static_cast<unsigned long>(carry/mantissa_bit_24()));
453   v.push_back(static_cast<unsigned long>(count    438   v.push_back(static_cast<unsigned long>(count24));
454   v.push_back(static_cast<unsigned long>(luxur    439   v.push_back(static_cast<unsigned long>(luxury));
455   v.push_back(static_cast<unsigned long>(nskip    440   v.push_back(static_cast<unsigned long>(nskip));
456   return v;                                       441   return v;
457 }                                                 442 }
458                                                   443 
459 std::istream & RanluxEngine::get ( std::istrea    444 std::istream & RanluxEngine::get ( std::istream& is )
460 {                                                 445 {
461   char beginMarker [MarkerLen];                   446   char beginMarker [MarkerLen];
462   is >> std::ws;                                  447   is >> std::ws;
463   is.width(MarkerLen);  // causes the next rea    448   is.width(MarkerLen);  // causes the next read to the char* to be <=
464       // that many bytes, INCLUDING A TERMINAT    449       // that many bytes, INCLUDING A TERMINATION \0 
465       // (Stroustrup, section 21.3.2)             450       // (Stroustrup, section 21.3.2)
466   is >> beginMarker;                              451   is >> beginMarker;
467   if (strcmp(beginMarker,"RanluxEngine-begin")    452   if (strcmp(beginMarker,"RanluxEngine-begin")) {
468      is.clear(std::ios::badbit | is.rdstate())    453      is.clear(std::ios::badbit | is.rdstate());
469      std::cerr << "\nInput stream mispositione    454      std::cerr << "\nInput stream mispositioned or"
470          << "\nRanluxEngine state description     455          << "\nRanluxEngine state description missing or"
471          << "\nwrong engine type found." << st    456          << "\nwrong engine type found." << std::endl;
472      return is;                                   457      return is;
473   }                                               458   }
474   return getState(is);                            459   return getState(is);
475 }                                                 460 }
476                                                   461 
477 std::string RanluxEngine::beginTag ( )  {         462 std::string RanluxEngine::beginTag ( )  { 
478   return "RanluxEngine-begin";                    463   return "RanluxEngine-begin"; 
479 }                                                 464 }
480                                                   465 
481 std::istream & RanluxEngine::getState ( std::i    466 std::istream & RanluxEngine::getState ( std::istream& is )
482 {                                                 467 {
483   if ( possibleKeywordInput ( is, "Uvec", theS    468   if ( possibleKeywordInput ( is, "Uvec", theSeed ) ) {
484     std::vector<unsigned long> v;                 469     std::vector<unsigned long> v;
485     unsigned long uu;                             470     unsigned long uu;
486     for (unsigned int ivec=0; ivec < VECTOR_ST    471     for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) {
487       is >> uu;                                   472       is >> uu;
488       if (!is) {                                  473       if (!is) {
489         is.clear(std::ios::badbit | is.rdstate    474         is.clear(std::ios::badbit | is.rdstate());
490         std::cerr << "\nRanluxEngine state (ve    475         std::cerr << "\nRanluxEngine state (vector) description improper."
491     << "\ngetState() has failed."                 476     << "\ngetState() has failed."
492          << "\nInput stream is probably mispos    477          << "\nInput stream is probably mispositioned now." << std::endl;
493         return is;                                478         return is;
494       }                                           479       }
495       v.push_back(uu);                            480       v.push_back(uu);
496     }                                             481     }
497     getState(v);                                  482     getState(v);
498     return (is);                                  483     return (is);
499   }                                               484   }
500                                                   485 
501 //  is >> theSeed;  Removed, encompassed by po    486 //  is >> theSeed;  Removed, encompassed by possibleKeywordInput()
502                                                   487 
503   char endMarker   [MarkerLen];                   488   char endMarker   [MarkerLen];
504   for (int i=0; i<24; ++i) {                      489   for (int i=0; i<24; ++i) {
505      is >> float_seed_table[i];                   490      is >> float_seed_table[i];
506   }                                               491   }
507   is >> i_lag; is >>  j_lag;                      492   is >> i_lag; is >>  j_lag;
508   is >> carry; is >> count24;                     493   is >> carry; is >> count24;
509   is >> luxury; is >> nskip;                      494   is >> luxury; is >> nskip;
510   is >> std::ws;                                  495   is >> std::ws;
511   is.width(MarkerLen);                            496   is.width(MarkerLen);  
512   is >> endMarker;                                497   is >> endMarker;
513   if (strcmp(endMarker,"RanluxEngine-end")) {     498   if (strcmp(endMarker,"RanluxEngine-end")) {
514      is.clear(std::ios::badbit | is.rdstate())    499      is.clear(std::ios::badbit | is.rdstate());
515      std::cerr << "\nRanluxEngine state descri    500      std::cerr << "\nRanluxEngine state description incomplete."
516          << "\nInput stream is probably mispos    501          << "\nInput stream is probably mispositioned now." << std::endl;
517      return is;                                   502      return is;
518   }                                               503   }
519   return is;                                      504   return is;
520 }                                                 505 }
521                                                   506 
522 bool RanluxEngine::get (const std::vector<unsi    507 bool RanluxEngine::get (const std::vector<unsigned long> & v) {
523   if ((v[0] & 0xffffffffUL) != engineIDulong<R    508   if ((v[0] & 0xffffffffUL) != engineIDulong<RanluxEngine>()) {
524     std::cerr <<                                  509     std::cerr << 
525       "\nRanluxEngine get:state vector has wro    510       "\nRanluxEngine get:state vector has wrong ID word - state unchanged\n";
526     return false;                                 511     return false;
527   }                                               512   }
528   return getState(v);                             513   return getState(v);
529 }                                                 514 }
530                                                   515 
531 bool RanluxEngine::getState (const std::vector    516 bool RanluxEngine::getState (const std::vector<unsigned long> & v) {
532   if (v.size() != VECTOR_STATE_SIZE ) {           517   if (v.size() != VECTOR_STATE_SIZE ) {
533     std::cerr <<                                  518     std::cerr << 
534       "\nRanluxEngine get:state vector has wro    519       "\nRanluxEngine get:state vector has wrong length - state unchanged\n";
535     return false;                                 520     return false;
536   }                                               521   }
537   for (int i=0; i<24; ++i) {                      522   for (int i=0; i<24; ++i) {
538     float_seed_table[i] = v[i+1]*mantissa_bit_    523     float_seed_table[i] = v[i+1]*mantissa_bit_24();
539   }                                               524   }
540   i_lag    = (int)v[25];                       << 525   i_lag    = v[25];
541   j_lag    = (int)v[26];                       << 526   j_lag    = v[26];
542   carry    = v[27]*mantissa_bit_24();             527   carry    = v[27]*mantissa_bit_24();
543   count24  = (int)v[28];                       << 528   count24  = v[28];
544   luxury   = (int)v[29];                       << 529   luxury   = v[29];
545   nskip    = (int)v[30];                       << 530   nskip    = v[30];
546   return true;                                    531   return true;
547 }                                                 532 }
548                                                   533 
549 }  // namespace CLHEP                             534 }  // namespace CLHEP
550                                                   535