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

MTwistEngine.cc
Go to the documentation of this file.
1 // $Id: MTwistEngine.cc,v 1.6 2010/06/16 17:24:53 garren Exp $
2 // -*- C++ -*-
3 //
4 // -----------------------------------------------------------------------
5 // HEP Random
6 // --- MTwistEngine ---
7 // class implementation file
8 // -----------------------------------------------------------------------
9 // A "fast, compact, huge-period generator" based on M. Matsumoto and
10 // T. Nishimura, "Mersenne Twister: A 623-dimensionally equidistributed
11 // uniform pseudorandom number generator", to appear in ACM Trans. on
12 // Modeling and Computer Simulation. It is a twisted GFSR generator
13 // with a Mersenne-prime period of 2^19937-1, uniform on open interval (0,1)
14 // =======================================================================
15 // Ken Smith - Started initial draft: 14th Jul 1998
16 // - Optimized to get std::pow() out of flat() method
17 // - Added conversion operators: 6th Aug 1998
18 // J. Marraffino - Added some explicit casts to deal with
19 // machines where sizeof(int) != sizeof(long) 22 Aug 1998
20 // M. Fischler - Modified constructors such that no two
21 // seeds will match sequences, no single/double
22 // seeds will match, explicit seeds give
23 // determined results, and default will not
24 // match any of the above or other defaults.
25 // - Modified use of the various exponents of 2
26 // to avoid per-instance space overhead and
27 // correct the rounding procedure 16 Sep 1998
28 // J. Marfaffino - Remove dependence on hepString class 13 May 1999
29 // M. Fischler - In restore, checkFile for file not found 03 Dec 2004
30 // M. Fischler - Methods for distrib. instacne save/restore 12/8/04
31 // M. Fischler - split get() into tag validation and
32 // getState() for anonymous restores 12/27/04
33 // M. Fischler - put/get for vectors of ulongs 3/14/05
34 // M. Fischler - State-saving using only ints, for portability 4/12/05
35 // M. Fischler - Improved seeding in setSeed (Savanah bug #17479) 11/15/06
36 // - (Possible optimization - now that the seeding is improved,
37 // is it still necessary for ctor to "warm up" by discarding
38 // 2000 iterations?)
39 //
40 // =======================================================================
41 
42 #include "CLHEP/Random/defs.h"
43 #include "CLHEP/Random/Random.h"
44 #include "CLHEP/Random/MTwistEngine.h"
45 #include "CLHEP/Random/engineIDulong.h"
46 #include <string.h> // for strcmp
47 #include <cstdlib> // for std::abs(int)
48 
49 namespace CLHEP {
50 
51 static const int MarkerLen = 64; // Enough room to hold a begin or end marker.
52 
53 std::string MTwistEngine::name() const {return "MTwistEngine";}
54 
55 int MTwistEngine::numEngines = 0;
56 int MTwistEngine::maxIndex = 215;
57 
60 {
61  int cycle = std::abs(int(numEngines/maxIndex));
62  int curIndex = std::abs(int(numEngines%maxIndex));
63  long mask = ((cycle & 0x007fffff) << 8);
64  long seedlist[2];
65  HepRandom::getTheTableSeeds( seedlist, curIndex );
66  seedlist[0] = (seedlist[0])^mask;
67  seedlist[1] = 0;
68  setSeeds( seedlist, numEngines );
69  count624=0;
70  ++numEngines;
71  for( int i=0; i < 2000; ++i ) flat(); // Warm up just a bit
72 }
73 
76 {
77  long seedlist[2]={seed,17587};
78  setSeeds( seedlist, 0 );
79  count624=0;
80  for( int i=0; i < 2000; ++i ) flat(); // Warm up just a bit
81 }
82 
83 MTwistEngine::MTwistEngine(int rowIndex, int colIndex)
85 {
86  int cycle = std::abs(int(rowIndex/maxIndex));
87  int row = std::abs(int(rowIndex%maxIndex));
88  int col = std::abs(int(colIndex%2));
89  long mask = (( cycle & 0x000007ff ) << 20 );
90  long seedlist[2];
91  HepRandom::getTheTableSeeds( seedlist, row );
92  seedlist[0] = (seedlist[col])^mask;
93  seedlist[1] = 690691;
94  setSeeds(seedlist, 4444772);
95  count624=0;
96  for( int i=0; i < 2000; ++i ) flat(); // Warm up just a bit
97 }
98 
99 MTwistEngine::MTwistEngine( std::istream& is )
100 : HepRandomEngine()
101 {
102  is >> *this;
103 }
104 
106 
108  unsigned int y;
109 
110  if( count624 >= N ) {
111  register int i;
112 
113  for( i=0; i < NminusM; ++i ) {
114  y = (mt[i] & 0x80000000) | (mt[i+1] & 0x7fffffff);
115  mt[i] = mt[i+M] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
116  }
117 
118  for( ; i < N-1 ; ++i ) {
119  y = (mt[i] & 0x80000000) | (mt[i+1] & 0x7fffffff);
120  mt[i] = mt[i-NminusM] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
121  }
122 
123  y = (mt[i] & 0x80000000) | (mt[0] & 0x7fffffff);
124  mt[i] = mt[M-1] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
125 
126  count624 = 0;
127  }
128 
129  y = mt[count624];
130  y ^= ( y >> 11);
131  y ^= ((y << 7 ) & 0x9d2c5680);
132  y ^= ((y << 15) & 0xefc60000);
133  y ^= ( y >> 18);
134 
135  return y * twoToMinus_32() + // Scale to range
136  (mt[count624++] >> 11) * twoToMinus_53() + // fill remaining bits
137  nearlyTwoToMinus_54(); // make sure non-zero
138 }
139 
140 void MTwistEngine::flatArray( const int size, double *vect ) {
141  for( int i=0; i < size; ++i) vect[i] = flat();
142 }
143 
144 void MTwistEngine::setSeed(long seed, int k) {
145 
146  // MF 11/15/06 - Change seeding algorithm to a newer one recommended
147  // by Matsumoto: The original algorithm was
148  // mt[i] = (69069 * mt[i-1]) & 0xffffffff and this gives
149  // problems when the seed bit pattern has lots of zeros
150  // (for example, 0x08000000). Savanah bug #17479.
151 
152  theSeed = seed ? seed : 4357;
153  int mti;
154  const int N1=624;
155  mt[0] = (unsigned int) (theSeed&0xffffffffUL);
156  for (mti=1; mti<N1; mti++) {
157  mt[mti] = (1812433253UL * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti);
158  /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
159  /* In the previous versions, MSBs of the seed affect */
160  /* only MSBs of the array mt[]. */
161  /* 2002/01/09 modified by Makoto Matsumoto */
162  mt[mti] &= 0xffffffffUL;
163  /* for >32 bit machines */
164  }
165  for( int i=1; i < 624; ++i ) {
166  mt[i] ^= k; // MF 9/16/98: distinguish starting points
167  }
168  // MF 11/15/06 This distinction of starting points based on values of k
169  // is kept even though the seeding algorithm itself is improved.
170 }
171 
172 void MTwistEngine::setSeeds(const long *seeds, int k) {
173  setSeed( (*seeds ? *seeds : 43571346), k );
174  int i;
175  for( i=1; i < 624; ++i ) {
176  mt[i] = ( seeds[1] + mt[i] ) & 0xffffffff; // MF 9/16/98
177  }
178  theSeeds = seeds;
179 }
180 
181 void MTwistEngine::saveStatus( const char filename[] ) const
182 {
183  std::ofstream outFile( filename, std::ios::out ) ;
184  if (!outFile.bad()) {
185  outFile << theSeed << std::endl;
186  for (int i=0; i<624; ++i) outFile <<std::setprecision(20) << mt[i] << " ";
187  outFile << std::endl;
188  outFile << count624 << std::endl;
189  }
190 }
191 
192 void MTwistEngine::restoreStatus( const char filename[] )
193 {
194  std::ifstream inFile( filename, std::ios::in);
195  if (!checkFile ( inFile, filename, engineName(), "restoreStatus" )) {
196  std::cerr << " -- Engine state remains unchanged\n";
197  return;
198  }
199 
200  if (!inFile.bad() && !inFile.eof()) {
201  inFile >> theSeed;
202  for (int i=0; i<624; ++i) inFile >> mt[i];
203  inFile >> count624;
204  }
205 }
206 
208 {
209  std::cout << std::endl;
210  std::cout << "--------- MTwist engine status ---------" << std::endl;
211  std::cout << std::setprecision(20);
212  std::cout << " Initial seed = " << theSeed << std::endl;
213  std::cout << " Current index = " << count624 << std::endl;
214  std::cout << " Array status mt[] = " << std::endl;
215  for (int i=0; i<624; i+=5) {
216  std::cout << mt[i] << " " << mt[i+1] << " " << mt[i+2] << " "
217  << mt[i+3] << " " << mt[i+4] << std::endl;
218  }
219  std::cout << "----------------------------------------" << std::endl;
220 }
221 
222 MTwistEngine::operator float() {
223  unsigned int y;
224 
225  if( count624 >= N ) {
226  register int i;
227 
228  for( i=0; i < NminusM; ++i ) {
229  y = (mt[i] & 0x80000000) | (mt[i+1] & 0x7fffffff);
230  mt[i] = mt[i+M] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
231  }
232 
233  for( ; i < N-1 ; ++i ) {
234  y = (mt[i] & 0x80000000) | (mt[i+1] & 0x7fffffff);
235  mt[i] = mt[i-NminusM] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
236  }
237 
238  y = (mt[i] & 0x80000000) | (mt[0] & 0x7fffffff);
239  mt[i] = mt[M-1] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
240 
241  count624 = 0;
242  }
243 
244  y = mt[count624++];
245  y ^= ( y >> 11);
246  y ^= ((y << 7 ) & 0x9d2c5680);
247  y ^= ((y << 15) & 0xefc60000);
248  y ^= ( y >> 18);
249 
250  return (float)(y * twoToMinus_32());
251 }
252 
253 MTwistEngine::operator unsigned int() {
254  unsigned int y;
255 
256  if( count624 >= N ) {
257  register int i;
258 
259  for( i=0; i < NminusM; ++i ) {
260  y = (mt[i] & 0x80000000) | (mt[i+1] & 0x7fffffff);
261  mt[i] = mt[i+M] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
262  }
263 
264  for( ; i < N-1 ; ++i ) {
265  y = (mt[i] & 0x80000000) | (mt[i+1] & 0x7fffffff);
266  mt[i] = mt[i-NminusM] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
267  }
268 
269  y = (mt[i] & 0x80000000) | (mt[0] & 0x7fffffff);
270  mt[i] = mt[M-1] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
271 
272  count624 = 0;
273  }
274 
275  y = mt[count624++];
276  y ^= ( y >> 11);
277  y ^= ((y << 7 ) & 0x9d2c5680);
278  y ^= ((y << 15) & 0xefc60000);
279  y ^= ( y >> 18);
280 
281  return y;
282 }
283 
284 std::ostream & MTwistEngine::put ( std::ostream& os ) const
285 {
286  char beginMarker[] = "MTwistEngine-begin";
287  char endMarker[] = "MTwistEngine-end";
288 
289  int pr = os.precision(20);
290  os << " " << beginMarker << " ";
291  os << theSeed << " ";
292  for (int i=0; i<624; ++i) {
293  os << mt[i] << "\n";
294  }
295  os << count624 << " ";
296  os << endMarker << "\n";
297  os.precision(pr);
298  return os;
299 }
300 
301 std::vector<unsigned long> MTwistEngine::put () const {
302  std::vector<unsigned long> v;
303  v.push_back (engineIDulong<MTwistEngine>());
304  for (int i=0; i<624; ++i) {
305  v.push_back(static_cast<unsigned long>(mt[i]));
306  }
307  v.push_back(count624);
308  return v;
309 }
310 
311 std::istream & MTwistEngine::get ( std::istream& is )
312 {
313  char beginMarker [MarkerLen];
314  is >> std::ws;
315  is.width(MarkerLen); // causes the next read to the char* to be <=
316  // that many bytes, INCLUDING A TERMINATION \0
317  // (Stroustrup, section 21.3.2)
318  is >> beginMarker;
319  if (strcmp(beginMarker,"MTwistEngine-begin")) {
320  is.clear(std::ios::badbit | is.rdstate());
321  std::cerr << "\nInput stream mispositioned or"
322  << "\nMTwistEngine state description missing or"
323  << "\nwrong engine type found." << std::endl;
324  return is;
325  }
326  return getState(is);
327 }
328 
329 std::string MTwistEngine::beginTag ( ) {
330  return "MTwistEngine-begin";
331 }
332 
333 std::istream & MTwistEngine::getState ( std::istream& is )
334 {
335  char endMarker [MarkerLen];
336  is >> theSeed;
337  for (int i=0; i<624; ++i) is >> mt[i];
338  is >> count624;
339  is >> std::ws;
340  is.width(MarkerLen);
341  is >> endMarker;
342  if (strcmp(endMarker,"MTwistEngine-end")) {
343  is.clear(std::ios::badbit | is.rdstate());
344  std::cerr << "\nMTwistEngine state description incomplete."
345  << "\nInput stream is probably mispositioned now." << std::endl;
346  return is;
347  }
348  return is;
349 }
350 
351 bool MTwistEngine::get (const std::vector<unsigned long> & v) {
352  if ((v[0] & 0xffffffffUL) != engineIDulong<MTwistEngine>()) {
353  std::cerr <<
354  "\nMTwistEngine get:state vector has wrong ID word - state unchanged\n";
355  return false;
356  }
357  return getState(v);
358 }
359 
360 bool MTwistEngine::getState (const std::vector<unsigned long> & v) {
361  if (v.size() != VECTOR_STATE_SIZE ) {
362  std::cerr <<
363  "\nMTwistEngine get:state vector has wrong length - state unchanged\n";
364  return false;
365  }
366  for (int i=0; i<624; ++i) {
367  mt[i]=v[i+1];
368  }
369  count624 = v[625];
370  return true;
371 }
372 
373 } // namespace CLHEP
CLHEP::HepRandomEngine
Definition: Matrix/CLHEP/Random/RandomEngine.h:55
CLHEP::MTwistEngine::restoreStatus
void restoreStatus(const char filename[]="MTwist.conf")
Definition: MTwistEngine.cc:192
CLHEP::HepRandomEngine::theSeed
long theSeed
Definition: Matrix/CLHEP/Random/RandomEngine.h:144
CLHEP::MTwistEngine::~MTwistEngine
virtual ~MTwistEngine()
Definition: MTwistEngine.cc:105
CLHEP::MTwistEngine::engineName
static std::string engineName()
Definition: Matrix/CLHEP/Random/MTwistEngine.h:78
CLHEP::HepRandom::getTheTableSeeds
static void getTheTableSeeds(long *seeds, int index)
Definition: Random.cc:152
is
HepRotation and so forth isNear() norm2() rectify() static Rotation row1 row4(To avoid bloat in the code pulled in for programs which don 't use all these features, we split the implementation .cc files. Only isNear() goes into the original Rotation.cc) --------------------------------------- HepAxisAngle and HepEulerAngles classes --------------------------------------- These classes are very useful and simple structures for holding the result of a nice intuituve decomposition of a rotation there is no longer much content in the distinct ZOOM PhysicsVectors library The only content left in the library is the object files representing the various Exception objects When we build the CLHEP classes for the ZOOM we will set up so as to use ZOOM SpaceVector is(but we can disable namespace usage and most of our users do so at this point). What I do is leave Hep3Vector in the global namespace
CLHEP::MTwistEngine::flat
double flat()
Definition: MTwistEngine.cc:107
CLHEP::HepRandomEngine::theSeeds
const long * theSeeds
Definition: Matrix/CLHEP/Random/RandomEngine.h:145
CLHEP::HepRandomEngine::nearlyTwoToMinus_54
static double nearlyTwoToMinus_54()
size
user code seldom needs to call this function directly ZMerrno whether or not they are still recorded ZMerrno size() Return the(integer) number of ZMthrow 'n exceptions currently recorded. 5) ZMerrno.clear() Set an internal counter to zero. This counter is available(see next function) to user code to track ZMthrow 'n exceptions that have occurred during any arbitrary time interval. 6) ZMerrno.countSinceCleared() Return the(integer) number of ZMthrow 'n exceptions that have been recorded via ZMerrno.write()
N
the goal is to keep the overall false rejection probability down at the to level For each validated we discuss which of course is by necessity relative timing We take the time for a single random via one of the fastest good and at any rate the ratios will vary by around depending on the processor and memory configuration used A timing for a distribution of units would mean no time used beyond the uniform random Summary Distribution Validated Validation Rejected Past N RandGauss N
Definition: validation.doc:48
CLHEP
Definition: ClhepVersion.h:13
CLHEP::HepRandomEngine::twoToMinus_53
static double twoToMinus_53()
v
they are gone ZOOM Features Discontinued The following features of the ZOOM package were felt to be extreme overkill These have been after checking that no existing user code was utilizing as in SpaceVector v
Definition: keyMergeIssues.doc:324
CLHEP::MTwistEngine::setSeeds
void setSeeds(const long *seeds, int)
Definition: MTwistEngine.cc:172
CLHEP::MTwistEngine::get
virtual std::istream & get(std::istream &is)
Definition: MTwistEngine.cc:311
CLHEP::MTwistEngine::MTwistEngine
MTwistEngine()
Definition: MTwistEngine.cc:58
CLHEP::MTwistEngine::VECTOR_STATE_SIZE
static const unsigned int VECTOR_STATE_SIZE
Definition: Matrix/CLHEP/Random/MTwistEngine.h:84
seeds
Technical Maintenance Note for CLHEP Random Consequences of seeding JamesRandom with positive seed values greater than In the source code JamesRandom The usual way of seeding a generator is via the default which makes use of the table of seeds(with some trickery to ensure that the values won 't repeat after the table rows are exhausted). The trickery preserves the fact that sees are never negative(because the table values are never negative
CLHEP::MTwistEngine::beginTag
static std::string beginTag()
Definition: MTwistEngine.cc:329
CLHEP::MTwistEngine::getState
virtual std::istream & getState(std::istream &is)
Definition: MTwistEngine.cc:333
CLHEP::MTwistEngine::setSeed
void setSeed(long seed, int)
Definition: MTwistEngine.cc:144
CLHEP::MTwistEngine::put
std::vector< unsigned long > put() const
Definition: MTwistEngine.cc:301
CLHEP::HepRandomEngine::twoToMinus_32
static double twoToMinus_32()
i
long i
Definition: JamesRandomSeeding.txt:27
CLHEP::MTwistEngine::flatArray
void flatArray(const int size, double *vect)
Definition: MTwistEngine.cc:140
CLHEP::MTwistEngine::showStatus
void showStatus() const
Definition: MTwistEngine.cc:207
CLHEP::HepRandomEngine::checkFile
static bool checkFile(std::istream &file, const std::string &filename, const std::string &classname, const std::string &methodname)
Definition: RandomEngine.cc:46
in
it has advantages For I leave the ZMthrows in
Definition: keyMergeIssues.doc:62
k
long k
Definition: JamesRandomSeeding.txt:29
CLHEP::MTwistEngine::name
std::string name() const
Definition: MTwistEngine.cc:53
CLHEP::MTwistEngine::saveStatus
void saveStatus(const char filename[]="MTwist.conf") const
Definition: MTwistEngine.cc:181