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
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056 #include "CLHEP/Random/Random.h"
00057 #include "CLHEP/Random/Ranlux64Engine.h"
00058 #include "CLHEP/Random/engineIDulong.h"
00059 #include "CLHEP/Random/DoubConv.h"
00060 #include <string.h>
00061 #include <cstdlib>
00062 #include <limits>
00063
00064 namespace CLHEP {
00065
00066 static const int MarkerLen = 64;
00067
00068
00069
00070 int Ranlux64Engine::numEngines = 0;
00071
00072
00073 int Ranlux64Engine::maxIndex = 215;
00074
00075 #ifndef WIN32
00076 namespace detail {
00077
00078 template< std::size_t n,
00079 bool = n < std::size_t(std::numeric_limits<unsigned long>::digits) >
00080 struct do_right_shift;
00081 template< std::size_t n >
00082 struct do_right_shift<n,true>
00083 {
00084 unsigned long operator()(unsigned long value) { return value >> n; }
00085 };
00086 template< std::size_t n >
00087 struct do_right_shift<n,false>
00088 {
00089 unsigned long operator()(unsigned long) { return 0ul; }
00090 };
00091
00092 template< std::size_t nbits >
00093 unsigned long rshift( unsigned long value )
00094 { return do_right_shift<nbits>()(value); }
00095
00096 }
00097 #endif
00098
00099 std::string Ranlux64Engine::name() const {return "Ranlux64Engine";}
00100
00101 Ranlux64Engine::Ranlux64Engine()
00102 : HepRandomEngine()
00103 {
00104 luxury = 1;
00105 int cycle = std::abs(int(numEngines/maxIndex));
00106 int curIndex = std::abs(int(numEngines%maxIndex));
00107 numEngines +=1;
00108 long mask = ((cycle & 0x007fffff) << 8);
00109 long seedlist[2];
00110 HepRandom::getTheTableSeeds( seedlist, curIndex );
00111 seedlist[0] ^= mask;
00112 seedlist[1] = 0;
00113
00114 setSeeds(seedlist, luxury);
00115 advance ( 8 );
00116
00117
00118 }
00119
00120 Ranlux64Engine::Ranlux64Engine(long seed, int lux)
00121 : HepRandomEngine()
00122 {
00123 luxury = lux;
00124 long seedlist[2]={seed,0};
00125 setSeeds(seedlist, lux);
00126 advance ( 2*lux + 1 );
00127
00128 }
00129
00130 Ranlux64Engine::Ranlux64Engine(int rowIndex, int, int lux)
00131 : HepRandomEngine()
00132 {
00133 luxury = lux;
00134 int cycle = std::abs(int(rowIndex/maxIndex));
00135 int row = std::abs(int(rowIndex%maxIndex));
00136 long mask = (( cycle & 0x000007ff ) << 20 );
00137 long seedlist[2];
00138 HepRandom::getTheTableSeeds( seedlist, row );
00139 seedlist[0] ^= mask;
00140 seedlist[1]= 0;
00141 setSeeds(seedlist, lux);
00142 }
00143
00144 Ranlux64Engine::Ranlux64Engine( std::istream& is )
00145 : HepRandomEngine()
00146 {
00147 is >> *this;
00148 }
00149
00150 Ranlux64Engine::~Ranlux64Engine() {}
00151
00152 double Ranlux64Engine::flat() {
00153
00154
00155
00156
00157
00158 if (index <= 0) update();
00159 return randoms[--index] + twoToMinus_49();
00160 }
00161
00162 void Ranlux64Engine::update() {
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180 advance(pDozens);
00181
00182
00183
00184
00185
00186
00187
00188
00189 register double y1;
00190
00191 if ( endIters == 1 ) {
00192 y1 = randoms[ 4] - randoms[11] - carry;
00193 if ( y1 < 0.0 ) {
00194 y1 += 1.0;
00195 carry = twoToMinus_48();
00196 } else {
00197 carry = 0.0;
00198 }
00199 randoms[11] = randoms[10];
00200 randoms[10] = randoms[ 9];
00201 randoms[ 9] = randoms[ 8];
00202 randoms[ 8] = randoms[ 7];
00203 randoms[ 7] = randoms[ 6];
00204 randoms[ 6] = randoms[ 5];
00205 randoms[ 5] = randoms[ 4];
00206 randoms[ 4] = randoms[ 3];
00207 randoms[ 3] = randoms[ 2];
00208 randoms[ 2] = randoms[ 1];
00209 randoms[ 1] = randoms[ 0];
00210 randoms[ 0] = y1;
00211
00212 } else {
00213
00214 int m, nr, ns;
00215 for ( m = 0, nr = 11, ns = 4; m < endIters; ++m, --nr ) {
00216 y1 = randoms [ns] - randoms[nr] - carry;
00217 if ( y1 < 0.0 ) {
00218 y1 += 1.0;
00219 carry = twoToMinus_48();
00220 } else {
00221 carry = 0.0;
00222 }
00223 randoms[nr] = y1;
00224 --ns;
00225 if ( ns < 0 ) {
00226 ns = 11;
00227 }
00228 }
00229
00230 double temp[12];
00231 for (m=0; m<12; m++) {
00232 temp[m]=randoms[m];
00233 }
00234
00235 ns = 11 - endIters;
00236 for (m=11; m>=0; --m) {
00237 randoms[m] = temp[ns];
00238 --ns;
00239 if ( ns < 0 ) {
00240 ns = 11;
00241 }
00242 }
00243
00244 }
00245
00246
00247
00248 index = 11;
00249
00250 }
00251
00252 void Ranlux64Engine::advance(int dozens) {
00253
00254 register double y1, y2, y3;
00255 register double cValue = twoToMinus_48();
00256 register double zero = 0.0;
00257 register double one = 1.0;
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277 int k;
00278 for ( k = dozens; k > 0; --k ) {
00279
00280 y1 = randoms[ 4] - randoms[11] - carry;
00281
00282 y2 = randoms[ 3] - randoms[10];
00283 if ( y1 < zero ) {
00284 y1 += one;
00285 y2 -= cValue;
00286 }
00287 randoms[11] = y1;
00288
00289 y3 = randoms[ 2] - randoms[ 9];
00290 if ( y2 < zero ) {
00291 y2 += one;
00292 y3 -= cValue;
00293 }
00294 randoms[10] = y2;
00295
00296 y1 = randoms[ 1] - randoms[ 8];
00297 if ( y3 < zero ) {
00298 y3 += one;
00299 y1 -= cValue;
00300 }
00301 randoms[ 9] = y3;
00302
00303 y2 = randoms[ 0] - randoms[ 7];
00304 if ( y1 < zero ) {
00305 y1 += one;
00306 y2 -= cValue;
00307 }
00308 randoms[ 8] = y1;
00309
00310 y3 = randoms[11] - randoms[ 6];
00311 if ( y2 < zero ) {
00312 y2 += one;
00313 y3 -= cValue;
00314 }
00315 randoms[ 7] = y2;
00316
00317 y1 = randoms[10] - randoms[ 5];
00318 if ( y3 < zero ) {
00319 y3 += one;
00320 y1 -= cValue;
00321 }
00322 randoms[ 6] = y3;
00323
00324 y2 = randoms[ 9] - randoms[ 4];
00325 if ( y1 < zero ) {
00326 y1 += one;
00327 y2 -= cValue;
00328 }
00329 randoms[ 5] = y1;
00330
00331 y3 = randoms[ 8] - randoms[ 3];
00332 if ( y2 < zero ) {
00333 y2 += one;
00334 y3 -= cValue;
00335 }
00336 randoms[ 4] = y2;
00337
00338 y1 = randoms[ 7] - randoms[ 2];
00339 if ( y3 < zero ) {
00340 y3 += one;
00341 y1 -= cValue;
00342 }
00343 randoms[ 3] = y3;
00344
00345 y2 = randoms[ 6] - randoms[ 1];
00346 if ( y1 < zero ) {
00347 y1 += one;
00348 y2 -= cValue;
00349 }
00350 randoms[ 2] = y1;
00351
00352 y3 = randoms[ 5] - randoms[ 0];
00353 if ( y2 < zero ) {
00354 y2 += one;
00355 y3 -= cValue;
00356 }
00357 randoms[ 1] = y2;
00358
00359 if ( y3 < zero ) {
00360 y3 += one;
00361 carry = cValue;
00362 }
00363 randoms[ 0] = y3;
00364
00365 }
00366
00367 }
00368
00369 void Ranlux64Engine::flatArray(const int size, double* vect) {
00370 for( int i=0; i < size; ++i ) {
00371 vect[i] = flat();
00372 }
00373 }
00374
00375 void Ranlux64Engine::setSeed(long seed, int lux) {
00376
00377
00378
00379
00380
00381
00382
00383 const int ecuyer_a(53668);
00384 const int ecuyer_b(40014);
00385 const int ecuyer_c(12211);
00386 const int ecuyer_d(2147483563);
00387
00388 const int lux_levels[3] = {109, 202, 397};
00389 theSeed = seed;
00390
00391 if( (lux > 2)||(lux < 0) ){
00392 pDiscard = (lux >= 12) ? (lux-12) : lux_levels[1];
00393 }else{
00394 pDiscard = lux_levels[luxury];
00395 }
00396 pDozens = pDiscard / 12;
00397 endIters = pDiscard % 12;
00398
00399 long init_table[24];
00400 long next_seed = seed;
00401 long k_multiple;
00402 int i;
00403 next_seed &= 0xffffffff;
00404 while( next_seed >= ecuyer_d ) {
00405 next_seed -= ecuyer_d;
00406 }
00407
00408 for(i = 0;i != 24;i++){
00409 k_multiple = next_seed / ecuyer_a;
00410 next_seed = ecuyer_b * (next_seed - k_multiple * ecuyer_a)
00411 - k_multiple * ecuyer_c;
00412 if(next_seed < 0) {
00413 next_seed += ecuyer_d;
00414 }
00415 next_seed &= 0xffffffff;
00416 init_table[i] = next_seed;
00417 }
00418
00419 if( sizeof(long) >= 8 ) {
00420 long topbits1, topbits2;
00421 #ifdef WIN32
00422 topbits1 = ( seed >> 32) & 0xffff ;
00423 topbits2 = ( seed >> 48) & 0xffff ;
00424 #else
00425 topbits1 = detail::rshift<32>(seed) & 0xffff ;
00426 topbits2 = detail::rshift<48>(seed) & 0xffff ;
00427 #endif
00428 init_table[0] ^= topbits1;
00429 init_table[2] ^= topbits2;
00430
00431
00432 }
00433
00434 for(i = 0;i < 12; i++){
00435 randoms[i] = (init_table[2*i ] ) * 2.0 * twoToMinus_32() +
00436 (init_table[2*i+1] >> 15) * twoToMinus_48();
00437
00438
00439
00440
00441
00442 }
00443
00444 carry = 0.0;
00445 if ( randoms[11] == 0. ) carry = twoToMinus_48();
00446 index = 11;
00447
00448 }
00449
00450 void Ranlux64Engine::setSeeds(const long * seeds, int lux) {
00451
00452
00453
00454
00455
00456
00457
00458 const int ecuyer_a = 53668;
00459 const int ecuyer_b = 40014;
00460 const int ecuyer_c = 12211;
00461 const int ecuyer_d = 2147483563;
00462
00463 const int lux_levels[3] = {109, 202, 397};
00464 const long *seedptr;
00465
00466 theSeeds = seeds;
00467 seedptr = seeds;
00468
00469 if(seeds == 0){
00470 setSeed(theSeed,lux);
00471 theSeeds = &theSeed;
00472 return;
00473 }
00474
00475 theSeed = *seeds;
00476
00477
00478
00479
00480 if( (lux > 2)||(lux < 0) ){
00481 pDiscard = (lux >= 12) ? (lux-12) : lux_levels[1];
00482 }else{
00483 pDiscard = lux_levels[luxury];
00484 }
00485 pDozens = pDiscard / 12;
00486 endIters = pDiscard % 12;
00487
00488 long init_table[24];
00489 long next_seed = *seeds;
00490 long k_multiple;
00491 int i;
00492
00493 for( i = 0;(i != 24)&&(*seedptr != 0);i++){
00494 init_table[i] = *seedptr & 0xffffffff;
00495 seedptr++;
00496 }
00497
00498 if(i != 24){
00499 next_seed = init_table[i-1];
00500 for(;i != 24;i++){
00501 k_multiple = next_seed / ecuyer_a;
00502 next_seed = ecuyer_b * (next_seed - k_multiple * ecuyer_a)
00503 - k_multiple * ecuyer_c;
00504 if(next_seed < 0) {
00505 next_seed += ecuyer_d;
00506 }
00507 next_seed &= 0xffffffff;
00508 init_table[i] = next_seed;
00509 }
00510 }
00511
00512 for(i = 0;i < 12; i++){
00513 randoms[i] = (init_table[2*i ] ) * 2.0 * twoToMinus_32() +
00514 (init_table[2*i+1] >> 15) * twoToMinus_48();
00515 }
00516
00517 carry = 0.0;
00518 if ( randoms[11] == 0. ) carry = twoToMinus_48();
00519 index = 11;
00520
00521 }
00522
00523 void Ranlux64Engine::saveStatus( const char filename[] ) const
00524 {
00525 std::ofstream outFile( filename, std::ios::out ) ;
00526 if (!outFile.bad()) {
00527 outFile << "Uvec\n";
00528 std::vector<unsigned long> v = put();
00529 for (unsigned int i=0; i<v.size(); ++i) {
00530 outFile << v[i] << "\n";
00531 }
00532 }
00533 }
00534
00535 void Ranlux64Engine::restoreStatus( const char filename[] )
00536 {
00537 std::ifstream inFile( filename, std::ios::in);
00538 if (!checkFile ( inFile, filename, engineName(), "restoreStatus" )) {
00539 std::cerr << " -- Engine state remains unchanged\n";
00540 return;
00541 }
00542 if ( possibleKeywordInput ( inFile, "Uvec", theSeed ) ) {
00543 std::vector<unsigned long> v;
00544 unsigned long xin;
00545 for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) {
00546 inFile >> xin;
00547 if (!inFile) {
00548 inFile.clear(std::ios::badbit | inFile.rdstate());
00549 std::cerr << "\nJamesRandom state (vector) description improper."
00550 << "\nrestoreStatus has failed."
00551 << "\nInput stream is probably mispositioned now." << std::endl;
00552 return;
00553 }
00554 v.push_back(xin);
00555 }
00556 getState(v);
00557 return;
00558 }
00559
00560 if (!inFile.bad() && !inFile.eof()) {
00561
00562 for (int i=0; i<12; ++i) {
00563 inFile >> randoms[i];
00564 }
00565 inFile >> carry; inFile >> index;
00566 inFile >> luxury; inFile >> pDiscard;
00567 pDozens = pDiscard / 12;
00568 endIters = pDiscard % 12;
00569 }
00570 }
00571
00572 void Ranlux64Engine::showStatus() const
00573 {
00574 std::cout << std::endl;
00575 std::cout << "--------- Ranlux engine status ---------" << std::endl;
00576 std::cout << " Initial seed = " << theSeed << std::endl;
00577 std::cout << " randoms[] = ";
00578 for (int i=0; i<12; ++i) {
00579 std::cout << randoms[i] << std::endl;
00580 }
00581 std::cout << std::endl;
00582 std::cout << " carry = " << carry << ", index = " << index << std::endl;
00583 std::cout << " luxury = " << luxury << " pDiscard = "
00584 << pDiscard << std::endl;
00585 std::cout << "----------------------------------------" << std::endl;
00586 }
00587
00588 std::ostream & Ranlux64Engine::put( std::ostream& os ) const
00589 {
00590 char beginMarker[] = "Ranlux64Engine-begin";
00591 os << beginMarker << "\nUvec\n";
00592 std::vector<unsigned long> v = put();
00593 for (unsigned int i=0; i<v.size(); ++i) {
00594 os << v[i] << "\n";
00595 }
00596 return os;
00597 }
00598
00599 std::vector<unsigned long> Ranlux64Engine::put () const {
00600 std::vector<unsigned long> v;
00601 v.push_back (engineIDulong<Ranlux64Engine>());
00602 std::vector<unsigned long> t;
00603 for (int i=0; i<12; ++i) {
00604 t = DoubConv::dto2longs(randoms[i]);
00605 v.push_back(t[0]); v.push_back(t[1]);
00606 }
00607 t = DoubConv::dto2longs(carry);
00608 v.push_back(t[0]); v.push_back(t[1]);
00609 v.push_back(static_cast<unsigned long>(index));
00610 v.push_back(static_cast<unsigned long>(luxury));
00611 v.push_back(static_cast<unsigned long>(pDiscard));
00612 return v;
00613 }
00614
00615 std::istream & Ranlux64Engine::get ( std::istream& is )
00616 {
00617 char beginMarker [MarkerLen];
00618 is >> std::ws;
00619 is.width(MarkerLen);
00620
00621
00622 is >> beginMarker;
00623 if (strcmp(beginMarker,"Ranlux64Engine-begin")) {
00624 is.clear(std::ios::badbit | is.rdstate());
00625 std::cerr << "\nInput stream mispositioned or"
00626 << "\nRanlux64Engine state description missing or"
00627 << "\nwrong engine type found." << std::endl;
00628 return is;
00629 }
00630 return getState(is);
00631 }
00632
00633 std::string Ranlux64Engine::beginTag ( ) {
00634 return "Ranlux64Engine-begin";
00635 }
00636
00637 std::istream & Ranlux64Engine::getState ( std::istream& is )
00638 {
00639 if ( possibleKeywordInput ( is, "Uvec", theSeed ) ) {
00640 std::vector<unsigned long> v;
00641 unsigned long uu;
00642 for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) {
00643 is >> uu;
00644 if (!is) {
00645 is.clear(std::ios::badbit | is.rdstate());
00646 std::cerr << "\nRanlux64Engine state (vector) description improper."
00647 << "\ngetState() has failed."
00648 << "\nInput stream is probably mispositioned now." << std::endl;
00649 return is;
00650 }
00651 v.push_back(uu);
00652 }
00653 getState(v);
00654 return (is);
00655 }
00656
00657
00658
00659 char endMarker [MarkerLen];
00660 for (int i=0; i<12; ++i) {
00661 is >> randoms[i];
00662 }
00663 is >> carry; is >> index;
00664 is >> luxury; is >> pDiscard;
00665 pDozens = pDiscard / 12;
00666 endIters = pDiscard % 12;
00667 is >> std::ws;
00668 is.width(MarkerLen);
00669 is >> endMarker;
00670 if (strcmp(endMarker,"Ranlux64Engine-end")) {
00671 is.clear(std::ios::badbit | is.rdstate());
00672 std::cerr << "\nRanlux64Engine state description incomplete."
00673 << "\nInput stream is probably mispositioned now." << std::endl;
00674 return is;
00675 }
00676 return is;
00677 }
00678
00679 bool Ranlux64Engine::get (const std::vector<unsigned long> & v) {
00680 if ((v[0] & 0xffffffffUL) != engineIDulong<Ranlux64Engine>()) {
00681 std::cerr <<
00682 "\nRanlux64Engine get:state vector has wrong ID word - state unchanged\n";
00683 return false;
00684 }
00685 return getState(v);
00686 }
00687
00688 bool Ranlux64Engine::getState (const std::vector<unsigned long> & v) {
00689 if (v.size() != VECTOR_STATE_SIZE ) {
00690 std::cerr <<
00691 "\nRanlux64Engine get:state vector has wrong length - state unchanged\n";
00692 return false;
00693 }
00694 std::vector<unsigned long> t(2);
00695 for (int i=0; i<12; ++i) {
00696 t[0] = v[2*i+1]; t[1] = v[2*i+2];
00697 randoms[i] = DoubConv::longs2double(t);
00698 }
00699 t[0] = v[25]; t[1] = v[26];
00700 carry = DoubConv::longs2double(t);
00701 index = v[27];
00702 luxury = v[28];
00703 pDiscard = v[29];
00704 return true;
00705 }
00706
00707 }