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/RanecuEngine.h"
00042 #include "CLHEP/Random/engineIDulong.h"
00043 #include <string.h>
00044 #include <cmath>
00045 #include <cstdlib>
00046
00047 namespace CLHEP {
00048
00049 static const int MarkerLen = 64;
00050
00051 static const double prec = 4.6566128E-10;
00052
00053 std::string RanecuEngine::name() const {return "RanecuEngine";}
00054
00055 void RanecuEngine::further_randomize (int seq1, int col, int index, int modulus)
00056 {
00057 table[seq1][col] -= (index&0x3FFFFFFF);
00058 while (table[seq1][col] <= 0) table[seq1][col] += (modulus-1);
00059 }
00060
00061
00062 int RanecuEngine::numEngines = 0;
00063
00064 RanecuEngine::RanecuEngine()
00065 : HepRandomEngine()
00066 {
00067 int cycle = std::abs(int(numEngines/maxSeq));
00068 seq = std::abs(int(numEngines%maxSeq));
00069 numEngines += 1;
00070 theSeed = seq;
00071 long mask = ((cycle & 0x007fffff) << 8);
00072 for (int i=0; i<2; ++i) {
00073 for (int j=0; j<maxSeq; ++j) {
00074 HepRandom::getTheTableSeeds(table[j],j);
00075 table[j][i] ^= mask;
00076 }
00077 }
00078 theSeeds = &table[seq][0];
00079 }
00080
00081 RanecuEngine::RanecuEngine(int index)
00082 : HepRandomEngine()
00083 {
00084 int cycle = std::abs(int(index/maxSeq));
00085 seq = std::abs(int(index%maxSeq));
00086 theSeed = seq;
00087 long mask = ((cycle & 0x000007ff) << 20);
00088 for (int j=0; j<maxSeq; ++j) {
00089 HepRandom::getTheTableSeeds(table[j],j);
00090 table[j][0] ^= mask;
00091 table[j][1] ^= mask;
00092 }
00093 theSeeds = &table[seq][0];
00094 further_randomize (seq, 0, index, shift1);
00095 }
00096
00097 RanecuEngine::RanecuEngine(std::istream& is)
00098 : HepRandomEngine()
00099 {
00100 is >> *this;
00101 }
00102
00103 RanecuEngine::~RanecuEngine() {}
00104
00105 void RanecuEngine::setSeed(long index, int dum)
00106 {
00107 seq = std::abs(int(index%maxSeq));
00108 theSeed = seq;
00109 HepRandom::getTheTableSeeds(table[seq],seq);
00110 theSeeds = &table[seq][0];
00111 further_randomize (seq, 0, index, shift1);
00112 further_randomize (seq, 1, dum, shift2);
00113 }
00114
00115 void RanecuEngine::setSeeds(const long* seeds, int pos)
00116 {
00117 if (pos != -1) {
00118 seq = std::abs(int(pos%maxSeq));
00119 theSeed = seq;
00120 }
00121
00122 table[seq][0] = std::abs(seeds[0])%shift1;
00123 table[seq][1] = std::abs(seeds[1])%shift2;
00124 theSeeds = &table[seq][0];
00125 }
00126
00127 void RanecuEngine::setIndex(long index)
00128 {
00129 seq = std::abs(int(index%maxSeq));
00130 theSeed = seq;
00131 theSeeds = &table[seq][0];
00132 }
00133
00134 void RanecuEngine::saveStatus( const char filename[] ) const
00135 {
00136 std::ofstream outFile( filename, std::ios::out ) ;
00137
00138 if (!outFile.bad()) {
00139 outFile << "Uvec\n";
00140 std::vector<unsigned long> v = put();
00141 for (unsigned int i=0; i<v.size(); ++i) {
00142 outFile << v[i] << "\n";
00143 }
00144 }
00145 }
00146
00147 void RanecuEngine::restoreStatus( const char filename[] )
00148 {
00149 std::ifstream inFile( filename, std::ios::in);
00150 if (!checkFile ( inFile, filename, engineName(), "restoreStatus" )) {
00151 std::cerr << " -- Engine state remains unchanged\n";
00152 return;
00153 }
00154 if ( possibleKeywordInput ( inFile, "Uvec", theSeed ) ) {
00155 std::vector<unsigned long> v;
00156 unsigned long xin;
00157 for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) {
00158 inFile >> xin;
00159 if (!inFile) {
00160 inFile.clear(std::ios::badbit | inFile.rdstate());
00161 std::cerr << "\nJamesRandom state (vector) description improper."
00162 << "\nrestoreStatus has failed."
00163 << "\nInput stream is probably mispositioned now." << std::endl;
00164 return;
00165 }
00166 v.push_back(xin);
00167 }
00168 getState(v);
00169 return;
00170 }
00171
00172 if (!inFile.bad() && !inFile.eof()) {
00173
00174 for (int i=0; i<2; ++i)
00175 inFile >> table[theSeed][i];
00176 seq = int(theSeed);
00177 }
00178 }
00179
00180 void RanecuEngine::showStatus() const
00181 {
00182 std::cout << std::endl;
00183 std::cout << "--------- Ranecu engine status ---------" << std::endl;
00184 std::cout << " Initial seed (index) = " << theSeed << std::endl;
00185 std::cout << " Current couple of seeds = "
00186 << table[theSeed][0] << ", "
00187 << table[theSeed][1] << std::endl;
00188 std::cout << "----------------------------------------" << std::endl;
00189 }
00190
00191 double RanecuEngine::flat()
00192 {
00193 const int index = seq;
00194 long seed1 = table[index][0];
00195 long seed2 = table[index][1];
00196
00197 int k1 = (int)(seed1/ecuyer_b);
00198 int k2 = (int)(seed2/ecuyer_e);
00199
00200 seed1 = ecuyer_a*(seed1-k1*ecuyer_b)-k1*ecuyer_c;
00201 if (seed1 < 0) seed1 += shift1;
00202 seed2 = ecuyer_d*(seed2-k2*ecuyer_e)-k2*ecuyer_f;
00203 if (seed2 < 0) seed2 += shift2;
00204
00205 table[index][0] = seed1;
00206 table[index][1] = seed2;
00207
00208 long diff = seed1-seed2;
00209
00210 if (diff <= 0) diff += (shift1-1);
00211 return (double)(diff*prec);
00212 }
00213
00214 void RanecuEngine::flatArray(const int size, double* vect)
00215 {
00216 const int index = seq;
00217 long seed1 = table[index][0];
00218 long seed2 = table[index][1];
00219 int k1, k2;
00220 register int i;
00221
00222 for (i=0; i<size; ++i)
00223 {
00224 k1 = (int)(seed1/ecuyer_b);
00225 k2 = (int)(seed2/ecuyer_e);
00226
00227 seed1 = ecuyer_a*(seed1-k1*ecuyer_b)-k1*ecuyer_c;
00228 if (seed1 < 0) seed1 += shift1;
00229 seed2 = ecuyer_d*(seed2-k2*ecuyer_e)-k2*ecuyer_f;
00230 if (seed2 < 0) seed2 += shift2;
00231
00232 long diff = seed1-seed2;
00233 if (diff <= 0) diff += (shift1-1);
00234
00235 vect[i] = (double)(diff*prec);
00236 }
00237 table[index][0] = seed1;
00238 table[index][1] = seed2;
00239 }
00240
00241 RanecuEngine::operator unsigned int() {
00242 const int index = seq;
00243 long seed1 = table[index][0];
00244 long seed2 = table[index][1];
00245
00246 int k1 = (int)(seed1/ecuyer_b);
00247 int k2 = (int)(seed2/ecuyer_e);
00248
00249 seed1 = ecuyer_a*(seed1-k1*ecuyer_b)-k1*ecuyer_c;
00250 if (seed1 < 0) seed1 += shift1;
00251 seed2 = ecuyer_d*(seed2-k2*ecuyer_e)-k2*ecuyer_f;
00252 if (seed2 < 0) seed2 += shift2;
00253
00254 table[index][0] = seed1;
00255 table[index][1] = seed2;
00256 long diff = seed1-seed2;
00257 if( diff <= 0 ) diff += (shift1-1);
00258
00259 return ((diff << 1) | (seed1&1))& 0xffffffff;
00260 }
00261
00262 std::ostream & RanecuEngine::put( std::ostream& os ) const
00263 {
00264 char beginMarker[] = "RanecuEngine-begin";
00265 os << beginMarker << "\nUvec\n";
00266 std::vector<unsigned long> v = put();
00267 for (unsigned int i=0; i<v.size(); ++i) {
00268 os << v[i] << "\n";
00269 }
00270 return os;
00271 }
00272
00273 std::vector<unsigned long> RanecuEngine::put () const {
00274 std::vector<unsigned long> v;
00275 v.push_back (engineIDulong<RanecuEngine>());
00276 v.push_back(static_cast<unsigned long>(theSeed));
00277 v.push_back(static_cast<unsigned long>(table[theSeed][0]));
00278 v.push_back(static_cast<unsigned long>(table[theSeed][1]));
00279 return v;
00280 }
00281
00282 std::istream & RanecuEngine::get ( std::istream& is )
00283 {
00284 char beginMarker [MarkerLen];
00285
00286 is >> std::ws;
00287 is.width(MarkerLen);
00288
00289
00290 is >> beginMarker;
00291 if (strcmp(beginMarker,"RanecuEngine-begin")) {
00292 is.clear(std::ios::badbit | is.rdstate());
00293 std::cerr << "\nInput stream mispositioned or"
00294 << "\nRanecuEngine state description missing or"
00295 << "\nwrong engine type found." << std::endl;
00296 return is;
00297 }
00298 return getState(is);
00299 }
00300
00301 std::string RanecuEngine::beginTag ( ) {
00302 return "RanecuEngine-begin";
00303 }
00304
00305 std::istream & RanecuEngine::getState ( std::istream& is )
00306 {
00307 if ( possibleKeywordInput ( is, "Uvec", theSeed ) ) {
00308 std::vector<unsigned long> v;
00309 unsigned long uu;
00310 for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) {
00311 is >> uu;
00312 if (!is) {
00313 is.clear(std::ios::badbit | is.rdstate());
00314 std::cerr << "\nRanecuEngine state (vector) description improper."
00315 << "\ngetState() has failed."
00316 << "\nInput stream is probably mispositioned now." << std::endl;
00317 return is;
00318 }
00319 v.push_back(uu);
00320 }
00321 getState(v);
00322 return (is);
00323 }
00324
00325
00326 char endMarker [MarkerLen];
00327 for (int i=0; i<2; ++i) {
00328 is >> table[theSeed][i];
00329 }
00330 is >> std::ws;
00331 is.width(MarkerLen);
00332 is >> endMarker;
00333 if (strcmp(endMarker,"RanecuEngine-end")) {
00334 is.clear(std::ios::badbit | is.rdstate());
00335 std::cerr << "\nRanecuEngine state description incomplete."
00336 << "\nInput stream is probably mispositioned now." << std::endl;
00337 return is;
00338 }
00339
00340 seq = int(theSeed);
00341 return is;
00342 }
00343
00344 bool RanecuEngine::get (const std::vector<unsigned long> & v) {
00345 if ((v[0] & 0xffffffffUL) != engineIDulong<RanecuEngine>()) {
00346 std::cerr <<
00347 "\nRanecuEngine get:state vector has wrong ID word - state unchanged\n";
00348 return false;
00349 }
00350 return getState(v);
00351 }
00352
00353 bool RanecuEngine::getState (const std::vector<unsigned long> & v) {
00354 if (v.size() != VECTOR_STATE_SIZE ) {
00355 std::cerr <<
00356 "\nRanecuEngine get:state vector has wrong length - state unchanged\n";
00357 return false;
00358 }
00359 theSeed = v[1];
00360 table[theSeed][0] = v[2];
00361 table[theSeed][1] = v[3];
00362 seq = int(theSeed);
00363 return true;
00364 }
00365
00366
00367 }