56 #include "CLHEP/Random/defs.h"
57 #include "CLHEP/Random/Random.h"
58 #include "CLHEP/Random/Ranlux64Engine.h"
59 #include "CLHEP/Random/engineIDulong.h"
60 #include "CLHEP/Random/DoubConv.hh"
67 static const int MarkerLen = 64;
71 int Ranlux64Engine::numEngines = 0;
74 int Ranlux64Engine::maxIndex = 215;
79 template< std::size_t
n,
80 bool = n < std::size_t(std::numeric_limits<unsigned long>::digits) >
81 struct do_right_shift;
82 template< std::
size_t n >
83 struct do_right_shift<n,true>
85 unsigned long operator()(
unsigned long value) {
return value >>
n; }
87 template< std::
size_t n >
88 struct do_right_shift<
n,false>
93 template< std::
size_t nbits >
94 unsigned long rshift(
unsigned long value )
95 {
return do_right_shift<nbits>()(value); }
106 int cycle = std::abs(
int(numEngines/maxIndex));
107 int curIndex = std::abs(
int(numEngines%maxIndex));
109 long mask = ((cycle & 0x007fffff) << 8);
125 long seedlist[2]={seed,0};
127 advance ( 2*lux + 1 );
135 int cycle = std::abs(
int(rowIndex/maxIndex));
136 int row = std::abs(
int(rowIndex%maxIndex));
137 long mask = (( cycle & 0x000007ff ) << 20 );
159 if (index <= 0) update();
163 void Ranlux64Engine::update() {
192 if ( endIters == 1 ) {
193 y1 = randoms[ 4] - randoms[11] - carry;
200 randoms[11] = randoms[10];
201 randoms[10] = randoms[ 9];
202 randoms[ 9] = randoms[ 8];
203 randoms[ 8] = randoms[ 7];
204 randoms[ 7] = randoms[ 6];
205 randoms[ 6] = randoms[ 5];
206 randoms[ 5] = randoms[ 4];
207 randoms[ 4] = randoms[ 3];
208 randoms[ 3] = randoms[ 2];
209 randoms[ 2] = randoms[ 1];
210 randoms[ 1] = randoms[ 0];
216 for ( m = 0, nr = 11, ns = 4; m < endIters; ++m, --nr ) {
217 y1 = randoms [ns] - randoms[nr] - carry;
232 for (m=0; m<12; m++) {
237 for (m=11; m>=0; --m) {
238 randoms[m] = temp[ns];
253 void Ranlux64Engine::advance(
int dozens) {
255 register double y1, y2, y3;
257 register double zero = 0.0;
258 register double one = 1.0;
279 for ( k = dozens; k > 0; --k ) {
281 y1 = randoms[ 4] - randoms[11] - carry;
283 y2 = randoms[ 3] - randoms[10];
290 y3 = randoms[ 2] - randoms[ 9];
297 y1 = randoms[ 1] - randoms[ 8];
304 y2 = randoms[ 0] - randoms[ 7];
311 y3 = randoms[11] - randoms[ 6];
318 y1 = randoms[10] - randoms[ 5];
325 y2 = randoms[ 9] - randoms[ 4];
332 y3 = randoms[ 8] - randoms[ 3];
339 y1 = randoms[ 7] - randoms[ 2];
346 y2 = randoms[ 6] - randoms[ 1];
353 y3 = randoms[ 5] - randoms[ 0];
371 for(
int i=0; i < size; ++i ) {
384 const int ecuyer_a(53668);
385 const int ecuyer_b(40014);
386 const int ecuyer_c(12211);
387 const int ecuyer_d(2147483563);
389 const int lux_levels[3] = {109, 202, 397};
392 if( (lux > 2)||(lux < 0) ){
393 pDiscard = (lux >= 12) ? (lux-12) : lux_levels[1];
395 pDiscard = lux_levels[luxury];
397 pDozens = pDiscard / 12;
398 endIters = pDiscard % 12;
401 long next_seed = seed;
404 next_seed &= 0xffffffff;
405 while( next_seed >= ecuyer_d ) {
406 next_seed -= ecuyer_d;
409 for(i = 0;i != 24;i++){
410 k_multiple = next_seed / ecuyer_a;
411 next_seed = ecuyer_b * (next_seed - k_multiple * ecuyer_a)
412 - k_multiple * ecuyer_c;
414 next_seed += ecuyer_d;
416 next_seed &= 0xffffffff;
417 init_table[i] = next_seed;
420 if(
sizeof(
long) >= 8 ) {
421 long topbits1, topbits2;
423 topbits1 = ( seed >> 32) & 0xffff ;
424 topbits2 = ( seed >> 48) & 0xffff ;
426 topbits1 = detail::rshift<32>(seed) & 0xffff ;
427 topbits2 = detail::rshift<48>(seed) & 0xffff ;
429 init_table[0] ^= topbits1;
430 init_table[2] ^= topbits2;
435 for(i = 0;i < 12; i++){
459 const int ecuyer_a = 53668;
460 const int ecuyer_b = 40014;
461 const int ecuyer_c = 12211;
462 const int ecuyer_d = 2147483563;
464 const int lux_levels[3] = {109, 202, 397};
481 if( (lux > 2)||(lux < 0) ){
482 pDiscard = (lux >= 12) ? (lux-12) : lux_levels[1];
484 pDiscard = lux_levels[luxury];
486 pDozens = pDiscard / 12;
487 endIters = pDiscard % 12;
490 long next_seed = *seeds;
494 for( i = 0;(i != 24)&&(*seedptr != 0);i++){
495 init_table[i] = *seedptr & 0xffffffff;
500 next_seed = init_table[i-1];
502 k_multiple = next_seed / ecuyer_a;
503 next_seed = ecuyer_b * (next_seed - k_multiple * ecuyer_a)
504 - k_multiple * ecuyer_c;
506 next_seed += ecuyer_d;
508 next_seed &= 0xffffffff;
509 init_table[i] = next_seed;
513 for(i = 0;i < 12; i++){
526 std::ofstream outFile( filename, std::ios::out ) ;
527 if (!outFile.bad()) {
529 std::vector<unsigned long> v =
put();
531 std::cout <<
"Result of v = put() is:\n";
533 for (
unsigned int i=0; i<v.size(); ++i) {
534 outFile << v[i] <<
"\n";
536 std::cout << v[i] <<
" ";
537 if (i%6==0) std::cout <<
"\n";
545 if (!outFile.bad()) {
546 outFile <<
theSeed << std::endl;
547 for (
int i=0; i<12; ++i) {
548 outFile <<std::setprecision(20) << randoms[i] << std::endl;
550 outFile << std::setprecision(20) << carry <<
" " << index << std::endl;
551 outFile << luxury <<
" " << pDiscard << std::endl;
558 std::ifstream inFile( filename, std::ios::in);
560 std::cerr <<
" -- Engine state remains unchanged\n";
564 std::vector<unsigned long> v;
569 std::cout <<
"ivec = " << ivec <<
" xin = " << xin <<
" ";
570 if (ivec%3 == 0) std::cout <<
"\n";
573 inFile.clear(std::ios::badbit | inFile.rdstate());
574 std::cerr <<
"\nJamesRandom state (vector) description improper."
575 <<
"\nrestoreStatus has failed."
576 <<
"\nInput stream is probably mispositioned now." << std::endl;
585 if (!inFile.bad() && !inFile.eof()) {
587 for (
int i=0; i<12; ++i) {
588 inFile >> randoms[i];
590 inFile >> carry; inFile >> index;
591 inFile >> luxury; inFile >> pDiscard;
592 pDozens = pDiscard / 12;
593 endIters = pDiscard % 12;
599 std::cout << std::endl;
600 std::cout <<
"--------- Ranlux engine status ---------" << std::endl;
601 std::cout <<
" Initial seed = " <<
theSeed << std::endl;
602 std::cout <<
" randoms[] = ";
603 for (
int i=0; i<12; ++i) {
604 std::cout << randoms[i] << std::endl;
606 std::cout << std::endl;
607 std::cout <<
" carry = " << carry <<
", index = " << index << std::endl;
608 std::cout <<
" luxury = " << luxury <<
" pDiscard = "
609 << pDiscard << std::endl;
610 std::cout <<
"----------------------------------------" << std::endl;
615 char beginMarker[] =
"Ranlux64Engine-begin";
616 os << beginMarker <<
"\nUvec\n";
617 std::vector<unsigned long> v =
put();
618 for (
unsigned int i=0; i<v.size(); ++i) {
623 char endMarker[] =
"Ranlux64Engine-end";
624 int pr = os.precision(20);
625 os <<
" " << beginMarker <<
" ";
627 for (
int i=0; i<12; ++i) {
628 os << randoms[i] << std::endl;
630 os << carry <<
" " << index <<
" ";
631 os << luxury <<
" " << pDiscard <<
"\n";
632 os << endMarker <<
" ";
639 std::vector<unsigned long> v;
640 v.push_back (engineIDulong<Ranlux64Engine>());
641 std::vector<unsigned long> t;
642 for (
int i=0; i<12; ++i) {
644 v.push_back(t[0]); v.push_back(t[1]);
647 v.push_back(t[0]); v.push_back(t[1]);
648 v.push_back(static_cast<unsigned long>(index));
649 v.push_back(static_cast<unsigned long>(luxury));
650 v.push_back(static_cast<unsigned long>(pDiscard));
656 char beginMarker [MarkerLen];
662 if (strcmp(beginMarker,
"Ranlux64Engine-begin")) {
663 is.clear(std::ios::badbit | is.rdstate());
664 std::cerr <<
"\nInput stream mispositioned or"
665 <<
"\nRanlux64Engine state description missing or"
666 <<
"\nwrong engine type found." << std::endl;
673 return "Ranlux64Engine-begin";
679 std::vector<unsigned long> v;
684 is.clear(std::ios::badbit | is.rdstate());
685 std::cerr <<
"\nRanlux64Engine state (vector) description improper."
686 <<
"\ngetState() has failed."
687 <<
"\nInput stream is probably mispositioned now." << std::endl;
698 char endMarker [MarkerLen];
699 for (
int i=0; i<12; ++i) {
702 is >> carry; is >> index;
703 is >> luxury; is >> pDiscard;
704 pDozens = pDiscard / 12;
705 endIters = pDiscard % 12;
709 if (strcmp(endMarker,
"Ranlux64Engine-end")) {
710 is.clear(std::ios::badbit | is.rdstate());
711 std::cerr <<
"\nRanlux64Engine state description incomplete."
712 <<
"\nInput stream is probably mispositioned now." << std::endl;
719 if ((v[0] & 0xffffffffUL) != engineIDulong<Ranlux64Engine>()) {
721 "\nRanlux64Engine get:state vector has wrong ID word - state unchanged\n";
730 "\nRanlux64Engine get:state vector has wrong length - state unchanged\n";
733 std::vector<unsigned long> t(2);
734 for (
int i=0; i<12; ++i) {
735 t[0] = v[2*i+1]; t[1] = v[2*i+2];
738 t[0] = v[25]; t[1] = v[26];