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

TripleRand.cc
Go to the documentation of this file.
1 // $Id: TripleRand.cc,v 1.6 2010/06/16 17:24:53 garren Exp $
2 // -*- C++ -*-
3 //
4 // -----------------------------------------------------------------------
5 // Hep Random
6 // --- TripleRand ---
7 // class implementation file
8 // -----------------------------------------------------------------------
9 // A canopy pseudo-random number generator. Using the Tausworthe
10 // exclusive-or shift register, a simple Integer Coungruence generator, and
11 // the Hurd 288 total bit shift register, all XOR'd with each other, we
12 // provide an engine that should be a fairly good "mother" generator.
13 // =======================================================================
14 // Ken Smith - Initial draft started: 23rd Jul 1998
15 // - Added conversion operators: 6th Aug 1998
16 // J. Marraffino - Added some explicit casts to deal with
17 // machines where sizeof(int) != sizeof(long) 22 Aug 1998
18 // M. Fischler - Modified use of the various exponents of 2
19 // to avoid per-instance space overhead and
20 // correct the rounding procedure 15 Sep 1998
21 // - modify constructors so that no sequence can
22 // ever accidentally be produced by differnt
23 // seeds, and so that exxcept for the default
24 // constructor, the seed fully determines the
25 // sequence. 15 Sep 1998
26 // M. Fischler - Eliminated dependency on hepString 13 May 1999
27 // M. Fischler - Eliminated Taus() and Cong() which were
28 // causing spurions warnings on SUN CC 27 May 1999
29 // M. Fischler - Put endl at end of puts of Tausworth and 10 Apr 2001
30 // integerCong
31 // M. Fischler - In restore, checkFile for file not found 03 Dec 2004
32 // M. Fischler - Methods put, get for instance save/restore 12/8/04
33 // M. Fischler - split get() into tag validation and
34 // getState() for anonymous restores 12/27/04
35 // M. Fischler - put/get for vectors of ulongs 3/14/05
36 // M. Fischler - State-saving using only ints, for portability 4/12/05
37 //
38 // =======================================================================
39 
40 #include "CLHEP/Random/TripleRand.h"
41 #include "CLHEP/Random/defs.h"
42 #include "CLHEP/Random/engineIDulong.h"
43 #include <string.h> // for strcmp
44 
45 namespace CLHEP {
46 
47 static const int MarkerLen = 64; // Enough room to hold a begin or end marker.
48 
49 //********************************************************************
50 // TripleRand
51 //********************************************************************
52 
53 // Number of instances with automatic seed selection
54 int TripleRand::numEngines = 0;
55 
56 std::string TripleRand::name() const {return "TripleRand";}
57 
59 : HepRandomEngine(),
60  tausworthe (1234567 + numEngines + 175321),
61  integerCong(69607 * tausworthe + 54329, numEngines),
62  hurd(19781127 + integerCong)
63 {
64  theSeed = 1234567;
65  ++numEngines;
66 }
67 
69 : HepRandomEngine(),
70  tausworthe ((unsigned int)seed + 175321),
71  integerCong(69607 * tausworthe + 54329, 1313),
72  hurd(19781127 + integerCong)
73 {
74  theSeed = seed;
75 }
76 
77 TripleRand::TripleRand(std::istream & is)
79 {
80  is >> *this;
81 }
82 
83 TripleRand::TripleRand(int rowIndex, int colIndex)
84 : HepRandomEngine(),
85  tausworthe (rowIndex + numEngines * colIndex + 175321),
86  integerCong(69607 * tausworthe + 54329, 19),
87  hurd(19781127 + integerCong)
88 {
89  theSeed = rowIndex;
90 }
91 
93 
94 double TripleRand::flat() {
95  unsigned int ic ( integerCong );
96  unsigned int t ( tausworthe );
97  unsigned int h ( hurd );
98  return ( (t ^ ic ^ h) * twoToMinus_32() + // most significant part
99  (h >> 11) * twoToMinus_53() + // fill in remaining bits
100  nearlyTwoToMinus_54() // make sure non-zero
101  );
102 }
103 
104 void TripleRand::flatArray(const int size, double* vect) {
105  for (int i = 0; i < size; ++i) {
106  vect[i] = flat();
107  }
108 }
109 
110 void TripleRand::setSeed(long seed, int) {
111  theSeed = seed;
112  tausworthe = Tausworthe((unsigned int)seed + numEngines + 175321);
113  integerCong = IntegerCong(69607 * tausworthe + 54329, numEngines);
114  hurd = Hurd288Engine( 19781127 + integerCong );
115 }
116 
117 void TripleRand::setSeeds(const long * seeds, int) {
118  setSeed(seeds ? *seeds : 1234567, 0);
119  theSeeds = seeds;
120 }
121 
122 void TripleRand::saveStatus(const char filename[]) const {
123  std::ofstream outFile(filename, std::ios::out);
124  if (!outFile.bad()) {
125  outFile << "Uvec\n";
126  std::vector<unsigned long> v = put();
127  #ifdef TRACE_IO
128  std::cout << "Result of v = put() is:\n";
129  #endif
130  for (unsigned int i=0; i<v.size(); ++i) {
131  outFile << v[i] << "\n";
132  #ifdef TRACE_IO
133  std::cout << v[i] << " ";
134  if (i%6==0) std::cout << "\n";
135  #endif
136  }
137  #ifdef TRACE_IO
138  std::cout << "\n";
139  #endif
140  }
141 #ifdef REMOVED
142  outFile << std::setprecision(20) << theSeed << " ";
143  tausworthe.put ( outFile );
144  integerCong.put( outFile);
145  outFile << ConstHurd() << std::endl;
146 #endif
147 }
148 
149 void TripleRand::restoreStatus(const char filename[]) {
150  std::ifstream inFile(filename, std::ios::in);
151  if (!checkFile ( inFile, filename, engineName(), "restoreStatus" )) {
152  std::cerr << " -- Engine state remains unchanged\n";
153  return;
154  }
155  if ( possibleKeywordInput ( inFile, "Uvec", theSeed ) ) {
156  std::vector<unsigned long> v;
157  unsigned long xin;
158  for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) {
159  inFile >> xin;
160  #ifdef TRACE_IO
161  std::cout << "ivec = " << ivec << " xin = " << xin << " ";
162  if (ivec%3 == 0) std::cout << "\n";
163  #endif
164  if (!inFile) {
165  inFile.clear(std::ios::badbit | inFile.rdstate());
166  std::cerr << "\nTripleRand state (vector) description improper."
167  << "\nrestoreStatus has failed."
168  << "\nInput stream is probably mispositioned now." << std::endl;
169  return;
170  }
171  v.push_back(xin);
172  }
173  getState(v);
174  return;
175  }
176 
177  if (!inFile.bad()) {
178 // inFile >> theSeed; removed -- encompased by possibleKeywordInput
179  tausworthe.get ( inFile );
180  integerCong.get( inFile );
181  inFile >> Hurd();
182  }
183 }
184 
186  std::cout << std::setprecision(20) << std::endl;
187  std::cout << "-------- TripleRand engine status ---------"
188  << std::endl;
189  std::cout << "Initial seed = " << theSeed << std::endl;
190  std::cout << "Tausworthe generator = " << std::endl;
191  tausworthe.put( std::cout );
192  std::cout << "IntegerCong generator = " << std::endl;
193  integerCong.put( std::cout );
194  std::cout << "Hurd288Engine generator= " << std::endl << ConstHurd();
195  std::cout << std::endl << "-----------------------------------------"
196  << std::endl;
197 }
198 
199 TripleRand::operator float() {
200  return (float)
201  ( ( integerCong ^ tausworthe ^ (unsigned int)hurd ) * twoToMinus_32()
202  + nearlyTwoToMinus_54() );
203  // make sure non-zero!
204 }
205 
206 TripleRand::operator unsigned int() {
207  return integerCong ^ tausworthe ^ (unsigned int)hurd;
208 }
209 
210 Hurd288Engine & TripleRand::Hurd() { return hurd; }
211 
212 const Hurd288Engine & TripleRand::ConstHurd() const
213  { return hurd; }
214 
215 std::ostream & TripleRand::put (std::ostream & os ) const {
216  char beginMarker[] = "TripleRand-begin";
217  os << beginMarker << "\nUvec\n";
218  std::vector<unsigned long> v = put();
219  for (unsigned int i=0; i<v.size(); ++i) {
220  os << v[i] << "\n";
221  }
222  return os;
223 #ifdef REMOVED
224  char endMarker[] = "TripleRand-end";
225  int pr=os.precision(20);
226  os << " " << beginMarker << "\n";
227  os << theSeed << "\n";
228  tausworthe.put( os );
229  integerCong.put( os );
230  os << ConstHurd();
231  os << " " << endMarker << "\n";
232  os.precision(pr);
233  return os;
234 #endif
235 }
236 
237 std::vector<unsigned long> TripleRand::put () const {
238  std::vector<unsigned long> v;
239  v.push_back (engineIDulong<TripleRand>());
240  tausworthe.put(v);
241  integerCong.put(v);
242  std::vector<unsigned long> vHurd = hurd.put();
243  for (unsigned int i = 0; i < vHurd.size(); ++i) {
244  v.push_back (vHurd[i]);
245  }
246  return v;
247 }
248 
249 std::istream & TripleRand::get (std::istream & is) {
250  char beginMarker [MarkerLen];
251  is >> std::ws;
252  is.width(MarkerLen); // causes the next read to the char* to be <=
253  // that many bytes, INCLUDING A TERMINATION \0
254  // (Stroustrup, section 21.3.2)
255  is >> beginMarker;
256  if (strcmp(beginMarker,"TripleRand-begin")) {
257  is.clear(std::ios::badbit | is.rdstate());
258  std::cerr << "\nInput mispositioned or"
259  << "\nTripleRand state description missing or"
260  << "\nwrong engine type found." << std::endl;
261  return is;
262  }
263  return getState(is);
264 }
265 
266 std::string TripleRand::beginTag ( ) {
267  return "TripleRand-begin";
268 }
269 
270 std::istream & TripleRand::getState (std::istream & is) {
271  if ( possibleKeywordInput ( is, "Uvec", theSeed ) ) {
272  std::vector<unsigned long> v;
273  unsigned long uu;
274  for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) {
275  is >> uu;
276  if (!is) {
277  is.clear(std::ios::badbit | is.rdstate());
278  std::cerr << "\nTripleRand state (vector) description improper."
279  << "\ngetState() has failed."
280  << "\nInput stream is probably mispositioned now." << std::endl;
281  return is;
282  }
283  v.push_back(uu);
284  }
285  getState(v);
286  return (is);
287  }
288 
289 // is >> theSeed; Removed, encompassed by possibleKeywordInput()
290 
291  char endMarker [MarkerLen];
292  tausworthe.get( is );
293  integerCong.get( is );
294  is >> Hurd();
295  is >> std::ws;
296  is.width(MarkerLen);
297  is >> endMarker;
298  if (strcmp(endMarker,"TripleRand-end")) {
299  is.clear(std::ios::badbit | is.rdstate());
300  std::cerr << "\nTripleRand state description incomplete."
301  << "\nInput stream is probably mispositioned now." << std::endl;
302  return is;
303  }
304  return is;
305 }
306 
307 bool TripleRand::get (const std::vector<unsigned long> & v) {
308  if ((v[0] & 0xffffffffUL) != engineIDulong<TripleRand>()) {
309  std::cerr <<
310  "\nTripleRand get:state vector has wrong ID word - state unchanged\n";
311  return false;
312  }
313  if (v.size() != VECTOR_STATE_SIZE) {
314  std::cerr << "\nTripleRand get:state vector has wrong size: "
315  << v.size() << " - state unchanged\n";
316  return false;
317  }
318  return getState(v);
319 }
320 
321 bool TripleRand::getState (const std::vector<unsigned long> & v) {
322  std::vector<unsigned long>::const_iterator iv = v.begin()+1;
323  if (!tausworthe.get(iv)) return false;
324  if (!integerCong.get(iv)) return false;
325  std::vector<unsigned long> vHurd;
326  while (iv != v.end()) {
327  vHurd.push_back(*iv++);
328  }
329  if (!hurd.get(vHurd)) {
330  std::cerr <<
331  "\nTripleRand get from vector: problem getting the hurd sub-engine state\n";
332  return false;
333  }
334  return true;
335 }
336 
337 //********************************************************************
338 // Tausworthe
339 //********************************************************************
340 
341 TripleRand::Tausworthe::Tausworthe() {
342  words[0] = 1234567;
343  for (wordIndex = 1; wordIndex < 4; ++wordIndex) {
344  words[wordIndex] = 69607 * words[wordIndex-1] + 54329;
345  }
346 }
347 
348 TripleRand::Tausworthe::Tausworthe(unsigned int seed) {
349  words[0] = seed;
350  for (wordIndex = 1; wordIndex < 4; ++wordIndex) {
351  words[wordIndex] = 69607 * words[wordIndex-1] + 54329;
352  }
353 }
354 
355 TripleRand::Tausworthe::operator unsigned int() {
356 
357 // Mathematically: Consider a sequence of bits b[n]. Repeatedly form
358 // b[0]' = b[127] ^ b[97]; b[n]' = b[n-1]. This sequence will have a very
359 // long period (2**127-1 according to Tausworthe's work).
360 
361 // The actual method used relies on the fact that the operations needed to
362 // form bit 0 for up to 96 iterations never depend on the results of the
363 // previous ones. So you can actually compute many bits at once. In fact
364 // you can compute 32 at once -- despite 127 - 97 < 32 -- but 24 was used in
365 // the method used in Canopy, where they wanted only single-precision float
366 // randoms. I will do 32 here.
367 
368 // When you do it this way, this looks disturbingly like the dread lagged XOR
369 // Fibonacci. And indeed, it is a lagged Fibonacii, F(4,3, op) with the op
370 // being the XOR of a combination of shifts of the two numbers. Although
371 // Tausworthe asserted excellent properties, I would be scared to death.
372 // However, the shifting and bit swapping really does randomize this in a
373 // serious way.
374 
375 // Statements have been made to the effect that shift register sequences fail
376 // the parking lot test because they achieve randomness by multiple foldings,
377 // and this produces a characteristic pattern. We observe that in this
378 // specific algorithm, which has a fairly long lever arm, the foldings become
379 // effectively random. This is evidenced by the fact that the generator
380 // does pass the Diehard tests, including the parking lot test.
381 
382 // To avoid shuffling of variables in memory, you either have to use circular
383 // pointers (and those give you ifs, which are also costly) or compute at least
384 // a few iterations at once. We do the latter. Although there is a possible
385 // trade of room for more speed, by computing and saving 256 instead of 128
386 // bits at once, I will stop at this level of optimization.
387 
388 // To remind: Each (32-bit) step takes the XOR of bits [127-96] with bits
389 // [95-64] and places it in bits [0-31]. But in the first step, we designate
390 // word0 as bits [0-31], in the second step, word 1 (since the bits it holds
391 // will no longer be needed), then word 2, then word 3. After this, the
392 // stream contains 128 random bits which we will use as 4 valid 32-bit
393 // random numbers.
394 
395 // Thus at the start of the first step, word[0] contains the newest (used)
396 // 32-bit random, and word[3] the oldest. After four steps, word[0] again
397 // contains the newest (now unused) random, and word[3] the oldest.
398 // Bit 0 of word[0] is logically the newest bit, and bit 31 of word[3]
399 // the oldest.
400 
401  if (wordIndex <= 0) {
402  for (wordIndex = 0; wordIndex < 4; ++wordIndex) {
403  words[wordIndex] = ( (words[(wordIndex+1) & 3] << 1 ) |
404  (words[wordIndex] >> 31) )
405  ^ ( (words[(wordIndex+1) & 3] << 31) |
406  (words[wordIndex] >> 1) );
407  }
408  }
409  return words[--wordIndex] & 0xffffffff;
410 }
411 
412 void TripleRand::Tausworthe::put( std::ostream & os ) const {
413  char beginMarker[] = "Tausworthe-begin";
414  char endMarker[] = "Tausworthe-end";
415 
416  int pr=os.precision(20);
417  os << " " << beginMarker << " ";
418  os << std::setprecision(20);
419  for (int i = 0; i < 4; ++i) {
420  os << words[i] << " ";
421  }
422  os << wordIndex;
423  os << " " << endMarker << " ";
424  os << std::endl;
425  os.precision(pr);
426 }
427 
428 void TripleRand::Tausworthe::put(std::vector<unsigned long> & v) const {
429  for (int i = 0; i < 4; ++i) {
430  v.push_back(static_cast<unsigned long>(words[i]));
431  }
432  v.push_back(static_cast<unsigned long>(wordIndex));
433 }
434 
435 void TripleRand::Tausworthe::get( std::istream & is ) {
436  char beginMarker [MarkerLen];
437  char endMarker [MarkerLen];
438 
439  is >> std::ws;
440  is.width(MarkerLen);
441  is >> beginMarker;
442  if (strcmp(beginMarker,"Tausworthe-begin")) {
443  is.clear(std::ios::badbit | is.rdstate());
444  std::cerr << "\nInput mispositioned or"
445  << "\nTausworthe state description missing or"
446  << "\nwrong engine type found." << std::endl;
447  }
448  for (int i = 0; i < 4; ++i) {
449  is >> words[i];
450  }
451  is >> wordIndex;
452  is >> std::ws;
453  is.width(MarkerLen);
454  is >> endMarker;
455  if (strcmp(endMarker,"Tausworthe-end")) {
456  is.clear(std::ios::badbit | is.rdstate());
457  std::cerr << "\nTausworthe state description incomplete."
458  << "\nInput stream is probably mispositioned now." << std::endl;
459  }
460 }
461 
462 bool
463 TripleRand::Tausworthe::get(std::vector<unsigned long>::const_iterator & iv){
464  for (int i = 0; i < 4; ++i) {
465  words[i] = *iv++;
466  }
467  wordIndex = *iv++;
468  return true;
469 }
470 
471 //********************************************************************
472 // IntegerCong
473 //********************************************************************
474 
475 TripleRand::IntegerCong::IntegerCong()
476 : state((unsigned int)3758656018U),
477  multiplier(66565),
478  addend(12341)
479 {
480 }
481 
482 TripleRand::IntegerCong::IntegerCong(unsigned int seed, int streamNumber)
483 : state(seed),
484  multiplier(65536 + 1024 + 5 + (8 * 1017 * streamNumber)),
485  addend(12341)
486 {
487  // As to the multiplier, the following comment was made:
488  // We want our multipliers larger than 2^16, and equal to
489  // 1 mod 4 (for max. period), but not equal to 1 mod 8
490  // (for max. potency -- the smaller and higher dimension the
491  // stripes, the better)
492 
493  // All of these will have fairly long periods. Depending on the choice
494  // of stream number, some of these will be quite decent when considered
495  // as independent random engines, while others will be poor. Thus these
496  // should not be used as stand-alone engines; but when combined with a
497  // generator known to be good, they allow creation of millions of good
498  // independent streams, without fear of two streams accidentally hitting
499  // nearby places in the good random sequence.
500 }
501 
502 TripleRand::IntegerCong::operator unsigned int() {
503  return state = (state * multiplier + addend) & 0xffffffff;
504 }
505 
506 void TripleRand::IntegerCong::put( std::ostream & os ) const {
507  char beginMarker[] = "IntegerCong-begin";
508  char endMarker[] = "IntegerCong-end";
509 
510  int pr=os.precision(20);
511  os << " " << beginMarker << " ";
512  os << state << " " << multiplier << " " << addend;
513  os << " " << endMarker << " ";
514  os << std::endl;
515  os.precision(pr);
516 }
517 
518 void TripleRand::IntegerCong::put(std::vector<unsigned long> & v) const {
519  v.push_back(static_cast<unsigned long>(state));
520  v.push_back(static_cast<unsigned long>(multiplier));
521  v.push_back(static_cast<unsigned long>(addend));
522 }
523 
524 void TripleRand::IntegerCong::get( std::istream & is ) {
525  char beginMarker [MarkerLen];
526  char endMarker [MarkerLen];
527 
528  is >> std::ws;
529  is.width(MarkerLen);
530  is >> beginMarker;
531  if (strcmp(beginMarker,"IntegerCong-begin")) {
532  is.clear(std::ios::badbit | is.rdstate());
533  std::cerr << "\nInput mispositioned or"
534  << "\nIntegerCong state description missing or"
535  << "\nwrong engine type found." << std::endl;
536  }
537  is >> state >> multiplier >> addend;
538  is >> std::ws;
539  is.width(MarkerLen);
540  is >> endMarker;
541  if (strcmp(endMarker,"IntegerCong-end")) {
542  is.clear(std::ios::badbit | is.rdstate());
543  std::cerr << "\nIntegerCong state description incomplete."
544  << "\nInput stream is probably mispositioned now." << std::endl;
545  }
546 }
547 
548 bool
549 TripleRand::IntegerCong::get(std::vector<unsigned long>::const_iterator & iv) {
550  state = *iv++;
551  multiplier = *iv++;
552  addend = *iv++;
553  return true;
554 }
555 
556 } // namespace CLHEP
static double twoToMinus_32()
static double twoToMinus_53()
static double nearlyTwoToMinus_54()
static bool checkFile(std::istream &file, const std::string &filename, const std::string &classname, const std::string &methodname)
Definition: RandomEngine.cc:46
virtual std::ostream & put(std::ostream &os) const
virtual std::istream & get(std::istream &is)
virtual std::istream & get(std::istream &is)
Definition: TripleRand.cc:249
virtual std::istream & getState(std::istream &is)
Definition: TripleRand.cc:270
std::vector< unsigned long > put() const
Definition: TripleRand.cc:237
void flatArray(const int size, double *vect)
Definition: TripleRand.cc:104
virtual ~TripleRand()
Definition: TripleRand.cc:92
std::string name() const
Definition: TripleRand.cc:56
void saveStatus(const char filename[]="TripleRand.conf") const
Definition: TripleRand.cc:122
void setSeed(long seed, int)
Definition: TripleRand.cc:110
void showStatus() const
Definition: TripleRand.cc:185
void restoreStatus(const char filename[]="TripleRand.conf")
Definition: TripleRand.cc:149
void setSeeds(const long *seeds, int)
Definition: TripleRand.cc:117
static std::string engineName()
static const unsigned int VECTOR_STATE_SIZE
static std::string beginTag()
Definition: TripleRand.cc:266
bool possibleKeywordInput(IS &is, const std::string &key, T &t)