| Geant4 Cross Reference |
1 // -*- C++ -*- 1 // -*- C++ -*-
2 // 2 //
3 // ------------------------------------------- 3 // -----------------------------------------------------------------------
4 // HEP Random 4 // HEP Random
5 // --- MTwistEngine --- 5 // --- MTwistEngine ---
6 // class implementation f 6 // class implementation file
7 // ------------------------------------------- 7 // -----------------------------------------------------------------------
8 // A "fast, compact, huge-period generator" ba 8 // A "fast, compact, huge-period generator" based on M. Matsumoto and
9 // T. Nishimura, "Mersenne Twister: A 623-dime 9 // T. Nishimura, "Mersenne Twister: A 623-dimensionally equidistributed
10 // uniform pseudorandom number generator", to 10 // uniform pseudorandom number generator", to appear in ACM Trans. on
11 // Modeling and Computer Simulation. It is a 11 // Modeling and Computer Simulation. It is a twisted GFSR generator
12 // with a Mersenne-prime period of 2^19937-1, 12 // with a Mersenne-prime period of 2^19937-1, uniform on open interval (0,1)
13 // =========================================== 13 // =======================================================================
14 // Ken Smith - Started initial draft: 14 // Ken Smith - Started initial draft: 14th Jul 1998
15 // - Optimized to get std::pow( 15 // - Optimized to get std::pow() out of flat() method
16 // - Added conversion operators 16 // - Added conversion operators: 6th Aug 1998
17 // J. Marraffino - Added some explicit casts 17 // J. Marraffino - Added some explicit casts to deal with
18 // machines where sizeof(int) 18 // machines where sizeof(int) != sizeof(long) 22 Aug 1998
19 // M. Fischler - Modified constructors such 19 // M. Fischler - Modified constructors such that no two
20 // seeds will match sequences, no singl 20 // seeds will match sequences, no single/double
21 // seeds will match, explicit seeds giv 21 // seeds will match, explicit seeds give
22 // determined results, and default will 22 // determined results, and default will not
23 // match any of the above or other defa 23 // match any of the above or other defaults.
24 // - Modified use of the variou 24 // - Modified use of the various exponents of 2
25 // to avoid per-instance spac 25 // to avoid per-instance space overhead and
26 // correct the rounding proce 26 // correct the rounding procedure 16 Sep 1998
27 // J. Marraffino - Remove dependence on hepSt 27 // J. Marraffino - Remove dependence on hepString class 13 May 1999
28 // M. Fischler - In restore, checkFile for 28 // M. Fischler - In restore, checkFile for file not found 03 Dec 2004
29 // M. Fischler - Methods for distrib. insta 29 // M. Fischler - Methods for distrib. instacne save/restore 12/8/04
30 // M. Fischler - split get() into tag valid 30 // M. Fischler - split get() into tag validation and
31 // getState() for anonymous r 31 // getState() for anonymous restores 12/27/04
32 // M. Fischler - put/get for vectors of ulo 32 // M. Fischler - put/get for vectors of ulongs 3/14/05
33 // M. Fischler - State-saving using only in 33 // M. Fischler - State-saving using only ints, for portability 4/12/05
34 // M. Fischler - Improved seeding in setSee 34 // M. Fischler - Improved seeding in setSeed (Savanah bug #17479) 11/15/06
35 // - (Possible optimization - now that th 35 // - (Possible optimization - now that the seeding is improved,
36 // is it still necessary for ctor to "w 36 // is it still necessary for ctor to "warm up" by discarding
37 // 2000 iterations?) 37 // 2000 iterations?)
38 // 38 //
39 // =========================================== 39 // =======================================================================
40 40
41 #include "CLHEP/Random/Random.h" 41 #include "CLHEP/Random/Random.h"
42 #include "CLHEP/Random/MTwistEngine.h" 42 #include "CLHEP/Random/MTwistEngine.h"
43 #include "CLHEP/Random/engineIDulong.h" 43 #include "CLHEP/Random/engineIDulong.h"
44 #include "CLHEP/Utility/atomic_int.h" 44 #include "CLHEP/Utility/atomic_int.h"
45 45
46 #include <atomic> 46 #include <atomic>
47 #include <cmath> 47 #include <cmath>
48 #include <iostream> 48 #include <iostream>
49 #include <string.h> // for strcmp 49 #include <string.h> // for strcmp
50 #include <vector> 50 #include <vector>
51 51
52 namespace CLHEP { 52 namespace CLHEP {
53 53
54 namespace { 54 namespace {
55 // Number of instances with automatic seed s 55 // Number of instances with automatic seed selection
56 CLHEP_ATOMIC_INT_TYPE numberOfEngines(0); 56 CLHEP_ATOMIC_INT_TYPE numberOfEngines(0);
57 57
58 // Maximum index into the seed table 58 // Maximum index into the seed table
59 const int maxIndex = 215; 59 const int maxIndex = 215;
60 } 60 }
61 61
62 static const int MarkerLen = 64; // Enough roo 62 static const int MarkerLen = 64; // Enough room to hold a begin or end marker.
63 63
64 std::string MTwistEngine::name() const {return 64 std::string MTwistEngine::name() const {return "MTwistEngine";}
65 65
66 MTwistEngine::MTwistEngine() 66 MTwistEngine::MTwistEngine()
67 : HepRandomEngine() 67 : HepRandomEngine()
68 { 68 {
69 int numEngines = numberOfEngines++; 69 int numEngines = numberOfEngines++;
70 int cycle = std::abs(int(numEngines/maxIndex 70 int cycle = std::abs(int(numEngines/maxIndex));
71 int curIndex = std::abs(int(numEngines%maxIn 71 int curIndex = std::abs(int(numEngines%maxIndex));
72 long mask = ((cycle & 0x007fffff) << 8); 72 long mask = ((cycle & 0x007fffff) << 8);
73 long seedlist[2]; 73 long seedlist[2];
74 HepRandom::getTheTableSeeds( seedlist, curIn 74 HepRandom::getTheTableSeeds( seedlist, curIndex );
75 seedlist[0] = (seedlist[0])^mask; 75 seedlist[0] = (seedlist[0])^mask;
76 seedlist[1] = 0; 76 seedlist[1] = 0;
77 setSeeds( seedlist, numEngines ); 77 setSeeds( seedlist, numEngines );
78 count624=0; 78 count624=0;
79 79
80 for( int i=0; i < 2000; ++i ) flat(); / 80 for( int i=0; i < 2000; ++i ) flat(); // Warm up just a bit
81 } 81 }
82 82
83 MTwistEngine::MTwistEngine(long seed) 83 MTwistEngine::MTwistEngine(long seed)
84 : HepRandomEngine() 84 : HepRandomEngine()
85 { 85 {
86 long seedlist[2]={seed,17587}; 86 long seedlist[2]={seed,17587};
87 setSeeds( seedlist, 0 ); 87 setSeeds( seedlist, 0 );
88 count624=0; 88 count624=0;
89 for( int i=0; i < 2000; ++i ) flat(); / 89 for( int i=0; i < 2000; ++i ) flat(); // Warm up just a bit
90 } 90 }
91 91
92 MTwistEngine::MTwistEngine(int rowIndex, int c 92 MTwistEngine::MTwistEngine(int rowIndex, int colIndex)
93 : HepRandomEngine() 93 : HepRandomEngine()
94 { 94 {
95 int cycle = std::abs(int(rowIndex/maxIndex)) 95 int cycle = std::abs(int(rowIndex/maxIndex));
96 int row = std::abs(int(rowIndex%maxIndex)); 96 int row = std::abs(int(rowIndex%maxIndex));
97 int col = std::abs(int(colIndex%2)); 97 int col = std::abs(int(colIndex%2));
98 long mask = (( cycle & 0x000007ff ) << 20 ); 98 long mask = (( cycle & 0x000007ff ) << 20 );
99 long seedlist[2]; 99 long seedlist[2];
100 HepRandom::getTheTableSeeds( seedlist, row ) 100 HepRandom::getTheTableSeeds( seedlist, row );
101 seedlist[0] = (seedlist[col])^mask; 101 seedlist[0] = (seedlist[col])^mask;
102 seedlist[1] = 690691; 102 seedlist[1] = 690691;
103 setSeeds(seedlist, 4444772); 103 setSeeds(seedlist, 4444772);
104 count624=0; 104 count624=0;
105 for( int i=0; i < 2000; ++i ) flat(); / 105 for( int i=0; i < 2000; ++i ) flat(); // Warm up just a bit
106 } 106 }
107 107
108 MTwistEngine::MTwistEngine( std::istream& is ) 108 MTwistEngine::MTwistEngine( std::istream& is )
109 : HepRandomEngine() 109 : HepRandomEngine()
110 { 110 {
111 is >> *this; 111 is >> *this;
112 } 112 }
113 113
114 MTwistEngine::~MTwistEngine() {} 114 MTwistEngine::~MTwistEngine() {}
115 115
116 double MTwistEngine::flat() { 116 double MTwistEngine::flat() {
117 unsigned int y; 117 unsigned int y;
118 118
119 if( count624 >= N ) { 119 if( count624 >= N ) {
120 int i; 120 int i;
121 121
122 for( i=0; i < NminusM; ++i ) { 122 for( i=0; i < NminusM; ++i ) {
123 y = (mt[i] & 0x80000000) | (mt[i+1] & 0x 123 y = (mt[i] & 0x80000000) | (mt[i+1] & 0x7fffffff);
124 mt[i] = mt[i+M] ^ (y >> 1) ^ ((y & 124 mt[i] = mt[i+M] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
125 } 125 }
126 126
127 for( ; i < N-1 ; ++i ) { 127 for( ; i < N-1 ; ++i ) {
128 y = (mt[i] & 0x80000000) | (mt[i+1] & 0x 128 y = (mt[i] & 0x80000000) | (mt[i+1] & 0x7fffffff);
129 mt[i] = mt[i-NminusM] ^ (y >> 1) ^ ((y & 129 mt[i] = mt[i-NminusM] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
130 } 130 }
131 131
132 y = (mt[i] & 0x80000000) | (mt[0] & 0x7fff 132 y = (mt[i] & 0x80000000) | (mt[0] & 0x7fffffff);
133 mt[i] = mt[M-1] ^ (y >> 1) ^ ((y & 0x1) ? 133 mt[i] = mt[M-1] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
134 134
135 count624 = 0; 135 count624 = 0;
136 } 136 }
137 137
138 y = mt[count624]; 138 y = mt[count624];
139 y ^= ( y >> 11); 139 y ^= ( y >> 11);
140 y ^= ((y << 7 ) & 0x9d2c5680); 140 y ^= ((y << 7 ) & 0x9d2c5680);
141 y ^= ((y << 15) & 0xefc60000); 141 y ^= ((y << 15) & 0xefc60000);
142 y ^= ( y >> 18); 142 y ^= ( y >> 18);
143 143
144 return y * twoToMinus_3 144 return y * twoToMinus_32() + // Scale to range
145 (mt[count624++] >> 11) * twoToMinus_5 145 (mt[count624++] >> 11) * twoToMinus_53() + // fill remaining bits
146 nearlyTwoToMinus_54(); 146 nearlyTwoToMinus_54(); // make sure non-zero
147 } 147 }
148 148
149 void MTwistEngine::flatArray( const int size, 149 void MTwistEngine::flatArray( const int size, double *vect ) {
150 for( int i=0; i < size; ++i) vect[i] = flat( 150 for( int i=0; i < size; ++i) vect[i] = flat();
151 } 151 }
152 152
153 void MTwistEngine::setSeed(long seed, int k) { 153 void MTwistEngine::setSeed(long seed, int k) {
154 154
155 // MF 11/15/06 - Change seeding algorithm to 155 // MF 11/15/06 - Change seeding algorithm to a newer one recommended
156 // by Matsumoto: The original 156 // by Matsumoto: The original algorithm was
157 // mt[i] = (69069 * mt[i-1]) & 0xfffff 157 // mt[i] = (69069 * mt[i-1]) & 0xffffffff and this gives
158 // problems when the seed bit pattern 158 // problems when the seed bit pattern has lots of zeros
159 // (for example, 0x08000000). Savanah 159 // (for example, 0x08000000). Savanah bug #17479.
160 160
161 theSeed = seed ? seed : 4357; 161 theSeed = seed ? seed : 4357;
162 int mti; 162 int mti;
163 const int N1=624; 163 const int N1=624;
164 mt[0] = (unsigned int) (theSeed&0xffffffffUL 164 mt[0] = (unsigned int) (theSeed&0xffffffffUL);
165 for (mti=1; mti<N1; mti++) { 165 for (mti=1; mti<N1; mti++) {
166 mt[mti] = (1812433253UL * (mt[mti-1] ^ (mt 166 mt[mti] = (1812433253UL * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti);
167 /* See Knuth TAOCP Vol2. 3rd Ed. P.106 167 /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
168 /* In the previous versions, MSBs of t 168 /* In the previous versions, MSBs of the seed affect */
169 /* only MSBs of the array mt[]. 169 /* only MSBs of the array mt[]. */
170 /* 2002/01/09 modified by Makoto Matsu 170 /* 2002/01/09 modified by Makoto Matsumoto */
171 mt[mti] &= 0xffffffffUL; 171 mt[mti] &= 0xffffffffUL;
172 /* for >32 bit machines */ 172 /* for >32 bit machines */
173 } 173 }
174 for( int i=1; i < 624; ++i ) { 174 for( int i=1; i < 624; ++i ) {
175 mt[i] ^= k; // MF 9/16/98: distinguish 175 mt[i] ^= k; // MF 9/16/98: distinguish starting points
176 } 176 }
177 // MF 11/15/06 This distinction of starting 177 // MF 11/15/06 This distinction of starting points based on values of k
178 // is kept even though the seedi 178 // is kept even though the seeding algorithm itself is improved.
179 } 179 }
180 180
181 void MTwistEngine::setSeeds(const long *seeds, 181 void MTwistEngine::setSeeds(const long *seeds, int k) {
182 setSeed( (*seeds ? *seeds : 43571346), k ); 182 setSeed( (*seeds ? *seeds : 43571346), k );
183 int i; 183 int i;
184 for( i=1; i < 624; ++i ) { 184 for( i=1; i < 624; ++i ) {
185 mt[i] = ( seeds[1] + mt[i] ) & 0xffffffff; 185 mt[i] = ( seeds[1] + mt[i] ) & 0xffffffff; // MF 9/16/98
186 } 186 }
187 theSeeds = seeds; 187 theSeeds = seeds;
188 } 188 }
189 189
190 void MTwistEngine::saveStatus( const char file 190 void MTwistEngine::saveStatus( const char filename[] ) const
191 { 191 {
192 std::ofstream outFile( filename, std::ios:: 192 std::ofstream outFile( filename, std::ios::out ) ;
193 if (!outFile.bad()) { 193 if (!outFile.bad()) {
194 outFile << theSeed << std::endl; 194 outFile << theSeed << std::endl;
195 for (int i=0; i<624; ++i) outFile <<std:: 195 for (int i=0; i<624; ++i) outFile <<std::setprecision(20) << mt[i] << " ";
196 outFile << std::endl; 196 outFile << std::endl;
197 outFile << count624 << std::endl; 197 outFile << count624 << std::endl;
198 } 198 }
199 } 199 }
200 200
201 void MTwistEngine::restoreStatus( const char f 201 void MTwistEngine::restoreStatus( const char filename[] )
202 { 202 {
203 std::ifstream inFile( filename, std::ios::i 203 std::ifstream inFile( filename, std::ios::in);
204 if (!checkFile ( inFile, filename, engineNa 204 if (!checkFile ( inFile, filename, engineName(), "restoreStatus" )) {
205 std::cerr << " -- Engine state remains u 205 std::cerr << " -- Engine state remains unchanged\n";
206 return; 206 return;
207 } 207 }
208 208
209 if (!inFile.bad() && !inFile.eof()) { 209 if (!inFile.bad() && !inFile.eof()) {
210 inFile >> theSeed; 210 inFile >> theSeed;
211 for (int i=0; i<624; ++i) inFile >> mt[i] 211 for (int i=0; i<624; ++i) inFile >> mt[i];
212 inFile >> count624; 212 inFile >> count624;
213 } 213 }
214 } 214 }
215 215
216 void MTwistEngine::showStatus() const 216 void MTwistEngine::showStatus() const
217 { 217 {
218 std::cout << std::endl; 218 std::cout << std::endl;
219 std::cout << "--------- MTwist engine statu 219 std::cout << "--------- MTwist engine status ---------" << std::endl;
220 std::cout << std::setprecision(20); 220 std::cout << std::setprecision(20);
221 std::cout << " Initial seed = " << the 221 std::cout << " Initial seed = " << theSeed << std::endl;
222 std::cout << " Current index = " << cou 222 std::cout << " Current index = " << count624 << std::endl;
223 std::cout << " Array status mt[] = " << std 223 std::cout << " Array status mt[] = " << std::endl;
224 // 2014/06/06 L Garren 224 // 2014/06/06 L Garren
225 // the final line has 4 elements, not 5 225 // the final line has 4 elements, not 5
226 for (int i=0; i<620; i+=5) { 226 for (int i=0; i<620; i+=5) {
227 std::cout << mt[i] << " " << mt[i+1] << 227 std::cout << mt[i] << " " << mt[i+1] << " " << mt[i+2] << " "
228 << mt[i+3] << " " << mt[i+4] << "\n"; 228 << mt[i+3] << " " << mt[i+4] << "\n";
229 } 229 }
230 std::cout << mt[620] << " " << mt[621] << 230 std::cout << mt[620] << " " << mt[621] << " " << mt[622] << " "
231 << mt[623] << std::endl; 231 << mt[623] << std::endl;
232 std::cout << "----------------------------- 232 std::cout << "----------------------------------------" << std::endl;
233 } 233 }
234 234
235 MTwistEngine::operator double() { 235 MTwistEngine::operator double() {
236 return flat(); 236 return flat();
237 } 237 }
238 238
239 MTwistEngine::operator float() { 239 MTwistEngine::operator float() {
240 unsigned int y; 240 unsigned int y;
241 241
242 if( count624 >= N ) { 242 if( count624 >= N ) {
243 int i; 243 int i;
244 244
245 for( i=0; i < NminusM; ++i ) { 245 for( i=0; i < NminusM; ++i ) {
246 y = (mt[i] & 0x80000000) | (mt[i+1] & 0x 246 y = (mt[i] & 0x80000000) | (mt[i+1] & 0x7fffffff);
247 mt[i] = mt[i+M] ^ (y >> 1) ^ ((y & 247 mt[i] = mt[i+M] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
248 } 248 }
249 249
250 for( ; i < N-1 ; ++i ) { 250 for( ; i < N-1 ; ++i ) {
251 y = (mt[i] & 0x80000000) | (mt[i+1] & 0x 251 y = (mt[i] & 0x80000000) | (mt[i+1] & 0x7fffffff);
252 mt[i] = mt[i-NminusM] ^ (y >> 1) ^ ((y & 252 mt[i] = mt[i-NminusM] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
253 } 253 }
254 254
255 y = (mt[i] & 0x80000000) | (mt[0] & 0x7fff 255 y = (mt[i] & 0x80000000) | (mt[0] & 0x7fffffff);
256 mt[i] = mt[M-1] ^ (y >> 1) ^ ((y & 0x1) ? 256 mt[i] = mt[M-1] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
257 257
258 count624 = 0; 258 count624 = 0;
259 } 259 }
260 260
261 y = mt[count624++]; 261 y = mt[count624++];
262 y ^= ( y >> 11); 262 y ^= ( y >> 11);
263 y ^= ((y << 7 ) & 0x9d2c5680); 263 y ^= ((y << 7 ) & 0x9d2c5680);
264 y ^= ((y << 15) & 0xefc60000); 264 y ^= ((y << 15) & 0xefc60000);
265 y ^= ( y >> 18); 265 y ^= ( y >> 18);
266 266
267 return (float)(y * twoToMinus_32()); 267 return (float)(y * twoToMinus_32());
268 } 268 }
269 269
270 MTwistEngine::operator unsigned int() { 270 MTwistEngine::operator unsigned int() {
271 unsigned int y; 271 unsigned int y;
272 272
273 if( count624 >= N ) { 273 if( count624 >= N ) {
274 int i; 274 int i;
275 275
276 for( i=0; i < NminusM; ++i ) { 276 for( i=0; i < NminusM; ++i ) {
277 y = (mt[i] & 0x80000000) | (mt[i+1] & 0x 277 y = (mt[i] & 0x80000000) | (mt[i+1] & 0x7fffffff);
278 mt[i] = mt[i+M] ^ (y >> 1) ^ ((y & 278 mt[i] = mt[i+M] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
279 } 279 }
280 280
281 for( ; i < N-1 ; ++i ) { 281 for( ; i < N-1 ; ++i ) {
282 y = (mt[i] & 0x80000000) | (mt[i+1] & 0x 282 y = (mt[i] & 0x80000000) | (mt[i+1] & 0x7fffffff);
283 mt[i] = mt[i-NminusM] ^ (y >> 1) ^ ((y & 283 mt[i] = mt[i-NminusM] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
284 } 284 }
285 285
286 y = (mt[i] & 0x80000000) | (mt[0] & 0x7fff 286 y = (mt[i] & 0x80000000) | (mt[0] & 0x7fffffff);
287 mt[i] = mt[M-1] ^ (y >> 1) ^ ((y & 0x1) ? 287 mt[i] = mt[M-1] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
288 288
289 count624 = 0; 289 count624 = 0;
290 } 290 }
291 291
292 y = mt[count624++]; 292 y = mt[count624++];
293 y ^= ( y >> 11); 293 y ^= ( y >> 11);
294 y ^= ((y << 7 ) & 0x9d2c5680); 294 y ^= ((y << 7 ) & 0x9d2c5680);
295 y ^= ((y << 15) & 0xefc60000); 295 y ^= ((y << 15) & 0xefc60000);
296 y ^= ( y >> 18); 296 y ^= ( y >> 18);
297 297
298 return y; 298 return y;
299 } 299 }
300 300
301 std::ostream & MTwistEngine::put ( std::ostrea 301 std::ostream & MTwistEngine::put ( std::ostream& os ) const
302 { 302 {
303 char beginMarker[] = "MTwistEngine-begin"; 303 char beginMarker[] = "MTwistEngine-begin";
304 char endMarker[] = "MTwistEngine-end"; 304 char endMarker[] = "MTwistEngine-end";
305 305
306 long pr = os.precision(20); << 306 int pr = os.precision(20);
307 os << " " << beginMarker << " "; 307 os << " " << beginMarker << " ";
308 os << theSeed << " "; 308 os << theSeed << " ";
309 for (int i=0; i<624; ++i) { 309 for (int i=0; i<624; ++i) {
310 os << mt[i] << "\n"; 310 os << mt[i] << "\n";
311 } 311 }
312 os << count624 << " "; 312 os << count624 << " ";
313 os << endMarker << "\n"; 313 os << endMarker << "\n";
314 os.precision(pr); 314 os.precision(pr);
315 return os; 315 return os;
316 } 316 }
317 317
318 std::vector<unsigned long> MTwistEngine::put ( 318 std::vector<unsigned long> MTwistEngine::put () const {
319 std::vector<unsigned long> v; 319 std::vector<unsigned long> v;
320 v.push_back (engineIDulong<MTwistEngine>()); 320 v.push_back (engineIDulong<MTwistEngine>());
321 for (int i=0; i<624; ++i) { 321 for (int i=0; i<624; ++i) {
322 v.push_back(static_cast<unsigned long>(mt 322 v.push_back(static_cast<unsigned long>(mt[i]));
323 } 323 }
324 v.push_back(count624); 324 v.push_back(count624);
325 return v; 325 return v;
326 } 326 }
327 327
328 std::istream & MTwistEngine::get ( std::istre 328 std::istream & MTwistEngine::get ( std::istream& is )
329 { 329 {
330 char beginMarker [MarkerLen]; 330 char beginMarker [MarkerLen];
331 is >> std::ws; 331 is >> std::ws;
332 is.width(MarkerLen); // causes the next rea 332 is.width(MarkerLen); // causes the next read to the char* to be <=
333 // that many bytes, INCLUDING A TERMINAT 333 // that many bytes, INCLUDING A TERMINATION \0
334 // (Stroustrup, section 21.3.2) 334 // (Stroustrup, section 21.3.2)
335 is >> beginMarker; 335 is >> beginMarker;
336 if (strcmp(beginMarker,"MTwistEngine-begin") 336 if (strcmp(beginMarker,"MTwistEngine-begin")) {
337 is.clear(std::ios::badbit | is.rdstate()) 337 is.clear(std::ios::badbit | is.rdstate());
338 std::cerr << "\nInput stream mispositione 338 std::cerr << "\nInput stream mispositioned or"
339 << "\nMTwistEngine state description 339 << "\nMTwistEngine state description missing or"
340 << "\nwrong engine type found." << st 340 << "\nwrong engine type found." << std::endl;
341 return is; 341 return is;
342 } 342 }
343 return getState(is); 343 return getState(is);
344 } 344 }
345 345
346 std::string MTwistEngine::beginTag ( ) { 346 std::string MTwistEngine::beginTag ( ) {
347 return "MTwistEngine-begin"; 347 return "MTwistEngine-begin";
348 } 348 }
349 349
350 std::istream & MTwistEngine::getState ( std:: 350 std::istream & MTwistEngine::getState ( std::istream& is )
351 { 351 {
352 char endMarker [MarkerLen]; 352 char endMarker [MarkerLen];
353 is >> theSeed; 353 is >> theSeed;
354 for (int i=0; i<624; ++i) is >> mt[i]; 354 for (int i=0; i<624; ++i) is >> mt[i];
355 is >> count624; 355 is >> count624;
356 is >> std::ws; 356 is >> std::ws;
357 is.width(MarkerLen); 357 is.width(MarkerLen);
358 is >> endMarker; 358 is >> endMarker;
359 if (strcmp(endMarker,"MTwistEngine-end")) { 359 if (strcmp(endMarker,"MTwistEngine-end")) {
360 is.clear(std::ios::badbit | is.rdstate()) 360 is.clear(std::ios::badbit | is.rdstate());
361 std::cerr << "\nMTwistEngine state descri 361 std::cerr << "\nMTwistEngine state description incomplete."
362 << "\nInput stream is probably mispos 362 << "\nInput stream is probably mispositioned now." << std::endl;
363 return is; 363 return is;
364 } 364 }
365 return is; 365 return is;
366 } 366 }
367 367
368 bool MTwistEngine::get (const std::vector<unsi 368 bool MTwistEngine::get (const std::vector<unsigned long> & v) {
369 if ((v[0] & 0xffffffffUL) != engineIDulong<M 369 if ((v[0] & 0xffffffffUL) != engineIDulong<MTwistEngine>()) {
370 std::cerr << 370 std::cerr <<
371 "\nMTwistEngine get:state vector has wro 371 "\nMTwistEngine get:state vector has wrong ID word - state unchanged\n";
372 return false; 372 return false;
373 } 373 }
374 return getState(v); 374 return getState(v);
375 } 375 }
376 376
377 bool MTwistEngine::getState (const std::vector 377 bool MTwistEngine::getState (const std::vector<unsigned long> & v) {
378 if (v.size() != VECTOR_STATE_SIZE ) { 378 if (v.size() != VECTOR_STATE_SIZE ) {
379 std::cerr << 379 std::cerr <<
380 "\nMTwistEngine get:state vector has wro 380 "\nMTwistEngine get:state vector has wrong length - state unchanged\n";
381 return false; 381 return false;
382 } 382 }
383 for (int i=0; i<624; ++i) { 383 for (int i=0; i<624; ++i) {
384 mt[i]=(unsigned int)v[i+1]; << 384 mt[i]=v[i+1];
385 } 385 }
386 count624 = (int)v[625]; << 386 count624 = v[625];
387 return true; 387 return true;
388 } 388 }
389 389
390 } // namespace CLHEP 390 } // namespace CLHEP
391 391