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 #include "CLHEP/Random/Random.h"
00041 #include "CLHEP/Random/RanluxEngine.h"
00042 #include "CLHEP/Random/engineIDulong.h"
00043 #include <string.h>
00044 #include <cstdlib>
00045
00046 namespace CLHEP {
00047
00048
00049 static const int MarkerLen = 64;
00050
00051 std::string RanluxEngine::name() const {return "RanluxEngine";}
00052
00053
00054 int RanluxEngine::numEngines = 0;
00055
00056
00057 int RanluxEngine::maxIndex = 215;
00058
00059 RanluxEngine::RanluxEngine(long seed, int lux)
00060 : HepRandomEngine()
00061 {
00062 long seedlist[2]={0,0};
00063
00064 luxury = lux;
00065 setSeed(seed, luxury);
00066
00067
00068 seedlist[0]=theSeed;
00069 seedlist[1]=0;
00070 setSeeds(seedlist, luxury);
00071 }
00072
00073 RanluxEngine::RanluxEngine()
00074 : HepRandomEngine()
00075 {
00076 long seed;
00077 long seedlist[2]={0,0};
00078
00079 luxury = 3;
00080 int cycle = std::abs(int(numEngines/maxIndex));
00081 int curIndex = std::abs(int(numEngines%maxIndex));
00082 numEngines +=1;
00083 long mask = ((cycle & 0x007fffff) << 8);
00084 HepRandom::getTheTableSeeds( seedlist, curIndex );
00085 seed = seedlist[0]^mask;
00086 setSeed(seed, luxury);
00087
00088
00089 seedlist[0]=theSeed;
00090 seedlist[1]=0;
00091 setSeeds(seedlist, luxury);
00092 }
00093
00094 RanluxEngine::RanluxEngine(int rowIndex, int colIndex, int lux)
00095 : HepRandomEngine()
00096 {
00097 long seed;
00098 long seedlist[2]={0,0};
00099
00100 luxury = lux;
00101 int cycle = std::abs(int(rowIndex/maxIndex));
00102 int row = std::abs(int(rowIndex%maxIndex));
00103 int col = std::abs(int(colIndex%2));
00104 long mask = (( cycle & 0x000007ff ) << 20 );
00105 HepRandom::getTheTableSeeds( seedlist, row );
00106 seed = ( seedlist[col] )^mask;
00107 setSeed(seed, luxury);
00108
00109
00110 seedlist[0]=theSeed;
00111 seedlist[1]=0;
00112 setSeeds(seedlist, luxury);
00113 }
00114
00115 RanluxEngine::RanluxEngine( std::istream& is )
00116 : HepRandomEngine()
00117 {
00118 is >> *this;
00119 }
00120
00121 RanluxEngine::~RanluxEngine() {}
00122
00123 void RanluxEngine::setSeed(long seed, int lux) {
00124
00125
00126
00127
00128
00129
00130
00131 const int ecuyer_a = 53668;
00132 const int ecuyer_b = 40014;
00133 const int ecuyer_c = 12211;
00134 const int ecuyer_d = 2147483563;
00135
00136 const int lux_levels[5] = {0,24,73,199,365};
00137
00138 long int_seed_table[24];
00139 long next_seed = seed;
00140 long k_multiple;
00141 int i;
00142
00143
00144
00145
00146 theSeed = seed;
00147 if( (lux > 4)||(lux < 0) ){
00148 if(lux >= 24){
00149 nskip = lux - 24;
00150 }else{
00151 nskip = lux_levels[3];
00152 }
00153 }else{
00154 luxury = lux;
00155 nskip = lux_levels[luxury];
00156 }
00157
00158
00159 for(i = 0;i != 24;i++){
00160 k_multiple = next_seed / ecuyer_a;
00161 next_seed = ecuyer_b * (next_seed - k_multiple * ecuyer_a)
00162 - k_multiple * ecuyer_c ;
00163 if(next_seed < 0)next_seed += ecuyer_d;
00164 int_seed_table[i] = next_seed % int_modulus;
00165 }
00166
00167 for(i = 0;i != 24;i++)
00168 float_seed_table[i] = int_seed_table[i] * mantissa_bit_24();
00169
00170 i_lag = 23;
00171 j_lag = 9;
00172 carry = 0. ;
00173
00174 if( float_seed_table[23] == 0. ) carry = mantissa_bit_24();
00175
00176 count24 = 0;
00177 }
00178
00179 void RanluxEngine::setSeeds(const long *seeds, int lux) {
00180
00181 const int ecuyer_a = 53668;
00182 const int ecuyer_b = 40014;
00183 const int ecuyer_c = 12211;
00184 const int ecuyer_d = 2147483563;
00185
00186 const int lux_levels[5] = {0,24,73,199,365};
00187 int i;
00188 long int_seed_table[24];
00189 long k_multiple,next_seed;
00190 const long *seedptr;
00191
00192 theSeeds = seeds;
00193 seedptr = seeds;
00194
00195 if(seeds == 0){
00196 setSeed(theSeed,lux);
00197 theSeeds = &theSeed;
00198 return;
00199 }
00200
00201 theSeed = *seeds;
00202
00203
00204
00205
00206 if( (lux > 4)||(lux < 0) ){
00207 if(lux >= 24){
00208 nskip = lux - 24;
00209 }else{
00210 nskip = lux_levels[3];
00211 }
00212 }else{
00213 luxury = lux;
00214 nskip = lux_levels[luxury];
00215 }
00216
00217 for( i = 0;(i != 24)&&(*seedptr != 0);i++){
00218 int_seed_table[i] = *seedptr % int_modulus;
00219 seedptr++;
00220 }
00221
00222 if(i != 24){
00223 next_seed = int_seed_table[i-1];
00224 for(;i != 24;i++){
00225 k_multiple = next_seed / ecuyer_a;
00226 next_seed = ecuyer_b * (next_seed - k_multiple * ecuyer_a)
00227 - k_multiple * ecuyer_c ;
00228 if(next_seed < 0)next_seed += ecuyer_d;
00229 int_seed_table[i] = next_seed % int_modulus;
00230 }
00231 }
00232
00233 for(i = 0;i != 24;i++)
00234 float_seed_table[i] = int_seed_table[i] * mantissa_bit_24();
00235
00236 i_lag = 23;
00237 j_lag = 9;
00238 carry = 0. ;
00239
00240 if( float_seed_table[23] == 0. ) carry = mantissa_bit_24();
00241
00242 count24 = 0;
00243 }
00244
00245 void RanluxEngine::saveStatus( const char filename[] ) const
00246 {
00247 std::ofstream outFile( filename, std::ios::out ) ;
00248 if (!outFile.bad()) {
00249 outFile << "Uvec\n";
00250 std::vector<unsigned long> v = put();
00251 for (unsigned int i=0; i<v.size(); ++i) {
00252 outFile << v[i] << "\n";
00253 }
00254 }
00255 }
00256
00257 void RanluxEngine::restoreStatus( const char filename[] )
00258 {
00259 std::ifstream inFile( filename, std::ios::in);
00260 if (!checkFile ( inFile, filename, engineName(), "restoreStatus" )) {
00261 std::cerr << " -- Engine state remains unchanged\n";
00262 return;
00263 }
00264 if ( possibleKeywordInput ( inFile, "Uvec", theSeed ) ) {
00265 std::vector<unsigned long> v;
00266 unsigned long xin;
00267 for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) {
00268 inFile >> xin;
00269 if (!inFile) {
00270 inFile.clear(std::ios::badbit | inFile.rdstate());
00271 std::cerr << "\nRanluxEngine state (vector) description improper."
00272 << "\nrestoreStatus has failed."
00273 << "\nInput stream is probably mispositioned now." << std::endl;
00274 return;
00275 }
00276 v.push_back(xin);
00277 }
00278 getState(v);
00279 return;
00280 }
00281
00282 if (!inFile.bad() && !inFile.eof()) {
00283
00284 for (int i=0; i<24; ++i)
00285 inFile >> float_seed_table[i];
00286 inFile >> i_lag; inFile >> j_lag;
00287 inFile >> carry; inFile >> count24;
00288 inFile >> luxury; inFile >> nskip;
00289 }
00290 }
00291
00292 void RanluxEngine::showStatus() const
00293 {
00294 std::cout << std::endl;
00295 std::cout << "--------- Ranlux engine status ---------" << std::endl;
00296 std::cout << " Initial seed = " << theSeed << std::endl;
00297 std::cout << " float_seed_table[] = ";
00298 for (int i=0; i<24; ++i)
00299 std::cout << float_seed_table[i] << " ";
00300 std::cout << std::endl;
00301 std::cout << " i_lag = " << i_lag << ", j_lag = " << j_lag << std::endl;
00302 std::cout << " carry = " << carry << ", count24 = " << count24 << std::endl;
00303 std::cout << " luxury = " << luxury << " nskip = " << nskip << std::endl;
00304 std::cout << "----------------------------------------" << std::endl;
00305 }
00306
00307 double RanluxEngine::flat() {
00308
00309 float next_random;
00310 float uni;
00311 int i;
00312
00313 uni = float_seed_table[j_lag] - float_seed_table[i_lag] - carry;
00314 if(uni < 0. ){
00315 uni += 1.0;
00316 carry = mantissa_bit_24();
00317 }else{
00318 carry = 0.;
00319 }
00320
00321 float_seed_table[i_lag] = uni;
00322 i_lag --;
00323 j_lag --;
00324 if(i_lag < 0) i_lag = 23;
00325 if(j_lag < 0) j_lag = 23;
00326
00327 if( uni < mantissa_bit_12() ){
00328 uni += mantissa_bit_24() * float_seed_table[j_lag];
00329 if( uni == 0) uni = mantissa_bit_24() * mantissa_bit_24();
00330 }
00331 next_random = uni;
00332 count24 ++;
00333
00334
00335
00336
00337 if(count24 == 24 ){
00338 count24 = 0;
00339 for( i = 0; i != nskip ; i++){
00340 uni = float_seed_table[j_lag] - float_seed_table[i_lag] - carry;
00341 if(uni < 0. ){
00342 uni += 1.0;
00343 carry = mantissa_bit_24();
00344 }else{
00345 carry = 0.;
00346 }
00347 float_seed_table[i_lag] = uni;
00348 i_lag --;
00349 j_lag --;
00350 if(i_lag < 0)i_lag = 23;
00351 if(j_lag < 0) j_lag = 23;
00352 }
00353 }
00354 return (double) next_random;
00355 }
00356
00357 void RanluxEngine::flatArray(const int size, double* vect)
00358 {
00359 float next_random;
00360 float uni;
00361 int i;
00362 int index;
00363
00364 for (index=0; index<size; ++index) {
00365 uni = float_seed_table[j_lag] - float_seed_table[i_lag] - carry;
00366 if(uni < 0. ){
00367 uni += 1.0;
00368 carry = mantissa_bit_24();
00369 }else{
00370 carry = 0.;
00371 }
00372
00373 float_seed_table[i_lag] = uni;
00374 i_lag --;
00375 j_lag --;
00376 if(i_lag < 0) i_lag = 23;
00377 if(j_lag < 0) j_lag = 23;
00378
00379 if( uni < mantissa_bit_12() ){
00380 uni += mantissa_bit_24() * float_seed_table[j_lag];
00381 if( uni == 0) uni = mantissa_bit_24() * mantissa_bit_24();
00382 }
00383 next_random = uni;
00384 vect[index] = (double)next_random;
00385 count24 ++;
00386
00387
00388
00389
00390 if(count24 == 24 ){
00391 count24 = 0;
00392 for( i = 0; i != nskip ; i++){
00393 uni = float_seed_table[j_lag] - float_seed_table[i_lag] - carry;
00394 if(uni < 0. ){
00395 uni += 1.0;
00396 carry = mantissa_bit_24();
00397 }else{
00398 carry = 0.;
00399 }
00400 float_seed_table[i_lag] = uni;
00401 i_lag --;
00402 j_lag --;
00403 if(i_lag < 0)i_lag = 23;
00404 if(j_lag < 0) j_lag = 23;
00405 }
00406 }
00407 }
00408 }
00409
00410 RanluxEngine::operator unsigned int() {
00411 return ((unsigned int)(flat() * exponent_bit_32()) & 0xffffffff) |
00412 (((unsigned int)(float_seed_table[i_lag]*exponent_bit_32())>>16) & 0xff);
00413
00414
00415 }
00416
00417 std::ostream & RanluxEngine::put ( std::ostream& os ) const
00418 {
00419 char beginMarker[] = "RanluxEngine-begin";
00420 os << beginMarker << "\nUvec\n";
00421 std::vector<unsigned long> v = put();
00422 for (unsigned int i=0; i<v.size(); ++i) {
00423 os << v[i] << "\n";
00424 }
00425 return os;
00426 }
00427
00428 std::vector<unsigned long> RanluxEngine::put () const {
00429 std::vector<unsigned long> v;
00430 v.push_back (engineIDulong<RanluxEngine>());
00431 for (int i=0; i<24; ++i) {
00432 v.push_back
00433 (static_cast<unsigned long>(float_seed_table[i]/mantissa_bit_24()));
00434 }
00435 v.push_back(static_cast<unsigned long>(i_lag));
00436 v.push_back(static_cast<unsigned long>(j_lag));
00437 v.push_back(static_cast<unsigned long>(carry/mantissa_bit_24()));
00438 v.push_back(static_cast<unsigned long>(count24));
00439 v.push_back(static_cast<unsigned long>(luxury));
00440 v.push_back(static_cast<unsigned long>(nskip));
00441 return v;
00442 }
00443
00444 std::istream & RanluxEngine::get ( std::istream& is )
00445 {
00446 char beginMarker [MarkerLen];
00447 is >> std::ws;
00448 is.width(MarkerLen);
00449
00450
00451 is >> beginMarker;
00452 if (strcmp(beginMarker,"RanluxEngine-begin")) {
00453 is.clear(std::ios::badbit | is.rdstate());
00454 std::cerr << "\nInput stream mispositioned or"
00455 << "\nRanluxEngine state description missing or"
00456 << "\nwrong engine type found." << std::endl;
00457 return is;
00458 }
00459 return getState(is);
00460 }
00461
00462 std::string RanluxEngine::beginTag ( ) {
00463 return "RanluxEngine-begin";
00464 }
00465
00466 std::istream & RanluxEngine::getState ( std::istream& is )
00467 {
00468 if ( possibleKeywordInput ( is, "Uvec", theSeed ) ) {
00469 std::vector<unsigned long> v;
00470 unsigned long uu;
00471 for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) {
00472 is >> uu;
00473 if (!is) {
00474 is.clear(std::ios::badbit | is.rdstate());
00475 std::cerr << "\nRanluxEngine state (vector) description improper."
00476 << "\ngetState() has failed."
00477 << "\nInput stream is probably mispositioned now." << std::endl;
00478 return is;
00479 }
00480 v.push_back(uu);
00481 }
00482 getState(v);
00483 return (is);
00484 }
00485
00486
00487
00488 char endMarker [MarkerLen];
00489 for (int i=0; i<24; ++i) {
00490 is >> float_seed_table[i];
00491 }
00492 is >> i_lag; is >> j_lag;
00493 is >> carry; is >> count24;
00494 is >> luxury; is >> nskip;
00495 is >> std::ws;
00496 is.width(MarkerLen);
00497 is >> endMarker;
00498 if (strcmp(endMarker,"RanluxEngine-end")) {
00499 is.clear(std::ios::badbit | is.rdstate());
00500 std::cerr << "\nRanluxEngine state description incomplete."
00501 << "\nInput stream is probably mispositioned now." << std::endl;
00502 return is;
00503 }
00504 return is;
00505 }
00506
00507 bool RanluxEngine::get (const std::vector<unsigned long> & v) {
00508 if ((v[0] & 0xffffffffUL) != engineIDulong<RanluxEngine>()) {
00509 std::cerr <<
00510 "\nRanluxEngine get:state vector has wrong ID word - state unchanged\n";
00511 return false;
00512 }
00513 return getState(v);
00514 }
00515
00516 bool RanluxEngine::getState (const std::vector<unsigned long> & v) {
00517 if (v.size() != VECTOR_STATE_SIZE ) {
00518 std::cerr <<
00519 "\nRanluxEngine get:state vector has wrong length - state unchanged\n";
00520 return false;
00521 }
00522 for (int i=0; i<24; ++i) {
00523 float_seed_table[i] = v[i+1]*mantissa_bit_24();
00524 }
00525 i_lag = v[25];
00526 j_lag = v[26];
00527 carry = v[27]*mantissa_bit_24();
00528 count24 = v[28];
00529 luxury = v[29];
00530 nskip = v[30];
00531 return true;
00532 }
00533
00534 }