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.h"
61#include "CLHEP/Utility/atomic_int.h"
77 const int maxIndex = 215;
80static const int MarkerLen = 64;
86template< std::size_t
n,
87 bool =
n < std::size_t(std::numeric_limits<unsigned long>::digits) >
88 struct do_right_shift;
89template< std::
size_t n >
90 struct do_right_shift<n,true>
92 unsigned long operator()(
unsigned long value) {
return value >>
n; }
94template< std::
size_t n >
95 struct do_right_shift<
n,false>
100template< std::
size_t nbits >
101 unsigned long rshift(
unsigned long value )
102{
return do_right_shift<nbits>()(value); }
113 int numEngines = numberOfEngines++;
114 int cycle = std::abs(
int(numEngines/maxIndex));
115 int curIndex = std::abs(
int(numEngines%maxIndex));
117 long mask = ((cycle & 0x007fffff) << 8);
133 long seedlist[2]={seed,0};
135 advance ( 2*lux + 1 );
143 int cycle = std::abs(
int(rowIndex/maxIndex));
144 int row = std::abs(
int(rowIndex%maxIndex));
145 long mask = (( cycle & 0x000007ff ) << 20 );
167 if (index <= 0) update();
171void Ranlux64Engine::update() {
200 if ( endIters == 1 ) {
201 y1 = randoms[ 4] - randoms[11] - carry;
208 randoms[11] = randoms[10];
209 randoms[10] = randoms[ 9];
210 randoms[ 9] = randoms[ 8];
211 randoms[ 8] = randoms[ 7];
212 randoms[ 7] = randoms[ 6];
213 randoms[ 6] = randoms[ 5];
214 randoms[ 5] = randoms[ 4];
215 randoms[ 4] = randoms[ 3];
216 randoms[ 3] = randoms[ 2];
217 randoms[ 2] = randoms[ 1];
218 randoms[ 1] = randoms[ 0];
224 for ( m = 0, nr = 11, ns = 4; m < endIters; ++m, --nr ) {
225 y1 = randoms [ns] - randoms[nr] - carry;
240 for (m=0; m<12; m++) {
245 for (m=11; m>=0; --m) {
246 randoms[m] = temp[ns];
261void Ranlux64Engine::advance(
int dozens) {
287 for ( k = dozens; k > 0; --k ) {
289 y1 = randoms[ 4] - randoms[11] - carry;
291 y2 = randoms[ 3] - randoms[10];
298 y3 = randoms[ 2] - randoms[ 9];
305 y1 = randoms[ 1] - randoms[ 8];
312 y2 = randoms[ 0] - randoms[ 7];
319 y3 = randoms[11] - randoms[ 6];
326 y1 = randoms[10] - randoms[ 5];
333 y2 = randoms[ 9] - randoms[ 4];
340 y3 = randoms[ 8] - randoms[ 3];
347 y1 = randoms[ 7] - randoms[ 2];
354 y2 = randoms[ 6] - randoms[ 1];
361 y3 = randoms[ 5] - randoms[ 0];
379 for(
int i=0; i < size; ++i ) {
392 const int ecuyer_a(53668);
393 const int ecuyer_b(40014);
394 const int ecuyer_c(12211);
395 const int ecuyer_d(2147483563);
397 const int lux_levels[3] = {109, 202, 397};
400 if( (lux > 2)||(lux < 0) ){
401 pDiscard = (lux >= 12) ? (lux-12) : lux_levels[1];
403 pDiscard = lux_levels[luxury];
405 pDozens = pDiscard / 12;
406 endIters = pDiscard % 12;
409 long next_seed = seed;
412 next_seed &= 0xffffffff;
413 while( next_seed >= ecuyer_d ) {
414 next_seed -= ecuyer_d;
417 for(i = 0;i != 24;i++){
418 k_multiple = next_seed / ecuyer_a;
419 next_seed = ecuyer_b * (next_seed - k_multiple * ecuyer_a)
420 - k_multiple * ecuyer_c;
422 next_seed += ecuyer_d;
424 next_seed &= 0xffffffff;
425 init_table[i] = next_seed;
428 if(
sizeof(
long) >= 8 ) {
429 int64_t topbits1, topbits2;
431 topbits1 = ( (int64_t) seed >> 32) & 0xffff ;
432 topbits2 = ( (int64_t) seed >> 48) & 0xffff ;
434 topbits1 = detail::rshift<32>(seed) & 0xffff ;
435 topbits2 = detail::rshift<48>(seed) & 0xffff ;
437 init_table[0] ^= topbits1;
438 init_table[2] ^= topbits2;
443 for(i = 0;i < 12; i++){
468 const int ecuyer_a = 53668;
469 const int ecuyer_b = 40014;
470 const int ecuyer_c = 12211;
471 const int ecuyer_d = 2147483563;
473 const int lux_levels[3] = {109, 202, 397};
490 if( (lux > 2)||(lux < 0) ){
491 pDiscard = (lux >= 12) ? (lux-12) : lux_levels[1];
493 pDiscard = lux_levels[luxury];
495 pDozens = pDiscard / 12;
496 endIters = pDiscard % 12;
499 long next_seed = *seeds;
503 for( i = 0;(i != 24)&&(*seedptr != 0);i++){
504 init_table[i] = *seedptr & 0xffffffff;
509 next_seed = init_table[i-1];
511 k_multiple = next_seed / ecuyer_a;
512 next_seed = ecuyer_b * (next_seed - k_multiple * ecuyer_a)
513 - k_multiple * ecuyer_c;
515 next_seed += ecuyer_d;
517 next_seed &= 0xffffffff;
518 init_table[i] = next_seed;
522 for(i = 0;i < 12; i++){
536 std::ofstream outFile( filename, std::ios::out ) ;
537 if (!outFile.bad()) {
539 std::vector<unsigned long> v =
put();
541 std::cout <<
"Result of v = put() is:\n";
543 for (
unsigned int i=0; i<v.size(); ++i) {
544 outFile << v[i] <<
"\n";
546 std::cout << v[i] <<
" ";
547 if (i%6==0) std::cout <<
"\n";
555 if (!outFile.bad()) {
556 outFile <<
theSeed << std::endl;
557 for (
int i=0; i<12; ++i) {
558 outFile <<std::setprecision(20) << randoms[i] << std::endl;
560 outFile << std::setprecision(20) << carry <<
" " << index << std::endl;
561 outFile << luxury <<
" " << pDiscard << std::endl;
568 std::ifstream inFile( filename, std::ios::in);
570 std::cerr <<
" -- Engine state remains unchanged\n";
574 std::vector<unsigned long> v;
579 std::cout <<
"ivec = " << ivec <<
" xin = " << xin <<
" ";
580 if (ivec%3 == 0) std::cout <<
"\n";
583 inFile.clear(std::ios::badbit | inFile.rdstate());
584 std::cerr <<
"\nJamesRandom state (vector) description improper."
585 <<
"\nrestoreStatus has failed."
586 <<
"\nInput stream is probably mispositioned now." << std::endl;
595 if (!inFile.bad() && !inFile.eof()) {
597 for (
int i=0; i<12; ++i) {
598 inFile >> randoms[i];
600 inFile >> carry; inFile >> index;
601 inFile >> luxury; inFile >> pDiscard;
602 pDozens = pDiscard / 12;
603 endIters = pDiscard % 12;
609 std::cout << std::endl;
610 std::cout <<
"--------- Ranlux engine status ---------" << std::endl;
611 std::cout <<
" Initial seed = " <<
theSeed << std::endl;
612 std::cout <<
" randoms[] = ";
613 for (
int i=0; i<12; ++i) {
614 std::cout << randoms[i] << std::endl;
616 std::cout << std::endl;
617 std::cout <<
" carry = " << carry <<
", index = " << index << std::endl;
618 std::cout <<
" luxury = " << luxury <<
" pDiscard = "
619 << pDiscard << std::endl;
620 std::cout <<
"----------------------------------------" << std::endl;
625 char beginMarker[] =
"Ranlux64Engine-begin";
626 os << beginMarker <<
"\nUvec\n";
627 std::vector<unsigned long> v =
put();
628 for (
unsigned int i=0; i<v.size(); ++i) {
633 char endMarker[] =
"Ranlux64Engine-end";
634 long pr = os.precision(20);
635 os <<
" " << beginMarker <<
" ";
637 for (
int i=0; i<12; ++i) {
638 os << randoms[i] << std::endl;
640 os << carry <<
" " << index <<
" ";
641 os << luxury <<
" " << pDiscard <<
"\n";
642 os << endMarker <<
" ";
649 std::vector<unsigned long> v;
650 v.push_back (engineIDulong<Ranlux64Engine>());
651 std::vector<unsigned long> t;
652 for (
int i=0; i<12; ++i) {
654 v.push_back(t[0]); v.push_back(t[1]);
657 v.push_back(t[0]); v.push_back(t[1]);
658 v.push_back(
static_cast<unsigned long>(index));
659 v.push_back(
static_cast<unsigned long>(luxury));
660 v.push_back(
static_cast<unsigned long>(pDiscard));
666 char beginMarker [MarkerLen];
672 if (strcmp(beginMarker,
"Ranlux64Engine-begin")) {
673 is.clear(std::ios::badbit | is.rdstate());
674 std::cerr <<
"\nInput stream mispositioned or"
675 <<
"\nRanlux64Engine state description missing or"
676 <<
"\nwrong engine type found." << std::endl;
683 return "Ranlux64Engine-begin";
689 std::vector<unsigned long> v;
694 is.clear(std::ios::badbit | is.rdstate());
695 std::cerr <<
"\nRanlux64Engine state (vector) description improper."
696 <<
"\ngetState() has failed."
697 <<
"\nInput stream is probably mispositioned now." << std::endl;
708 char endMarker [MarkerLen];
709 for (
int i=0; i<12; ++i) {
712 is >> carry; is >> index;
713 is >> luxury; is >> pDiscard;
714 pDozens = pDiscard / 12;
715 endIters = pDiscard % 12;
719 if (strcmp(endMarker,
"Ranlux64Engine-end")) {
720 is.clear(std::ios::badbit | is.rdstate());
721 std::cerr <<
"\nRanlux64Engine state description incomplete."
722 <<
"\nInput stream is probably mispositioned now." << std::endl;
729 if ((v[0] & 0xffffffffUL) != engineIDulong<Ranlux64Engine>()) {
731 "\nRanlux64Engine get:state vector has wrong ID word - state unchanged\n";
740 "\nRanlux64Engine get:state vector has wrong length - state unchanged\n";
743 std::vector<unsigned long> t(2);
744 for (
int i=0; i<12; ++i) {
745 t[0] = v[2*i+1]; t[1] = v[2*i+2];
748 t[0] = v[25]; t[1] = v[26];
752 pDiscard = (int)v[29];
#define CLHEP_ATOMIC_INT_TYPE
static double longs2double(const std::vector< unsigned long > &v)
static std::vector< unsigned long > dto2longs(double d)
static double twoToMinus_32()
static double twoToMinus_49()
static double twoToMinus_48()
static bool checkFile(std::istream &file, const std::string &filename, const std::string &classname, const std::string &methodname)
static void getTheTableSeeds(long *seeds, int index)
void setSeed(long seed, int lxr=1)
void restoreStatus(const char filename[]="Ranlux64.conf")
std::vector< unsigned long > put() const
virtual std::istream & getState(std::istream &is)
static std::string engineName()
void setSeeds(const long *seeds, int lxr=1)
void saveStatus(const char filename[]="Ranlux64.conf") const
static std::string beginTag()
void flatArray(const int size, double *vect)
static const unsigned int VECTOR_STATE_SIZE
virtual std::istream & get(std::istream &is)
virtual ~Ranlux64Engine()
unsigned long rshift(unsigned long value)
bool possibleKeywordInput(IS &is, const std::string &key, T &t)
unsigned long operator()(unsigned long)