Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/externals/clhep/src/DualRand.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/DualRand.cc (Version 11.3.0) and /externals/clhep/src/DualRand.cc (Version 9.5.p2)


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