00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042 #include "CLHEP/Random/Random.h"
00043 #include "CLHEP/Random/MTwistEngine.h"
00044 #include "CLHEP/Random/engineIDulong.h"
00045 #include <string.h>
00046 #include <cstdlib>
00047
00048 namespace CLHEP {
00049
00050 static const int MarkerLen = 64;
00051
00052 std::string MTwistEngine::name() const {return "MTwistEngine";}
00053
00054 int MTwistEngine::numEngines = 0;
00055 int MTwistEngine::maxIndex = 215;
00056
00057 MTwistEngine::MTwistEngine()
00058 : HepRandomEngine()
00059 {
00060 int cycle = std::abs(int(numEngines/maxIndex));
00061 int curIndex = std::abs(int(numEngines%maxIndex));
00062 long mask = ((cycle & 0x007fffff) << 8);
00063 long seedlist[2];
00064 HepRandom::getTheTableSeeds( seedlist, curIndex );
00065 seedlist[0] = (seedlist[0])^mask;
00066 seedlist[1] = 0;
00067 setSeeds( seedlist, numEngines );
00068 count624=0;
00069 ++numEngines;
00070 for( int i=0; i < 2000; ++i ) flat();
00071 }
00072
00073 MTwistEngine::MTwistEngine(long seed)
00074 : HepRandomEngine()
00075 {
00076 long seedlist[2]={seed,17587};
00077 setSeeds( seedlist, 0 );
00078 count624=0;
00079 for( int i=0; i < 2000; ++i ) flat();
00080 }
00081
00082 MTwistEngine::MTwistEngine(int rowIndex, int colIndex)
00083 : HepRandomEngine()
00084 {
00085 int cycle = std::abs(int(rowIndex/maxIndex));
00086 int row = std::abs(int(rowIndex%maxIndex));
00087 int col = std::abs(int(colIndex%2));
00088 long mask = (( cycle & 0x000007ff ) << 20 );
00089 long seedlist[2];
00090 HepRandom::getTheTableSeeds( seedlist, row );
00091 seedlist[0] = (seedlist[col])^mask;
00092 seedlist[1] = 690691;
00093 setSeeds(seedlist, 4444772);
00094 count624=0;
00095 for( int i=0; i < 2000; ++i ) flat();
00096 }
00097
00098 MTwistEngine::MTwistEngine( std::istream& is )
00099 : HepRandomEngine()
00100 {
00101 is >> *this;
00102 }
00103
00104 MTwistEngine::~MTwistEngine() {}
00105
00106 double MTwistEngine::flat() {
00107 unsigned int y;
00108
00109 if( count624 >= N ) {
00110 register int i;
00111
00112 for( i=0; i < NminusM; ++i ) {
00113 y = (mt[i] & 0x80000000) | (mt[i+1] & 0x7fffffff);
00114 mt[i] = mt[i+M] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
00115 }
00116
00117 for( ; i < N-1 ; ++i ) {
00118 y = (mt[i] & 0x80000000) | (mt[i+1] & 0x7fffffff);
00119 mt[i] = mt[i-NminusM] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
00120 }
00121
00122 y = (mt[i] & 0x80000000) | (mt[0] & 0x7fffffff);
00123 mt[i] = mt[M-1] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
00124
00125 count624 = 0;
00126 }
00127
00128 y = mt[count624];
00129 y ^= ( y >> 11);
00130 y ^= ((y << 7 ) & 0x9d2c5680);
00131 y ^= ((y << 15) & 0xefc60000);
00132 y ^= ( y >> 18);
00133
00134 return y * twoToMinus_32() +
00135 (mt[count624++] >> 11) * twoToMinus_53() +
00136 nearlyTwoToMinus_54();
00137 }
00138
00139 void MTwistEngine::flatArray( const int size, double *vect ) {
00140 for( int i=0; i < size; ++i) vect[i] = flat();
00141 }
00142
00143 void MTwistEngine::setSeed(long seed, int k) {
00144
00145
00146
00147
00148
00149
00150
00151 theSeed = seed ? seed : 4357;
00152 int mti;
00153 const int N1=624;
00154 mt[0] = (unsigned int) (theSeed&0xffffffffUL);
00155 for (mti=1; mti<N1; mti++) {
00156 mt[mti] = (1812433253UL * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti);
00157
00158
00159
00160
00161 mt[mti] &= 0xffffffffUL;
00162
00163 }
00164 for( int i=1; i < 624; ++i ) {
00165 mt[i] ^= k;
00166 }
00167
00168
00169 }
00170
00171 void MTwistEngine::setSeeds(const long *seeds, int k) {
00172 setSeed( (*seeds ? *seeds : 43571346), k );
00173 int i;
00174 for( i=1; i < 624; ++i ) {
00175 mt[i] = ( seeds[1] + mt[i] ) & 0xffffffff;
00176 }
00177 theSeeds = seeds;
00178 }
00179
00180 void MTwistEngine::saveStatus( const char filename[] ) const
00181 {
00182 std::ofstream outFile( filename, std::ios::out ) ;
00183 if (!outFile.bad()) {
00184 outFile << theSeed << std::endl;
00185 for (int i=0; i<624; ++i) outFile <<std::setprecision(20) << mt[i] << " ";
00186 outFile << std::endl;
00187 outFile << count624 << std::endl;
00188 }
00189 }
00190
00191 void MTwistEngine::restoreStatus( const char filename[] )
00192 {
00193 std::ifstream inFile( filename, std::ios::in);
00194 if (!checkFile ( inFile, filename, engineName(), "restoreStatus" )) {
00195 std::cerr << " -- Engine state remains unchanged\n";
00196 return;
00197 }
00198
00199 if (!inFile.bad() && !inFile.eof()) {
00200 inFile >> theSeed;
00201 for (int i=0; i<624; ++i) inFile >> mt[i];
00202 inFile >> count624;
00203 }
00204 }
00205
00206 void MTwistEngine::showStatus() const
00207 {
00208 std::cout << std::endl;
00209 std::cout << "--------- MTwist engine status ---------" << std::endl;
00210 std::cout << std::setprecision(20);
00211 std::cout << " Initial seed = " << theSeed << std::endl;
00212 std::cout << " Current index = " << count624 << std::endl;
00213 std::cout << " Array status mt[] = " << std::endl;
00214 for (int i=0; i<624; i+=5) {
00215 std::cout << mt[i] << " " << mt[i+1] << " " << mt[i+2] << " "
00216 << mt[i+3] << " " << mt[i+4] << std::endl;
00217 }
00218 std::cout << "----------------------------------------" << std::endl;
00219 }
00220
00221 MTwistEngine::operator float() {
00222 unsigned int y;
00223
00224 if( count624 >= N ) {
00225 register int i;
00226
00227 for( i=0; i < NminusM; ++i ) {
00228 y = (mt[i] & 0x80000000) | (mt[i+1] & 0x7fffffff);
00229 mt[i] = mt[i+M] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
00230 }
00231
00232 for( ; i < N-1 ; ++i ) {
00233 y = (mt[i] & 0x80000000) | (mt[i+1] & 0x7fffffff);
00234 mt[i] = mt[i-NminusM] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
00235 }
00236
00237 y = (mt[i] & 0x80000000) | (mt[0] & 0x7fffffff);
00238 mt[i] = mt[M-1] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
00239
00240 count624 = 0;
00241 }
00242
00243 y = mt[count624++];
00244 y ^= ( y >> 11);
00245 y ^= ((y << 7 ) & 0x9d2c5680);
00246 y ^= ((y << 15) & 0xefc60000);
00247 y ^= ( y >> 18);
00248
00249 return (float)(y * twoToMinus_32());
00250 }
00251
00252 MTwistEngine::operator unsigned int() {
00253 unsigned int y;
00254
00255 if( count624 >= N ) {
00256 register int i;
00257
00258 for( i=0; i < NminusM; ++i ) {
00259 y = (mt[i] & 0x80000000) | (mt[i+1] & 0x7fffffff);
00260 mt[i] = mt[i+M] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
00261 }
00262
00263 for( ; i < N-1 ; ++i ) {
00264 y = (mt[i] & 0x80000000) | (mt[i+1] & 0x7fffffff);
00265 mt[i] = mt[i-NminusM] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
00266 }
00267
00268 y = (mt[i] & 0x80000000) | (mt[0] & 0x7fffffff);
00269 mt[i] = mt[M-1] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
00270
00271 count624 = 0;
00272 }
00273
00274 y = mt[count624++];
00275 y ^= ( y >> 11);
00276 y ^= ((y << 7 ) & 0x9d2c5680);
00277 y ^= ((y << 15) & 0xefc60000);
00278 y ^= ( y >> 18);
00279
00280 return y;
00281 }
00282
00283 std::ostream & MTwistEngine::put ( std::ostream& os ) const
00284 {
00285 char beginMarker[] = "MTwistEngine-begin";
00286 char endMarker[] = "MTwistEngine-end";
00287
00288 int pr = os.precision(20);
00289 os << " " << beginMarker << " ";
00290 os << theSeed << " ";
00291 for (int i=0; i<624; ++i) {
00292 os << mt[i] << "\n";
00293 }
00294 os << count624 << " ";
00295 os << endMarker << "\n";
00296 os.precision(pr);
00297 return os;
00298 }
00299
00300 std::vector<unsigned long> MTwistEngine::put () const {
00301 std::vector<unsigned long> v;
00302 v.push_back (engineIDulong<MTwistEngine>());
00303 for (int i=0; i<624; ++i) {
00304 v.push_back(static_cast<unsigned long>(mt[i]));
00305 }
00306 v.push_back(count624);
00307 return v;
00308 }
00309
00310 std::istream & MTwistEngine::get ( std::istream& is )
00311 {
00312 char beginMarker [MarkerLen];
00313 is >> std::ws;
00314 is.width(MarkerLen);
00315
00316
00317 is >> beginMarker;
00318 if (strcmp(beginMarker,"MTwistEngine-begin")) {
00319 is.clear(std::ios::badbit | is.rdstate());
00320 std::cerr << "\nInput stream mispositioned or"
00321 << "\nMTwistEngine state description missing or"
00322 << "\nwrong engine type found." << std::endl;
00323 return is;
00324 }
00325 return getState(is);
00326 }
00327
00328 std::string MTwistEngine::beginTag ( ) {
00329 return "MTwistEngine-begin";
00330 }
00331
00332 std::istream & MTwistEngine::getState ( std::istream& is )
00333 {
00334 char endMarker [MarkerLen];
00335 is >> theSeed;
00336 for (int i=0; i<624; ++i) is >> mt[i];
00337 is >> count624;
00338 is >> std::ws;
00339 is.width(MarkerLen);
00340 is >> endMarker;
00341 if (strcmp(endMarker,"MTwistEngine-end")) {
00342 is.clear(std::ios::badbit | is.rdstate());
00343 std::cerr << "\nMTwistEngine state description incomplete."
00344 << "\nInput stream is probably mispositioned now." << std::endl;
00345 return is;
00346 }
00347 return is;
00348 }
00349
00350 bool MTwistEngine::get (const std::vector<unsigned long> & v) {
00351 if ((v[0] & 0xffffffffUL) != engineIDulong<MTwistEngine>()) {
00352 std::cerr <<
00353 "\nMTwistEngine get:state vector has wrong ID word - state unchanged\n";
00354 return false;
00355 }
00356 return getState(v);
00357 }
00358
00359 bool MTwistEngine::getState (const std::vector<unsigned long> & v) {
00360 if (v.size() != VECTOR_STATE_SIZE ) {
00361 std::cerr <<
00362 "\nMTwistEngine get:state vector has wrong length - state unchanged\n";
00363 return false;
00364 }
00365 for (int i=0; i<624; ++i) {
00366 mt[i]=v[i+1];
00367 }
00368 count624 = v[625];
00369 return true;
00370 }
00371
00372 }