56 #include "CLHEP/Random/Random.h"
57 #include "CLHEP/Random/Ranlux64Engine.h"
58 #include "CLHEP/Random/engineIDulong.h"
59 #include "CLHEP/Random/DoubConv.h"
60 #include "CLHEP/Utility/atomic_int.h"
70 CLHEP_ATOMIC_INT_TYPE numberOfEngines(0);
73 const int maxIndex = 215;
82 template< std::size_t
n,
83 bool = n < std::size_t(std::numeric_limits<unsigned long>::digits) >
84 struct do_right_shift;
85 template< std::
size_t n >
86 struct do_right_shift<n,true>
88 unsigned long operator()(
unsigned long value) {
return value >>
n; }
90 template< std::
size_t n >
91 struct do_right_shift<n,false>
96 template< std::
size_t nbits >
97 unsigned long rshift(
unsigned long value )
98 {
return do_right_shift<nbits>()(value); }
105 Ranlux64Engine::Ranlux64Engine()
109 int numEngines = numberOfEngines++;
110 int cycle = std::abs(
int(numEngines/maxIndex));
111 int curIndex = std::abs(
int(numEngines%maxIndex));
113 long mask = ((cycle & 0x007fffff) << 8);
115 HepRandom::getTheTableSeeds( seedlist, curIndex );
125 Ranlux64Engine::Ranlux64Engine(
long seed,
int lux)
129 long seedlist[2]={seed,0};
131 advance ( 2*lux + 1 );
135 Ranlux64Engine::Ranlux64Engine(
int rowIndex,
int,
int lux)
139 int cycle = std::abs(
int(rowIndex/maxIndex));
140 int row = std::abs(
int(rowIndex%maxIndex));
141 long mask = (( cycle & 0x000007ff ) << 20 );
143 HepRandom::getTheTableSeeds( seedlist, row );
149 Ranlux64Engine::Ranlux64Engine( std::istream& is )
155 Ranlux64Engine::~Ranlux64Engine() {}
163 if (index <= 0) update();
164 return randoms[--index] + twoToMinus_49();
167 void Ranlux64Engine::update() {
196 if ( endIters == 1 ) {
197 y1 = randoms[ 4] - randoms[11] - carry;
200 carry = twoToMinus_48();
204 randoms[11] = randoms[10];
205 randoms[10] = randoms[ 9];
206 randoms[ 9] = randoms[ 8];
207 randoms[ 8] = randoms[ 7];
208 randoms[ 7] = randoms[ 6];
209 randoms[ 6] = randoms[ 5];
210 randoms[ 5] = randoms[ 4];
211 randoms[ 4] = randoms[ 3];
212 randoms[ 3] = randoms[ 2];
213 randoms[ 2] = randoms[ 1];
214 randoms[ 1] = randoms[ 0];
220 for ( m = 0, nr = 11, ns = 4; m < endIters; ++
m, --nr ) {
221 y1 = randoms [
ns] - randoms[nr] - carry;
224 carry = twoToMinus_48();
236 for (m=0; m<12; m++) {
241 for (m=11; m>=0; --
m) {
242 randoms[
m] = temp[
ns];
257 void Ranlux64Engine::advance(
int dozens) {
260 double cValue = twoToMinus_48();
283 for ( k = dozens; k > 0; --k ) {
285 y1 = randoms[ 4] - randoms[11] - carry;
287 y2 = randoms[ 3] - randoms[10];
294 y3 = randoms[ 2] - randoms[ 9];
301 y1 = randoms[ 1] - randoms[ 8];
308 y2 = randoms[ 0] - randoms[ 7];
315 y3 = randoms[11] - randoms[ 6];
322 y1 = randoms[10] - randoms[ 5];
329 y2 = randoms[ 9] - randoms[ 4];
336 y3 = randoms[ 8] - randoms[ 3];
343 y1 = randoms[ 7] - randoms[ 2];
350 y2 = randoms[ 6] - randoms[ 1];
357 y3 = randoms[ 5] - randoms[ 0];
374 void Ranlux64Engine::flatArray(
const int size,
double* vect) {
375 for(
int i=0; i < size; ++i ) {
380 void Ranlux64Engine::setSeed(
long seed,
int lux) {
388 const int ecuyer_a(53668);
389 const int ecuyer_b(40014);
390 const int ecuyer_c(12211);
391 const int ecuyer_d(2147483563);
393 const int lux_levels[3] = {109, 202, 397};
396 if( (lux > 2)||(lux < 0) ){
397 pDiscard = (lux >= 12) ? (lux-12) : lux_levels[1];
399 pDiscard = lux_levels[luxury];
401 pDozens = pDiscard / 12;
402 endIters = pDiscard % 12;
405 long next_seed = seed;
408 next_seed &= 0xffffffff;
409 while( next_seed >= ecuyer_d ) {
410 next_seed -= ecuyer_d;
413 for(i = 0;i != 24;i++){
414 k_multiple = next_seed / ecuyer_a;
415 next_seed = ecuyer_b * (next_seed - k_multiple * ecuyer_a)
416 - k_multiple * ecuyer_c;
418 next_seed += ecuyer_d;
420 next_seed &= 0xffffffff;
421 init_table[i] = next_seed;
424 if(
sizeof(
long) >= 8 ) {
425 long topbits1, topbits2;
427 topbits1 = ( seed >> 32) & 0xffff ;
428 topbits2 = ( seed >> 48) & 0xffff ;
430 topbits1 = detail::rshift<32>(seed) & 0xffff ;
431 topbits2 = detail::rshift<48>(seed) & 0xffff ;
433 init_table[0] ^= topbits1;
434 init_table[2] ^= topbits2;
439 for(i = 0;i < 12; i++){
440 randoms[i] = (init_table[2*i ] ) * 2.0 * twoToMinus_32() +
441 (init_table[2*i+1] >> 15) * twoToMinus_48();
450 if ( randoms[11] == 0. ) carry = twoToMinus_48();
463 const int ecuyer_a = 53668;
464 const int ecuyer_b = 40014;
465 const int ecuyer_c = 12211;
466 const int ecuyer_d = 2147483563;
468 const int lux_levels[3] = {109, 202, 397};
475 setSeed(theSeed,lux);
485 if( (lux > 2)||(lux < 0) ){
486 pDiscard = (lux >= 12) ? (lux-12) : lux_levels[1];
488 pDiscard = lux_levels[luxury];
490 pDozens = pDiscard / 12;
491 endIters = pDiscard % 12;
494 long next_seed = *seeds;
498 for( i = 0;(i != 24)&&(*seedptr != 0);i++){
499 init_table[i] = *seedptr & 0xffffffff;
504 next_seed = init_table[i-1];
506 k_multiple = next_seed / ecuyer_a;
507 next_seed = ecuyer_b * (next_seed - k_multiple * ecuyer_a)
508 - k_multiple * ecuyer_c;
510 next_seed += ecuyer_d;
512 next_seed &= 0xffffffff;
513 init_table[i] = next_seed;
517 for(i = 0;i < 12; i++){
518 randoms[i] = (init_table[2*i ] ) * 2.0 * twoToMinus_32() +
519 (init_table[2*i+1] >> 15) * twoToMinus_48();
523 if ( randoms[11] == 0. ) carry = twoToMinus_48();
528 void Ranlux64Engine::saveStatus(
const char filename[] )
const
530 std::ofstream outFile( filename, std::ios::out ) ;
531 if (!outFile.bad()) {
533 std::vector<unsigned long> v = put();
534 for (
unsigned int i=0; i<v.size(); ++i) {
535 outFile << v[i] <<
"\n";
540 void Ranlux64Engine::restoreStatus(
const char filename[] )
542 std::ifstream inFile( filename, std::ios::in);
543 if (!checkFile ( inFile, filename, engineName(),
"restoreStatus" )) {
544 std::cerr <<
" -- Engine state remains unchanged\n";
547 if ( possibleKeywordInput ( inFile,
"Uvec", theSeed ) ) {
548 std::vector<unsigned long> v;
550 for (
unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) {
553 inFile.clear(std::ios::badbit | inFile.rdstate());
554 std::cerr <<
"\nJamesRandom state (vector) description improper."
555 <<
"\nrestoreStatus has failed."
556 <<
"\nInput stream is probably mispositioned now." << std::endl;
565 if (!inFile.bad() && !inFile.eof()) {
567 for (
int i=0; i<12; ++i) {
568 inFile >> randoms[i];
570 inFile >> carry; inFile >> index;
571 inFile >> luxury; inFile >> pDiscard;
572 pDozens = pDiscard / 12;
573 endIters = pDiscard % 12;
577 void Ranlux64Engine::showStatus()
const
579 std::cout << std::endl;
580 std::cout <<
"--------- Ranlux engine status ---------" << std::endl;
581 std::cout <<
" Initial seed = " << theSeed << std::endl;
582 std::cout <<
" randoms[] = ";
583 for (
int i=0; i<12; ++i) {
584 std::cout << randoms[i] << std::endl;
586 std::cout << std::endl;
587 std::cout <<
" carry = " << carry <<
", index = " << index << std::endl;
588 std::cout <<
" luxury = " << luxury <<
" pDiscard = "
589 << pDiscard << std::endl;
590 std::cout <<
"----------------------------------------" << std::endl;
593 std::ostream & Ranlux64Engine::put( std::ostream& os )
const
595 char beginMarker[] =
"Ranlux64Engine-begin";
596 os << beginMarker <<
"\nUvec\n";
597 std::vector<unsigned long> v = put();
598 for (
unsigned int i=0; i<v.size(); ++i) {
604 std::vector<unsigned long> Ranlux64Engine::put ()
const {
605 std::vector<unsigned long> v;
606 v.push_back (engineIDulong<Ranlux64Engine>());
607 std::vector<unsigned long> t;
608 for (
int i=0; i<12; ++i) {
609 t = DoubConv::dto2longs(randoms[i]);
610 v.push_back(t[0]); v.push_back(t[1]);
612 t = DoubConv::dto2longs(carry);
613 v.push_back(t[0]); v.push_back(t[1]);
614 v.push_back(static_cast<unsigned long>(index));
615 v.push_back(static_cast<unsigned long>(luxury));
616 v.push_back(static_cast<unsigned long>(pDiscard));
620 std::istream & Ranlux64Engine::get ( std::istream& is )
628 if (strcmp(beginMarker,
"Ranlux64Engine-begin")) {
629 is.clear(std::ios::badbit | is.rdstate());
630 std::cerr <<
"\nInput stream mispositioned or"
631 <<
"\nRanlux64Engine state description missing or"
632 <<
"\nwrong engine type found." << std::endl;
638 std::string Ranlux64Engine::beginTag ( ) {
639 return "Ranlux64Engine-begin";
642 std::istream & Ranlux64Engine::getState ( std::istream& is )
644 if ( possibleKeywordInput ( is,
"Uvec", theSeed ) ) {
645 std::vector<unsigned long> v;
647 for (
unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) {
650 is.clear(std::ios::badbit | is.rdstate());
651 std::cerr <<
"\nRanlux64Engine state (vector) description improper."
652 <<
"\ngetState() has failed."
653 <<
"\nInput stream is probably mispositioned now." << std::endl;
665 for (
int i=0; i<12; ++i) {
668 is >> carry; is >> index;
669 is >> luxury; is >> pDiscard;
670 pDozens = pDiscard / 12;
671 endIters = pDiscard % 12;
675 if (strcmp(endMarker,
"Ranlux64Engine-end")) {
676 is.clear(std::ios::badbit | is.rdstate());
677 std::cerr <<
"\nRanlux64Engine state description incomplete."
678 <<
"\nInput stream is probably mispositioned now." << std::endl;
684 bool Ranlux64Engine::get (
const std::vector<unsigned long> & v) {
685 if ((v[0] & 0xffffffffUL) != engineIDulong<Ranlux64Engine>()) {
687 "\nRanlux64Engine get:state vector has wrong ID word - state unchanged\n";
693 bool Ranlux64Engine::getState (
const std::vector<unsigned long> & v) {
694 if (v.size() != VECTOR_STATE_SIZE ) {
696 "\nRanlux64Engine get:state vector has wrong length - state unchanged\n";
699 std::vector<unsigned long> t(2);
700 for (
int i=0; i<12; ++i) {
701 t[0] = v[2*i+1]; t[1] = v[2*i+2];
702 randoms[i] = DoubConv::longs2double(t);
704 t[0] = v[25]; t[1] = v[26];
705 carry = DoubConv::longs2double(t);
unsigned long rshift(unsigned long value)
static const int MarkerLen
unsigned long operator()(unsigned long)
void setSeeds(const SeedVector &sv)
Set the seeds of the current generator.