42 #include "CLHEP/Random/defs.h"
43 #include "CLHEP/Random/Random.h"
44 #include "CLHEP/Random/MTwistEngine.h"
45 #include "CLHEP/Random/engineIDulong.h"
51 static const int MarkerLen = 64;
55 int MTwistEngine::numEngines = 0;
56 int MTwistEngine::maxIndex = 215;
61 int cycle = std::abs(
int(numEngines/maxIndex));
62 int curIndex = std::abs(
int(numEngines%maxIndex));
63 long mask = ((cycle & 0x007fffff) << 8);
66 seedlist[0] = (seedlist[0])^mask;
71 for(
int i=0;
i < 2000; ++
i )
flat();
77 long seedlist[2]={seed,17587};
80 for(
int i=0;
i < 2000; ++
i )
flat();
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 );
92 seedlist[0] = (seedlist[col])^mask;
96 for(
int i=0;
i < 2000; ++
i )
flat();
110 if( count624 >= N ) {
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 );
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 );
123 y = (mt[
i] & 0x80000000) | (mt[0] & 0x7fffffff);
124 mt[
i] = mt[M-1] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
131 y ^= ((y << 7 ) & 0x9d2c5680);
132 y ^= ((y << 15) & 0xefc60000);
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);
162 mt[mti] &= 0xffffffffUL;
165 for(
int i=1;
i < 624; ++
i ) {
175 for(
i=1;
i < 624; ++
i ) {
176 mt[
i] = (
seeds[1] + mt[
i] ) & 0xffffffff;
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;
196 std::cerr <<
" -- Engine state remains unchanged\n";
200 if (!inFile.bad() && !inFile.eof()) {
202 for (
int i=0; i<624; ++i) inFile >> mt[
i];
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;
219 std::cout <<
"----------------------------------------" << std::endl;
222 MTwistEngine::operator float() {
225 if( count624 >=
N ) {
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 );
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 );
238 y = (mt[
i] & 0x80000000) | (mt[0] & 0x7fffffff);
239 mt[
i] = mt[M-1] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
246 y ^= ((y << 7 ) & 0x9d2c5680);
247 y ^= ((y << 15) & 0xefc60000);
250 return (
float)(y * twoToMinus_32());
253 MTwistEngine::operator
unsigned int() {
256 if( count624 >=
N ) {
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 );
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 );
269 y = (mt[
i] & 0x80000000) | (mt[0] & 0x7fffffff);
270 mt[
i] = mt[M-1] ^ (y >> 1) ^ ((y & 0x1) ? 0x9908b0df : 0x0 );
277 y ^= ((y << 7 ) & 0x9d2c5680);
278 y ^= ((y << 15) & 0xefc60000);
286 char beginMarker[] =
"MTwistEngine-begin";
287 char endMarker[] =
"MTwistEngine-end";
289 int pr = os.precision(20);
290 os <<
" " << beginMarker <<
" ";
292 for (
int i=0;
i<624; ++
i) {
295 os << count624 <<
" ";
296 os << endMarker <<
"\n";
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]));
307 v.push_back(count624);
313 char beginMarker [MarkerLen];
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;
330 return "MTwistEngine-begin";
335 char endMarker [MarkerLen];
337 for (
int i=0; i<624; ++i) is >> mt[
i];
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;
352 if ((
v[0] & 0xffffffffUL) != engineIDulong<MTwistEngine>()) {
354 "\nMTwistEngine get:state vector has wrong ID word - state unchanged\n";
363 "\nMTwistEngine get:state vector has wrong length - state unchanged\n";
366 for (
int i=0;
i<624; ++
i) {