CLHEP VERSION Reference Documentation
   
CLHEP Home Page     CLHEP Documentation     CLHEP Bug Reports

RandEngine.cc
Go to the documentation of this file.
1// $Id: RandEngine.cc,v 1.8 2010/06/16 17:24:53 garren Exp $
2// -*- C++ -*-
3//
4// -----------------------------------------------------------------------
5// HEP Random
6// --- RandEngine ---
7// class implementation file
8// -----------------------------------------------------------------------
9// This file is part of Geant4 (simulation toolkit for HEP).
10
11// =======================================================================
12// Gabriele Cosmo - Created: 5th September 1995
13// - Minor corrections: 31st October 1996
14// - Added methods for engine status: 19th November 1996
15// - mx is initialised to RAND_MAX: 2nd April 1997
16// - Fixed bug in setSeeds(): 15th September 1997
17// - Private copy constructor and operator=: 26th Feb 1998
18// J.Marraffino - Added stream operators and related constructor.
19// Added automatic seed selection from seed table and
20// engine counter. Removed RAND_MAX and replaced by
21// std::pow(0.5,32.). Flat() returns now 32 bits values
22// obtained by concatenation: 15th Feb 1998
23// Ken Smith - Added conversion operators: 6th Aug 1998
24// J. Marraffino - Remove dependence on hepString class 13 May 1999
25// M. Fischler - Rapaired bug that in flat() that relied on rand() to
26// deliver 15-bit results. This bug was causing flat()
27// on Linux systems to yield randoms with mean of 5/8(!)
28// - Modified algorithm such that on systems with 31-bit rand()
29// it will call rand() only once instead of twice. Feb 2004
30// M. Fischler - Modified the general-case template for RandEngineBuilder
31// such that when RAND_MAX is an unexpected value the routine
32// will still deliver a sensible flat() random.
33// M. Fischler - Methods for distrib. instance save/restore 12/8/04
34// M. Fischler - split get() into tag validation and
35// getState() for anonymous restores 12/27/04
36// M. Fischler - put/get for vectors of ulongs 3/14/05
37// M. Fischler - State-saving using only ints, for portability 4/12/05
38//
39// =======================================================================
40
41#include "CLHEP/Random/defs.h"
42#include "CLHEP/Random/RandEngine.h"
43#include "CLHEP/Random/Random.h"
44#include "CLHEP/Random/engineIDulong.h"
45#include <string.h> // for strcmp
46#include <cstdlib> // for int()
47
48namespace CLHEP {
49
50#ifdef NOTDEF
51// The way to test for proper behavior of the RandEngineBuilder
52// for arbitrary RAND_MAX, on a system where RAND_MAX is some
53// fixed specialized value and rand() behaves accordingly, is
54// to set up a fake RAND_MAX and a fake version of rand()
55// by enabling this block.
56#undef RAND_MAX
57#define RAND_MAX 9382956
58#include "CLHEP/Random/MTwistEngine.h"
59#include "CLHEP/Random/RandFlat.h"
60MTwistEngine * fakeFlat = new MTwistEngine;
61RandFlat rflat (fakeFlat, 0, RAND_MAX+1);
62int rand() { return (int)rflat(); }
63#endif
64
65static const int MarkerLen = 64; // Enough room to hold a begin or end marker.
66
67// number of instances with automatic seed selection
68int RandEngine::numEngines = 0;
69
70// Maximum index into the seed table
71int RandEngine::maxIndex = 215;
72
73std::string RandEngine::name() const {return "RandEngine";}
74
77{
78 setSeed(seed,0);
79 setSeeds(&theSeed,0);
80 seq = 0;
81}
82
83RandEngine::RandEngine(int rowIndex, int colIndex)
85{
86 long seeds[2];
87 long seed;
88
89 int cycle = std::abs(int(rowIndex/maxIndex));
90 int row = std::abs(int(rowIndex%maxIndex));
91 int col = std::abs(int(colIndex%2));
92 long mask = ((cycle & 0x000007ff) << 20 );
93 HepRandom::getTheTableSeeds( seeds, row );
94 seed = (seeds[col])^mask;
95 setSeed(seed,0);
96 setSeeds(&theSeed,0);
97 seq = 0;
98}
99
102{
103 long seeds[2];
104 long seed;
105 int cycle,curIndex;
106
107 cycle = std::abs(int(numEngines/maxIndex));
108 curIndex = std::abs(int(numEngines%maxIndex));
109 numEngines += 1;
110 long mask = ((cycle & 0x007fffff) << 8);
111 HepRandom::getTheTableSeeds( seeds, curIndex );
112 seed = seeds[0]^mask;
113 setSeed(seed,0);
114 setSeeds(&theSeed,0);
115 seq = 0;
116}
117
118RandEngine::RandEngine(std::istream& is)
120{
121 is >> *this;
122}
123
125
126void RandEngine::setSeed(long seed, int)
127{
128 theSeed = seed;
129 srand( int(seed) );
130 seq = 0;
131}
132
133void RandEngine::setSeeds(const long* seeds, int)
134{
135 setSeed(seeds ? *seeds : 19780503L, 0);
136 theSeeds = seeds;
137}
138
139void RandEngine::saveStatus( const char filename[] ) const
140{
141 std::ofstream outFile( filename, std::ios::out ) ;
142
143 if (!outFile.bad()) {
144 outFile << "Uvec\n";
145 std::vector<unsigned long> v = put();
146 #ifdef TRACE_IO
147 std::cout << "Result of v = put() is:\n";
148 #endif
149 for (unsigned int i=0; i<v.size(); ++i) {
150 outFile << v[i] << "\n";
151 #ifdef TRACE_IO
152 std::cout << v[i] << " ";
153 if (i%6==0) std::cout << "\n";
154 #endif
155 }
156 #ifdef TRACE_IO
157 std::cout << "\n";
158 #endif
159 }
160#ifdef REMOVED
161 if (!outFile.bad()) {
162 outFile << theSeed << std::endl;
163 outFile << seq << std::endl;
164 }
165#endif
166}
167
168void RandEngine::restoreStatus( const char filename[] )
169{
170 // The only way to restore the status of RandEngine is to
171 // keep track of the number of shooted random sequences, reset
172 // the engine and re-shoot them again. The Rand algorithm does
173 // not provide any way of getting its internal status.
174
175 std::ifstream inFile( filename, std::ios::in);
176 if (!checkFile ( inFile, filename, engineName(), "restoreStatus" )) {
177 std::cout << " -- Engine state remains unchanged\n";
178 return;
179 }
180 if ( possibleKeywordInput ( inFile, "Uvec", theSeed ) ) {
181 std::vector<unsigned long> v;
182 unsigned long xin;
183 for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) {
184 inFile >> xin;
185 #ifdef TRACE_IO
186 std::cout << "ivec = " << ivec << " xin = " << xin << " ";
187 if (ivec%3 == 0) std::cout << "\n";
188 #endif
189 if (!inFile) {
190 inFile.clear(std::ios::badbit | inFile.rdstate());
191 std::cerr << "\nRandEngine state (vector) description improper."
192 << "\nrestoreStatus has failed."
193 << "\nInput stream is probably mispositioned now." << std::endl;
194 return;
195 }
196 v.push_back(xin);
197 }
198 getState(v);
199 return;
200 }
201
202 long count;
203
204 if (!inFile.bad() && !inFile.eof()) {
205// inFile >> theSeed; removed -- encompased by possibleKeywordInput
206 inFile >> count;
207 setSeed(theSeed,0);
208 seq = 0;
209 while (seq < count) flat();
210 }
211}
212
214{
215 std::cout << std::endl;
216 std::cout << "---------- Rand engine status ----------" << std::endl;
217 std::cout << " Initial seed = " << theSeed << std::endl;
218 std::cout << " Shooted sequences = " << seq << std::endl;
219 std::cout << "----------------------------------------" << std::endl;
220}
221
222// ====================================================
223// Implementation of flat() (along with needed helpers)
224// ====================================================
225
226// Here we set up such that **at compile time**, the compiler decides based on
227// RAND_MAX how to form a random double with 32 leading random bits, using
228// one or two calls to rand(). Some of the lowest order bits of 32 are allowed
229// to be as weak as mere XORs of some higher bits, but not to be always fixed.
230//
231// The method decision is made at compile time, rather than using a run-time
232// if on the value of RAND_MAX. Although it is possible to cope with arbitrary
233// values of RAND_MAX of the form 2**N-1, with the same efficiency, the
234// template techniques needed would look mysterious and perhaps over-stress
235// some older compilers. We therefore only treat RAND_MAX = 2**15-1 (as on
236// most older systems) and 2**31-1 (as on the later Linux systems) as special
237// and super-efficient cases. We detect any different value, and use an
238// algorithm which is correct even if RAND_MAX is not one less than a power
239// of 2.
240
241 template <int> struct RandEngineBuilder { // RAND_MAX any arbitrary value
242 static unsigned int thirtyTwoRandomBits(long& seq) {
243
244 static bool prepared = false;
245 static unsigned int iT;
246 static unsigned int iK;
247 static unsigned int iS;
248 static int iN;
249 static double fS;
250 static double fT;
251
252 if ( (RAND_MAX >> 31) > 0 )
253 {
254 // Here, we know that integer arithmetic is 64 bits.
255 if ( !prepared ) {
256 iS = (unsigned long)RAND_MAX + 1;
257 iK = 1;
258// int StoK = S;
259 int StoK = iS;
260 // The two statements below are equivalent, but some compilers
261 // are getting too smart and complain about the first statement.
262 //if ( (RAND_MAX >> 32) == 0) {
263 if( (unsigned long) (RAND_MAX) <= (( (1uL) << 31 ) - 1 ) ) {
264 iK = 2;
265// StoK = S*S;
266 StoK = iS*iS;
267 }
268 int m;
269 for ( m = 0; m < 64; ++m ) {
270 StoK >>= 1;
271 if (StoK == 0) break;
272 }
273 iT = 1 << m;
274 prepared = true;
275 }
276 int v = 0;
277 do {
278 for ( int i = 0; i < iK; ++i ) {
279 v = iS*v+rand(); ++seq;
280 }
281 } while (v < iT);
282 return v & 0xFFFFFFFF;
283
284 }
285
286 else if ( (RAND_MAX >> 26) == 0 )
287 {
288 // Here, we know it is safe to work in terms of doubles without loss
289 // of precision, but we have no idea how many randoms we will need to
290 // generate 32 bits.
291 if ( !prepared ) {
292 fS = (unsigned long)RAND_MAX + 1;
293 double twoTo32 = std::ldexp(1.0,32);
294 double StoK = fS;
295 for ( iK = 1; StoK < twoTo32; StoK *= fS, iK++ ) { }
296 int m;
297 fT = 1.0;
298 for ( m = 0; m < 64; ++m ) {
299 StoK *= .5;
300 if (StoK < 1.0) break;
301 fT *= 2.0;
302 }
303 prepared = true;
304 }
305 double v = 0;
306 do {
307 for ( int i = 0; i < iK; ++i ) {
308 v = fS*v+rand(); ++seq;
309 }
310 } while (v < fT);
311 return ((unsigned int)v) & 0xFFFFFFFF;
312
313 }
314 else
315 {
316 // Here, we know that 16 random bits are available from each of
317 // two random numbers.
318 if ( !prepared ) {
319 iS = (unsigned long)RAND_MAX + 1;
320 int SshiftN = iS;
321 for (iN = 0; SshiftN > 1; SshiftN >>= 1, iN++) { }
322 iN -= 17;
323 prepared = true;
324 }
325 unsigned int x1, x2;
326 do { x1 = rand(); ++seq;} while (x1 < (1<<16) );
327 do { x2 = rand(); ++seq;} while (x2 < (1<<16) );
328 return x1 | (x2 << 16);
329 }
330
331 }
332};
333
334template <> struct RandEngineBuilder<2147483647> { // RAND_MAX = 2**31 - 1
335 inline static unsigned int thirtyTwoRandomBits(long& seq) {
336 unsigned int x = rand() << 1; ++seq; // bits 31-1
337 x ^= ( (x>>23) ^ (x>>7) ) ^1; // bit 0 (weakly pseudo-random)
338 return x & 0xFFFFFFFF; // mask in case int is 64 bits
339 }
340};
341
342
343template <> struct RandEngineBuilder<32767> { // RAND_MAX = 2**15 - 1
344 inline static unsigned int thirtyTwoRandomBits(long& seq) {
345 unsigned int x = rand() << 17; ++seq; // bits 31-17
346 x ^= rand() << 2; ++seq; // bits 16-2
347 x ^= ( (x>>23) ^ (x>>7) ) ^3; // bits 1-0 (weakly pseudo-random)
348 return x & 0xFFFFFFFF; // mask in case int is 64 bits
349 }
350};
351
353{
354 double r;
356 } while ( r == 0 );
357 return r/4294967296.0;
358}
359
360void RandEngine::flatArray(const int size, double* vect)
361{
362 int i;
363
364 for (i=0; i<size; ++i)
365 vect[i]=flat();
366}
367
368RandEngine::operator unsigned int() {
370}
371
372std::ostream & RandEngine::put ( std::ostream& os ) const
373{
374 char beginMarker[] = "RandEngine-begin";
375 char endMarker[] = "RandEngine-end";
376
377 os << " " << beginMarker << "\n";
378 os << theSeed << " " << seq << " ";
379 os << endMarker << "\n";
380 return os;
381}
382
383std::vector<unsigned long> RandEngine::put () const {
384 std::vector<unsigned long> v;
385 v.push_back (engineIDulong<RandEngine>());
386 v.push_back(static_cast<unsigned long>(theSeed));
387 v.push_back(static_cast<unsigned long>(seq));
388 return v;
389}
390
391std::istream & RandEngine::get ( std::istream& is )
392{
393 // The only way to restore the status of RandEngine is to
394 // keep track of the number of shooted random sequences, reset
395 // the engine and re-shoot them again. The Rand algorithm does
396 // not provide any way of getting its internal status.
397 char beginMarker [MarkerLen];
398 is >> std::ws;
399 is.width(MarkerLen); // causes the next read to the char* to be <=
400 // that many bytes, INCLUDING A TERMINATION \0
401 // (Stroustrup, section 21.3.2)
402 is >> beginMarker;
403 if (strcmp(beginMarker,"RandEngine-begin")) {
404 is.clear(std::ios::badbit | is.rdstate());
405 std::cout << "\nInput stream mispositioned or"
406 << "\nRandEngine state description missing or"
407 << "\nwrong engine type found." << std::endl;
408 return is;
409 }
410 return getState(is);
411}
412
413std::string RandEngine::beginTag ( ) {
414 return "RandEngine-begin";
415}
416
417std::istream & RandEngine::getState ( std::istream& is )
418{
419 if ( possibleKeywordInput ( is, "Uvec", theSeed ) ) {
420 std::vector<unsigned long> v;
421 unsigned long uu;
422 for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) {
423 is >> uu;
424 if (!is) {
425 is.clear(std::ios::badbit | is.rdstate());
426 std::cerr << "\nRandEngine state (vector) description improper."
427 << "\ngetState() has failed."
428 << "\nInput stream is probably mispositioned now." << std::endl;
429 return is;
430 }
431 v.push_back(uu);
432 }
433 getState(v);
434 return (is);
435 }
436
437// is >> theSeed; Removed, encompassed by possibleKeywordInput()
438
439 char endMarker [MarkerLen];
440 long count;
441 is >> count;
442 is >> std::ws;
443 is.width(MarkerLen);
444 is >> endMarker;
445 if (strcmp(endMarker,"RandEngine-end")) {
446 is.clear(std::ios::badbit | is.rdstate());
447 std::cerr << "\nRandEngine state description incomplete."
448 << "\nInput stream is probably mispositioned now." << std::endl;
449 return is;
450 }
451 setSeed(theSeed,0);
452 while (seq < count) flat();
453 return is;
454}
455
456bool RandEngine::get (const std::vector<unsigned long> & v) {
457 if ((v[0] & 0xffffffffUL) != engineIDulong<RandEngine>()) {
458 std::cerr <<
459 "\nRandEngine get:state vector has wrong ID word - state unchanged\n";
460 return false;
461 }
462 return getState(v);
463}
464
465bool RandEngine::getState (const std::vector<unsigned long> & v) {
466 if (v.size() != VECTOR_STATE_SIZE ) {
467 std::cerr <<
468 "\nRandEngine get:state vector has wrong length - state unchanged\n";
469 return false;
470 }
471 theSeed = v[1];
472 int count = v[2];
473 setSeed(theSeed,0);
474 while (seq < count) flat();
475 return true;
476}
477} // namespace CLHEP
static bool checkFile(std::istream &file, const std::string &filename, const std::string &classname, const std::string &methodname)
static void getTheTableSeeds(long *seeds, int index)
Definition Random.cc:152
virtual std::istream & get(std::istream &is)
void flatArray(const int size, double *vect)
virtual std::istream & getState(std::istream &is)
void saveStatus(const char filename[]="Rand.conf") const
static const unsigned int VECTOR_STATE_SIZE
void setSeed(long seed, int dum=0)
static std::string beginTag()
std::vector< unsigned long > put() const
std::string name() const
Definition RandEngine.cc:73
void restoreStatus(const char filename[]="Rand.conf")
virtual ~RandEngine()
void showStatus() const
void setSeeds(const long *seeds, int dum=0)
bool possibleKeywordInput(IS &is, const std::string &key, T &t)
static unsigned int thirtyTwoRandomBits(long &seq)
static unsigned int thirtyTwoRandomBits(long &seq)
static unsigned int thirtyTwoRandomBits(long &seq)