Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/externals/clhep/src/JamesRandom.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 ]

  1 // -*- C++ -*-
  2 //
  3 // -----------------------------------------------------------------------
  4 //                             HEP Random
  5 //                       --- HepJamesRandom ---
  6 //                      class implementation file
  7 // -----------------------------------------------------------------------
  8 // This file is part of Geant4 (simulation toolkit for HEP).
  9 //
 10 // This algorithm implements the original universal random number generator
 11 // as proposed by Marsaglia & Zaman in report FSU-SCRI-87-50 and coded
 12 // in FORTRAN77 by Fred James as the RANMAR generator, part of the MATHLIB
 13 // HEP library.
 14 
 15 // =======================================================================
 16 // Gabriele Cosmo - Created: 5th September 1995
 17 //                - Fixed a bug in setSeed(): 26th February 1996
 18 //                - Minor corrections: 31st October 1996
 19 //                - Added methods for engine status: 19th November 1996
 20 //                - Fixed bug in setSeeds(): 15th September 1997
 21 // J.Marraffino   - Added stream operators and related constructor.
 22 //                  Added automatic seed selection from seed table and
 23 //                  engine counter: 16th Feb 1998
 24 // Ken Smith      - Added conversion operators:  6th Aug 1998
 25 // J. Marraffino  - Remove dependence on hepString class  13 May 1999
 26 // V. Innocente   - changed pointers to indices     3 may 2000
 27 // M. Fischler    - In restore, checkFile for file not found    03 Dec 2004
 28 // M. Fischler    - Methods for distrib. instacne save/restore  12/8/04    
 29 // M. Fischler    - split get() into tag validation and 
 30 //                  getState() for anonymous restores           12/27/04    
 31 // M. Fischler    - Enforcement that seeds be non-negative
 32 //        (lest the sequence be non-random)          2/14/05    
 33 // M. Fischler    - put/get for vectors of ulongs   3/14/05
 34 // M. Fischler    - State-saving using only ints, for portability 4/12/05
 35 //        
 36 // =======================================================================
 37 
 38 #include "CLHEP/Random/Random.h"
 39 #include "CLHEP/Random/JamesRandom.h"
 40 #include "CLHEP/Random/engineIDulong.h"
 41 #include "CLHEP/Random/DoubConv.h"
 42 #include "CLHEP/Utility/atomic_int.h"
 43 
 44 #include <atomic>
 45 #include <cmath>
 46 #include <cstdlib>
 47 #include <iostream>
 48 #include <string.h> // for strcmp
 49 #include <string>
 50 #include <vector>
 51 
 52 namespace CLHEP {
 53 
 54 namespace {
 55   // Number of instances with automatic seed selection
 56   CLHEP_ATOMIC_INT_TYPE numberOfEngines(0);
 57 
 58   // Maximum index into the seed table
 59   const int maxIndex = 215;
 60 }
 61 
 62 static const int MarkerLen = 64; // Enough room to hold a begin or end marker. 
 63 
 64 std::string HepJamesRandom::name() const {return "HepJamesRandom";}
 65 
 66 HepJamesRandom::HepJamesRandom(long seed)
 67 : HepRandomEngine()
 68 {
 69   setSeed(seed,0);
 70   setSeeds(&theSeed,0);
 71 }
 72 
 73 HepJamesRandom::HepJamesRandom()      // 15 Feb. 1998  JMM
 74 : HepRandomEngine()
 75 {
 76   long seeds[2];
 77   long seed;
 78 
 79   int numEngines = numberOfEngines++;
 80   int cycle = std::abs(int(numEngines/maxIndex));
 81   int curIndex = std::abs(int(numEngines%maxIndex));
 82 
 83   long mask = ((cycle & 0x007fffff) << 8);
 84   HepRandom::getTheTableSeeds( seeds, curIndex );
 85   seed = seeds[0]^mask;
 86   setSeed(seed,0);
 87   setSeeds(&theSeed,0);
 88 }
 89 
 90 HepJamesRandom::HepJamesRandom(int rowIndex, int colIndex) // 15 Feb. 1998  JMM
 91 : HepRandomEngine()
 92 {
 93   long seed;
 94    long seeds[2];
 95 
 96   int cycle = std::abs(int(rowIndex/maxIndex));
 97   int row = std::abs(int(rowIndex%maxIndex));
 98   int col = std::abs(int(colIndex%2));
 99   long mask = ((cycle & 0x000007ff) << 20);
100   HepRandom::getTheTableSeeds( seeds, row );
101   seed = (seeds[col])^mask;
102   setSeed(seed,0);
103   setSeeds(&theSeed,0);
104 }
105 
106 HepJamesRandom::HepJamesRandom(std::istream& is)
107 : HepRandomEngine()
108 {
109   is >> *this;
110 }
111 
112 HepJamesRandom::~HepJamesRandom() {}
113 
114 void HepJamesRandom::saveStatus( const char filename[] ) const
115 {
116   std::ofstream outFile( filename, std::ios::out ) ;
117 
118   if (!outFile.bad()) {
119     outFile << "Uvec\n";
120     std::vector<unsigned long> v = put();
121     for (unsigned int i=0; i<v.size(); ++i) {
122       outFile << v[i] << "\n";
123     }
124   }
125 }
126 
127 void HepJamesRandom::restoreStatus( const char filename[] )
128 {
129    int ipos, jpos;
130    std::ifstream inFile( filename, std::ios::in);
131    if (!checkFile ( inFile, filename, engineName(), "restoreStatus" )) {
132      std::cerr << "  -- Engine state remains unchanged\n";
133      return;
134    }
135   if ( possibleKeywordInput ( inFile, "Uvec", theSeed ) ) {
136     std::vector<unsigned long> v;
137     unsigned long xin;
138     for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) {
139       inFile >> xin;
140       if (!inFile) {
141         inFile.clear(std::ios::badbit | inFile.rdstate());
142         std::cerr << "\nJamesRandom state (vector) description improper."
143          << "\nrestoreStatus has failed."
144          << "\nInput stream is probably mispositioned now." << std::endl;
145         return;
146       }
147       v.push_back(xin);
148     }
149     getState(v);
150     return;
151   }
152 
153    if (!inFile.bad() && !inFile.eof()) {
154 //     inFile >> theSeed;  removed -- encompased by possibleKeywordInput
155      for (int i=0; i<97; ++i)
156        inFile >> u[i];
157      inFile >> c; inFile >> cd; inFile >> cm;
158      inFile >> jpos;
159      ipos = (64+jpos)%97;
160      i97 = ipos;
161      j97 = jpos;
162    }
163 }
164 
165 void HepJamesRandom::showStatus() const
166 {
167    std::cout << std::endl;
168    std::cout << "----- HepJamesRandom engine status -----" << std::endl;
169    std::cout << " Initial seed = " << theSeed << std::endl;
170    std::cout << " u[] = ";
171    for (int i=0; i<97; ++i)
172      std::cout << u[i] << " ";
173    std::cout << std::endl;
174    std::cout << " c = " << c << ", cd = " << cd << ", cm = " << cm
175        << std::endl;
176    std::cout << " i97 = " << i97 << ", u[i97] = " << u[i97] << std::endl;
177    std::cout << " j97 = " << j97 << ", u[j97] = " << u[j97] << std::endl;
178    std::cout << "----------------------------------------" << std::endl;
179 }
180 
181 void HepJamesRandom::setSeed(long seed, int)
182 {
183   // The input value for "seed" should be within the range [0,900000000]
184   //
185   // Negative seeds result in serious flaws in the randomness;
186   // seeds above 900000000 are OK because of the %177 in the expression for i,
187   // but may have the same effect as other seeds below 900000000.
188 
189   int m, n;
190   float s, t;
191   long mm;
192 
193   if (seed < 0) {
194     std::cout << "Seed for HepJamesRandom must be non-negative\n" 
195       << "Seed value supplied was " << seed  
196   << "\nUsing its absolute value instead\n";
197     seed = -seed;
198   }
199   
200   long ij = seed/30082;
201   long kl = seed - 30082*ij;
202   long i = (ij/177) % 177 + 2;
203   long j = ij % 177 + 2;
204   long k = (kl/169) % 178 + 1;
205   long l = kl % 169;
206 
207   theSeed = seed;
208 
209   for ( n = 1 ; n < 98 ; n++ ) {
210     s = 0.0;
211     t = 0.5;
212     for ( m = 1 ; m < 25 ; m++) {
213       mm = ( ( (i*j) % 179 ) * k ) % 179;
214       i = j;
215       j = k;
216       k = mm;
217       l = ( 53 * l + 1 ) % 169;
218       if ( (l*mm % 64 ) >= 32 )
219         s += t;
220       t *= 0.5;
221     }
222     u[n-1] = s;
223   }
224   c = 362436.0 / 16777216.0;
225   cd = 7654321.0 / 16777216.0;
226   cm = 16777213.0 / 16777216.0;
227 
228   i97 = 96;
229   j97 = 32;
230 
231 }
232 
233 void HepJamesRandom::setSeeds(const long* seeds, int)
234 {
235   setSeed(seeds ? *seeds : 19780503L, 0);
236   theSeeds = seeds;
237 }
238 
239 double HepJamesRandom::flat()
240 {
241    double uni;
242 
243    do {
244       uni = u[i97] - u[j97];
245       if ( uni < 0.0 ) uni++;
246       u[i97] = uni;
247       
248       if (i97 == 0) i97 = 96;
249       else i97--;
250       
251       if (j97 == 0) j97 = 96;
252       else j97--;
253       
254       c -= cd;
255       if (c < 0.0) c += cm;
256       
257       uni -= c;
258       if (uni < 0.0) uni += 1.0;
259    } while ( uni <= 0.0 || uni >= 1.0 );
260    
261    return uni;
262 }
263 
264 void HepJamesRandom::flatArray(const int size, double* vect)
265 {
266 //   double uni;
267    int i;
268 
269    for (i=0; i<size; ++i) {
270      vect[i] = flat();
271    }   
272 }
273 
274 HepJamesRandom::operator double() {
275    return flat();
276 }
277 
278 HepJamesRandom::operator float() {
279    return float( flat() );
280 }
281 
282 HepJamesRandom::operator unsigned int() {
283    return ((unsigned int)(flat() * exponent_bit_32()) & 0xffffffff )  |
284          (((unsigned int)( u[i97] * exponent_bit_32())>>16)  & 0xff);
285 }
286 
287 std::ostream & HepJamesRandom::put ( std::ostream& os ) const {
288   char beginMarker[] = "JamesRandom-begin";
289   os << beginMarker << "\nUvec\n";
290   std::vector<unsigned long> v = put();
291   for (unsigned int i=0; i<v.size(); ++i) {
292      os <<  v[i] <<  "\n";
293   }
294   return os;  
295 }
296 
297 std::vector<unsigned long> HepJamesRandom::put () const {
298   std::vector<unsigned long> v;
299   v.push_back (engineIDulong<HepJamesRandom>());
300   std::vector<unsigned long> t;
301   for (int i=0; i<97; ++i) {
302     t = DoubConv::dto2longs(u[i]);
303     v.push_back(t[0]); v.push_back(t[1]);
304   }
305   t = DoubConv::dto2longs(c);
306   v.push_back(t[0]); v.push_back(t[1]);
307   t = DoubConv::dto2longs(cd);
308   v.push_back(t[0]); v.push_back(t[1]);
309   t = DoubConv::dto2longs(cm);
310   v.push_back(t[0]); v.push_back(t[1]);
311   v.push_back(static_cast<unsigned long>(j97));
312   return v;
313 }
314 
315 
316 std::istream & HepJamesRandom::get  ( std::istream& is) {
317   char beginMarker [MarkerLen];
318   is >> std::ws;
319   is.width(MarkerLen);  // causes the next read to the char* to be <=
320       // that many bytes, INCLUDING A TERMINATION \0 
321       // (Stroustrup, section 21.3.2)
322   is >> beginMarker;
323   if (strcmp(beginMarker,"JamesRandom-begin")) {
324      is.clear(std::ios::badbit | is.rdstate());
325      std::cerr << "\nInput stream mispositioned or"
326          << "\nJamesRandom state description missing or"
327          << "\nwrong engine type found." << std::endl;
328      return is;
329   }
330   return getState(is);
331 }
332 
333 std::string HepJamesRandom::beginTag ( )  { 
334   return "JamesRandom-begin"; 
335 }
336 
337 std::istream & HepJamesRandom::getState  ( std::istream& is) {
338   if ( possibleKeywordInput ( is, "Uvec", theSeed ) ) {
339     std::vector<unsigned long> v;
340     unsigned long uu;
341     for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) {
342       is >> uu;
343       if (!is) {
344         is.clear(std::ios::badbit | is.rdstate());
345         std::cerr << "\nJamesRandom state (vector) description improper."
346     << "\ngetState() has failed."
347          << "\nInput stream is probably mispositioned now." << std::endl;
348         return is;
349       }
350       v.push_back(uu);
351     }
352     getState(v);
353     return (is);
354   }
355 
356 //  is >> theSeed;  Removed, encompassed by possibleKeywordInput()
357 
358   int ipos, jpos;
359   char   endMarker [MarkerLen];
360   for (int i=0; i<97; ++i) {
361      is >> u[i];
362   }
363   is >> c; is >> cd; is >> cm;
364   is >> jpos;
365   is >> std::ws;
366   is.width(MarkerLen);
367   is >> endMarker;
368   if(strcmp(endMarker,"JamesRandom-end")) {
369      is.clear(std::ios::badbit | is.rdstate());
370      std::cerr << "\nJamesRandom state description incomplete."
371          << "\nInput stream is probably mispositioned now." << std::endl;
372      return is;
373   }
374 
375   ipos = (64+jpos)%97;
376   i97 = ipos;
377   j97 = jpos;
378   return is;
379 }
380 
381 bool HepJamesRandom::get (const std::vector<unsigned long> & v) {
382   if ( (v[0] & 0xffffffffUL) != engineIDulong<HepJamesRandom>()) {
383     std::cerr << 
384       "\nHepJamesRandom get:state vector has wrong ID word - state unchanged\n";
385     return false;
386   }
387   return getState(v);
388 }
389 
390 bool HepJamesRandom::getState (const std::vector<unsigned long> & v) {
391   if (v.size() != VECTOR_STATE_SIZE ) {
392     std::cerr << 
393       "\nHepJamesRandom get:state vector has wrong length - state unchanged\n";
394     return false;
395   }
396   std::vector<unsigned long> t(2);
397   for (int i=0; i<97; ++i) {
398     t[0] = v[2*i+1]; t[1] = v[2*i+2];
399     u[i] = DoubConv::longs2double(t);
400   }
401   t[0] = v[195]; t[1] = v[196]; c  = DoubConv::longs2double(t);
402   t[0] = v[197]; t[1] = v[198]; cd = DoubConv::longs2double(t);
403   t[0] = v[199]; t[1] = v[200]; cm = DoubConv::longs2double(t);
404   j97  = (int)v[201];
405   i97  = (64+j97)%97; 
406   return true;
407 }
408 
409 }  // namespace CLHEP
410