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

MTwistEngine.cc
Go to the documentation of this file.
00001 // $Id: MTwistEngine.cc,v 1.6 2010/06/16 17:24:53 garren Exp $
00002 // -*- C++ -*-
00003 //
00004 // -----------------------------------------------------------------------
00005 //                             HEP Random
00006 //                        --- MTwistEngine ---
00007 //                      class implementation file
00008 // -----------------------------------------------------------------------
00009 // A "fast, compact, huge-period generator" based on M. Matsumoto and 
00010 // T. Nishimura, "Mersenne Twister: A 623-dimensionally equidistributed 
00011 // uniform pseudorandom number generator", to appear in ACM Trans. on
00012 // Modeling and Computer Simulation.  It is a twisted GFSR generator
00013 // with a Mersenne-prime period of 2^19937-1, uniform on open interval (0,1)
00014 // =======================================================================
00015 // Ken Smith      - Started initial draft:                    14th Jul 1998
00016 //                - Optimized to get pow() out of flat() method
00017 //                - Added conversion operators:                6th Aug 1998
00018 // J. Marraffino  - Added some explicit casts to deal with
00019 //                  machines where sizeof(int) != sizeof(long)  22 Aug 1998
00020 // M. Fischler    - Modified constructors such that no two
00021 //                  seeds will match sequences, no single/double
00022 //                  seeds will match, explicit seeds give 
00023 //                  determined results, and default will not 
00024 //                  match any of the above or other defaults.
00025 //                - Modified use of the various exponents of 2
00026 //                  to avoid per-instance space overhead and
00027 //                  correct the rounding procedure              16 Sep 1998
00028 // J. Marfaffino  - Remove dependence on hepString class        13 May 1999
00029 // M. Fischler    - In restore, checkFile for file not found    03 Dec 2004
00030 // M. Fischler    - Methods for distrib. instacne save/restore  12/8/04    
00031 // M. Fischler    - split get() into tag validation and 
00032 //                  getState() for anonymous restores           12/27/04    
00033 // M. Fischler    - put/get for vectors of ulongs               3/14/05
00034 // M. Fischler    - State-saving using only ints, for portability 4/12/05
00035 // M. Fischler    - Improved seeding in setSeed  (Savanah bug #17479) 11/15/06
00036 //                - (Possible optimization - now that the seeding is improved,
00037 //                  is it still necessary for ctor to "warm up" by discarding
00038 //                  2000 iterations?)
00039 //                  
00040 // =======================================================================
00041 
00042 #include "CLHEP/Random/defs.h"
00043 #include "CLHEP/Random/Random.h"
00044 #include "CLHEP/Random/MTwistEngine.h"
00045 #include "CLHEP/Random/engineIDulong.h"
00046 #include <string.h>     // for strcmp
00047 #include <cstdlib>      // for abs(int)
00048 
00049 using namespace std;
00050 
00051 namespace CLHEP {
00052 
00053 static const int MarkerLen = 64; // Enough room to hold a begin or end marker. 
00054 
00055 std::string MTwistEngine::name() const {return "MTwistEngine";}
00056 
00057 int MTwistEngine::numEngines = 0;
00058 int MTwistEngine::maxIndex = 215;
00059 
00060 MTwistEngine::MTwistEngine() 
00061 : HepRandomEngine()
00062 {
00063   int cycle = abs(int(numEngines/maxIndex));
00064   int curIndex = abs(int(numEngines%maxIndex));
00065   long mask = ((cycle & 0x007fffff) << 8);
00066   long seedlist[2];
00067   HepRandom::getTheTableSeeds( seedlist, curIndex );
00068   seedlist[0] = (seedlist[0])^mask;
00069   seedlist[1] = 0;
00070   setSeeds( seedlist, numEngines );
00071   count624=0;
00072   ++numEngines;
00073   for( int i=0; i < 2000; ++i ) flat();      // Warm up just a bit
00074 }
00075 
00076 MTwistEngine::MTwistEngine(long seed)  
00077 : HepRandomEngine()
00078 {
00079   long seedlist[2]={seed,17587};
00080   setSeeds( seedlist, 0 );
00081   count624=0;
00082   for( int i=0; i < 2000; ++i ) flat();      // Warm up just a bit
00083 }
00084 
00085 MTwistEngine::MTwistEngine(int rowIndex, int colIndex) 
00086 : HepRandomEngine()
00087 {
00088   int cycle = abs(int(rowIndex/maxIndex));
00089   int row = abs(int(rowIndex%maxIndex));
00090   int col = abs(int(colIndex%2));
00091   long mask = (( cycle & 0x000007ff ) << 20 );
00092   long seedlist[2];
00093   HepRandom::getTheTableSeeds( seedlist, row );
00094   seedlist[0] = (seedlist[col])^mask;
00095   seedlist[1] = 690691;
00096   setSeeds(seedlist, 4444772);
00097   count624=0;
00098   for( int i=0; i < 2000; ++i ) flat();      // Warm up just a bit
00099 }
00100 
00101 MTwistEngine::MTwistEngine( std::istream& is )  
00102 : HepRandomEngine()
00103 {
00104   is >> *this;
00105 }
00106 
00107 MTwistEngine::~MTwistEngine() {}
00108 
00109 double MTwistEngine::flat() {
00110   unsigned int y;
00111 
00112    if( count624 >= N ) {
00113     register int i;
00114 
00115     for( i=0; i < NminusM; ++i ) {
00116       y = (mt[i] & 0x80000000) | (mt[i+1] & 0x7fffffff);
00117       mt[i] = mt[i+M]       ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
00118     }
00119 
00120     for(    ; i < N-1    ; ++i ) {
00121       y = (mt[i] & 0x80000000) | (mt[i+1] & 0x7fffffff);
00122       mt[i] = mt[i-NminusM] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
00123     }
00124 
00125     y = (mt[i] & 0x80000000) | (mt[0] & 0x7fffffff);
00126     mt[i] = mt[M-1] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
00127 
00128     count624 = 0;
00129   }
00130 
00131   y = mt[count624];
00132   y ^= ( y >> 11);
00133   y ^= ((y << 7 ) & 0x9d2c5680);
00134   y ^= ((y << 15) & 0xefc60000);
00135   y ^= ( y >> 18);
00136 
00137   return                      y * twoToMinus_32()  +    // Scale to range 
00138          (mt[count624++] >> 11) * twoToMinus_53()  +    // fill remaining bits
00139                             nearlyTwoToMinus_54();      // make sure non-zero
00140 }
00141 
00142 void MTwistEngine::flatArray( const int size, double *vect ) {
00143   for( int i=0; i < size; ++i) vect[i] = flat();
00144 }
00145 
00146 void MTwistEngine::setSeed(long seed, int k) {
00147 
00148   // MF 11/15/06 - Change seeding algorithm to a newer one recommended 
00149   //               by Matsumoto: The original algorithm was 
00150   //               mt[i] = (69069 * mt[i-1]) & 0xffffffff and this gives
00151   //               problems when the seed bit pattern has lots of zeros
00152   //               (for example, 0x08000000).  Savanah bug #17479.
00153 
00154   theSeed = seed ? seed : 4357;
00155   int mti;
00156   const int N1=624;
00157   mt[0] = (unsigned int) (theSeed&0xffffffffUL);
00158   for (mti=1; mti<N1; mti++) {
00159     mt[mti] = (1812433253UL * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti); 
00160         /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
00161         /* In the previous versions, MSBs of the seed affect   */
00162         /* only MSBs of the array mt[].                        */
00163         /* 2002/01/09 modified by Makoto Matsumoto             */
00164     mt[mti] &= 0xffffffffUL;
00165         /* for >32 bit machines */
00166   }
00167   for( int i=1; i < 624; ++i ) {
00168     mt[i] ^= k;                 // MF 9/16/98: distinguish starting points
00169   }
00170   // MF 11/15/06 This distinction of starting points based on values of k
00171   //             is kept even though the seeding algorithm itself is improved.
00172 }
00173 
00174 void MTwistEngine::setSeeds(const long *seeds, int k) {
00175   setSeed( (*seeds ? *seeds : 43571346), k );
00176   int i;
00177   for( i=1; i < 624; ++i ) {
00178     mt[i] = ( seeds[1] + mt[i] ) & 0xffffffff; // MF 9/16/98 
00179   }
00180   theSeeds = seeds;
00181 }
00182 
00183 void MTwistEngine::saveStatus( const char filename[] ) const
00184 {
00185    std::ofstream outFile( filename, std::ios::out ) ;
00186    if (!outFile.bad()) {
00187      outFile << theSeed << std::endl;
00188      for (int i=0; i<624; ++i) outFile <<std::setprecision(20) << mt[i] << " ";
00189      outFile << std::endl;
00190      outFile << count624 << std::endl;
00191    }
00192 }
00193 
00194 void MTwistEngine::restoreStatus( const char filename[] )
00195 {
00196    std::ifstream inFile( filename, std::ios::in);
00197    if (!checkFile ( inFile, filename, engineName(), "restoreStatus" )) {
00198      std::cerr << "  -- Engine state remains unchanged\n";
00199      return;
00200    }
00201 
00202    if (!inFile.bad() && !inFile.eof()) {
00203      inFile >> theSeed;
00204      for (int i=0; i<624; ++i) inFile >> mt[i];
00205      inFile >> count624;
00206    }
00207 }
00208 
00209 void MTwistEngine::showStatus() const
00210 {
00211    std::cout << std::endl;
00212    std::cout << "--------- MTwist engine status ---------" << std::endl;
00213    std::cout << std::setprecision(20);
00214    std::cout << " Initial seed      = " << theSeed << std::endl;
00215    std::cout << " Current index     = " << count624 << std::endl;
00216    std::cout << " Array status mt[] = " << std::endl;
00217    for (int i=0; i<624; i+=5) {
00218      std::cout << mt[i]   << " " << mt[i+1] << " " << mt[i+2] << " " 
00219                << mt[i+3] << " " << mt[i+4] << std::endl;
00220    }
00221    std::cout << "----------------------------------------" << std::endl;
00222 }
00223 
00224 MTwistEngine::operator float() {
00225   unsigned int y;
00226 
00227   if( count624 >= N ) {
00228     register int i;
00229 
00230     for( i=0; i < NminusM; ++i ) {
00231       y = (mt[i] & 0x80000000) | (mt[i+1] & 0x7fffffff);
00232       mt[i] = mt[i+M]       ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
00233     }
00234 
00235     for(    ; i < N-1    ; ++i ) {
00236       y = (mt[i] & 0x80000000) | (mt[i+1] & 0x7fffffff);
00237       mt[i] = mt[i-NminusM] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
00238     }
00239 
00240     y = (mt[i] & 0x80000000) | (mt[0] & 0x7fffffff);
00241     mt[i] = mt[M-1] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
00242 
00243     count624 = 0;
00244   }
00245 
00246   y = mt[count624++];
00247   y ^= ( y >> 11);
00248   y ^= ((y << 7 ) & 0x9d2c5680);
00249   y ^= ((y << 15) & 0xefc60000);
00250   y ^= ( y >> 18);
00251 
00252   return (float)(y * twoToMinus_32());
00253 }
00254 
00255 MTwistEngine::operator unsigned int() {
00256   unsigned int y;
00257 
00258   if( count624 >= N ) {
00259     register int i;
00260 
00261     for( i=0; i < NminusM; ++i ) {
00262       y = (mt[i] & 0x80000000) | (mt[i+1] & 0x7fffffff);
00263       mt[i] = mt[i+M]       ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
00264     }
00265 
00266     for(    ; i < N-1    ; ++i ) {
00267       y = (mt[i] & 0x80000000) | (mt[i+1] & 0x7fffffff);
00268       mt[i] = mt[i-NminusM] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
00269     }
00270 
00271     y = (mt[i] & 0x80000000) | (mt[0] & 0x7fffffff);
00272     mt[i] = mt[M-1] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
00273 
00274     count624 = 0;
00275   }
00276 
00277   y = mt[count624++];
00278   y ^= ( y >> 11);
00279   y ^= ((y << 7 ) & 0x9d2c5680);
00280   y ^= ((y << 15) & 0xefc60000);
00281   y ^= ( y >> 18);
00282 
00283   return y;
00284 }
00285 
00286 std::ostream & MTwistEngine::put ( std::ostream& os ) const
00287 {
00288    char beginMarker[] = "MTwistEngine-begin";
00289    char endMarker[]   = "MTwistEngine-end";
00290 
00291    int pr = os.precision(20);
00292    os << " " << beginMarker << " ";
00293    os << theSeed << " ";
00294    for (int i=0; i<624; ++i) {
00295      os << mt[i] << "\n";
00296    }
00297    os << count624 << " ";
00298    os << endMarker << "\n";
00299    os.precision(pr);
00300    return os;
00301 }
00302 
00303 std::vector<unsigned long> MTwistEngine::put () const {
00304   std::vector<unsigned long> v;
00305   v.push_back (engineIDulong<MTwistEngine>());
00306   for (int i=0; i<624; ++i) {
00307      v.push_back(static_cast<unsigned long>(mt[i]));
00308   }
00309   v.push_back(count624); 
00310   return v;
00311 }
00312 
00313 std::istream &  MTwistEngine::get ( std::istream& is )
00314 {
00315   char beginMarker [MarkerLen];
00316   is >> std::ws;
00317   is.width(MarkerLen);  // causes the next read to the char* to be <=
00318                         // that many bytes, INCLUDING A TERMINATION \0 
00319                         // (Stroustrup, section 21.3.2)
00320   is >> beginMarker;
00321   if (strcmp(beginMarker,"MTwistEngine-begin")) {
00322      is.clear(std::ios::badbit | is.rdstate());
00323      std::cerr << "\nInput stream mispositioned or"
00324                << "\nMTwistEngine state description missing or"
00325                << "\nwrong engine type found." << std::endl;
00326      return is;
00327    }
00328   return getState(is);
00329 }
00330 
00331 std::string MTwistEngine::beginTag ( )  { 
00332   return "MTwistEngine-begin"; 
00333 }
00334 
00335 std::istream &  MTwistEngine::getState ( std::istream& is )
00336 {
00337   char endMarker   [MarkerLen];
00338   is >> theSeed;
00339   for (int i=0; i<624; ++i)  is >> mt[i];
00340   is >> count624;
00341   is >> std::ws;
00342   is.width(MarkerLen);  
00343   is >> endMarker;
00344   if (strcmp(endMarker,"MTwistEngine-end")) {
00345      is.clear(std::ios::badbit | is.rdstate());
00346      std::cerr << "\nMTwistEngine state description incomplete."
00347                << "\nInput stream is probably mispositioned now." << std::endl;
00348      return is;
00349    }
00350    return is;
00351 }
00352 
00353 bool MTwistEngine::get (const std::vector<unsigned long> & v) {
00354   if ((v[0] & 0xffffffffUL) != engineIDulong<MTwistEngine>()) {
00355     std::cerr << 
00356         "\nMTwistEngine get:state vector has wrong ID word - state unchanged\n";
00357     return false;
00358   }
00359   return getState(v);
00360 }
00361 
00362 bool MTwistEngine::getState (const std::vector<unsigned long> & v) {
00363   if (v.size() != VECTOR_STATE_SIZE ) {
00364     std::cerr << 
00365         "\nMTwistEngine get:state vector has wrong length - state unchanged\n";
00366     return false;
00367   }
00368   for (int i=0; i<624; ++i) {
00369      mt[i]=v[i+1];
00370   }
00371   count624 = v[625];
00372   return true;
00373 }
00374 
00375 }  // namespace CLHEP