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
45namespace CLHEP {
46
47static 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
54int TripleRand::numEngines = 0;
55
56std::string TripleRand::name() const {return "TripleRand";}
57
60 tausworthe (1234567 + numEngines + 175321),
61 integerCong(69607 * tausworthe + 54329, numEngines),
62 hurd(19781127 + integerCong)
63{
64 theSeed = 1234567;
65 ++numEngines;
66}
67
70 tausworthe ((unsigned int)seed + 175321),
71 integerCong(69607 * tausworthe + 54329, 1313),
72 hurd(19781127 + integerCong)
73{
74 theSeed = seed;
75}
76
77TripleRand::TripleRand(std::istream & is)
79{
80 is >> *this;
81}
82
83TripleRand::TripleRand(int rowIndex, int colIndex)
85 tausworthe (rowIndex + numEngines * colIndex + 175321),
86 integerCong(69607 * tausworthe + 54329, 19),
87 hurd(19781127 + integerCong)
88{
89 theSeed = rowIndex;
90}
91
93
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
104void TripleRand::flatArray(const int size, double* vect) {
105 for (int i = 0; i < size; ++i) {
106 vect[i] = flat();
107 }
108}
109
110void 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
117void TripleRand::setSeeds(const long * seeds, int) {
118 setSeed(seeds ? *seeds : 1234567, 0);
119 theSeeds = seeds;
120}
121
122void 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
149void 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
199TripleRand::operator float() {
200 return (float)
201 ( ( integerCong ^ tausworthe ^ (unsigned int)hurd ) * twoToMinus_32()
202 + nearlyTwoToMinus_54() );
203 // make sure non-zero!
204}
205
206TripleRand::operator unsigned int() {
207 return integerCong ^ tausworthe ^ (unsigned int)hurd;
208}
209
210Hurd288Engine & TripleRand::Hurd() { return hurd; }
211
212const Hurd288Engine & TripleRand::ConstHurd() const
213 { return hurd; }
214
215std::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
237std::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
249std::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
266std::string TripleRand::beginTag ( ) {
267 return "TripleRand-begin";
268}
269
270std::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
307bool 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
321bool 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
341TripleRand::Tausworthe::Tausworthe() {
342 words[0] = 1234567;
343 for (wordIndex = 1; wordIndex < 4; ++wordIndex) {
344 words[wordIndex] = 69607 * words[wordIndex-1] + 54329;
345 }
346}
347
348TripleRand::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
355TripleRand::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
412void 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
428void 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
435void 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
462bool
463TripleRand::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
475TripleRand::IntegerCong::IntegerCong()
476: state((unsigned int)3758656018U),
477 multiplier(66565),
478 addend(12341)
479{
480}
481
482TripleRand::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
502TripleRand::IntegerCong::operator unsigned int() {
503 return state = (state * multiplier + addend) & 0xffffffff;
504}
505
506void 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
518void 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
524void 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
548bool
549TripleRand::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)
virtual std::ostream & put(std::ostream &os) const
virtual std::istream & get(std::istream &is)
virtual std::istream & get(std::istream &is)
virtual std::istream & getState(std::istream &is)
std::vector< unsigned long > put() const
void flatArray(const int size, double *vect)
virtual ~TripleRand()
Definition TripleRand.cc:92
std::string name() const
Definition TripleRand.cc:56
void saveStatus(const char filename[]="TripleRand.conf") const
void setSeed(long seed, int)
void showStatus() const
void restoreStatus(const char filename[]="TripleRand.conf")
void setSeeds(const long *seeds, int)
static const unsigned int VECTOR_STATE_SIZE
static std::string beginTag()
bool possibleKeywordInput(IS &is, const std::string &key, T &t)