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 #include "CLHEP/Random/DualRand.h"
00055 #include "CLHEP/Random/engineIDulong.h"
00056 #include <string.h>
00057
00058 namespace CLHEP {
00059
00060 static const int MarkerLen = 64;
00061
00062 std::string DualRand::name() const {return "DualRand";}
00063
00064
00065 int DualRand::numEngines = 0;
00066
00067
00068
00069
00070
00071
00072
00073 DualRand::DualRand()
00074 : HepRandomEngine(),
00075 tausworthe (1234567 + numEngines + 175321),
00076 integerCong(69607 * tausworthe + 54329, numEngines)
00077 {
00078 theSeed = 1234567;
00079 ++numEngines;
00080 }
00081
00082 DualRand::DualRand(long seed)
00083 : HepRandomEngine(),
00084 tausworthe ((unsigned int)seed + 175321),
00085 integerCong(69607 * tausworthe + 54329, 8043)
00086 {
00087 theSeed = seed;
00088 }
00089
00090 DualRand::DualRand(std::istream & is)
00091 : HepRandomEngine()
00092 {
00093 is >> *this;
00094 }
00095
00096 DualRand::DualRand(int rowIndex, int colIndex)
00097 : HepRandomEngine(),
00098 tausworthe (rowIndex + 1000 * colIndex + 85329),
00099 integerCong(69607 * tausworthe + 54329, 1123)
00100 {
00101 theSeed = rowIndex;
00102 }
00103
00104 DualRand::~DualRand() { }
00105
00106 double DualRand::flat() {
00107 unsigned int ic ( integerCong );
00108 unsigned int t ( tausworthe );
00109 return ( (t ^ ic) * twoToMinus_32() +
00110 (t >> 11) * twoToMinus_53() +
00111 nearlyTwoToMinus_54()
00112 );
00113 }
00114
00115 void DualRand::flatArray(const int size, double* vect) {
00116 for (int i = 0; i < size; ++i) {
00117 vect[i] = flat();
00118 }
00119 }
00120
00121 void DualRand::setSeed(long seed, int) {
00122 theSeed = seed;
00123 tausworthe = Tausworthe((unsigned int)seed + numEngines + 175321);
00124 integerCong = IntegerCong(69607 * tausworthe + 54329, numEngines);
00125 }
00126
00127 void DualRand::setSeeds(const long * seeds, int) {
00128 setSeed(seeds ? *seeds : 1234567, 0);
00129 theSeeds = seeds;
00130 }
00131
00132 void DualRand::saveStatus(const char filename[]) const {
00133 std::ofstream outFile(filename, std::ios::out);
00134 if (!outFile.bad()) {
00135 outFile << "Uvec\n";
00136 std::vector<unsigned long> v = put();
00137 for (unsigned int i=0; i<v.size(); ++i) {
00138 outFile << v[i] << "\n";
00139 }
00140 }
00141 }
00142
00143 void DualRand::restoreStatus(const char filename[]) {
00144 std::ifstream inFile(filename, std::ios::in);
00145 if (!checkFile ( inFile, filename, engineName(), "restoreStatus" )) {
00146 std::cerr << " -- Engine state remains unchanged\n";
00147 return;
00148 }
00149 if ( possibleKeywordInput ( inFile, "Uvec", theSeed ) ) {
00150 std::vector<unsigned long> v;
00151 unsigned long xin;
00152 for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) {
00153 inFile >> xin;
00154 if (!inFile) {
00155 inFile.clear(std::ios::badbit | inFile.rdstate());
00156 std::cerr << "\nDualRand state (vector) description improper."
00157 << "\nrestoreStatus has failed."
00158 << "\nInput stream is probably mispositioned now." << std::endl;
00159 return;
00160 }
00161 v.push_back(xin);
00162 }
00163 getState(v);
00164 return;
00165 }
00166
00167 if (!inFile.bad()) {
00168
00169 tausworthe.get(inFile);
00170 integerCong.get(inFile);
00171 }
00172 }
00173
00174 void DualRand::showStatus() const {
00175 int pr=std::cout.precision(20);
00176 std::cout << std::endl;
00177 std::cout << "-------- DualRand engine status ---------"
00178 << std::endl;
00179 std::cout << "Initial seed = " << theSeed << std::endl;
00180 std::cout << "Tausworthe generator = " << std::endl;
00181 tausworthe.put(std::cout);
00182 std::cout << "\nIntegerCong generator = " << std::endl;
00183 integerCong.put(std::cout);
00184 std::cout << std::endl << "-----------------------------------------"
00185 << std::endl;
00186 std::cout.precision(pr);
00187 }
00188
00189 DualRand::operator float() {
00190 return (float) ( (integerCong ^ tausworthe) * twoToMinus_32()
00191 + nearlyTwoToMinus_54() );
00192
00193 }
00194
00195 DualRand::operator unsigned int() {
00196 return (integerCong ^ tausworthe) & 0xffffffff;
00197 }
00198
00199 std::ostream & DualRand::put(std::ostream & os) const {
00200 char beginMarker[] = "DualRand-begin";
00201 os << beginMarker << "\nUvec\n";
00202 std::vector<unsigned long> v = put();
00203 for (unsigned int i=0; i<v.size(); ++i) {
00204 os << v[i] << "\n";
00205 }
00206 return os;
00207 }
00208
00209 std::vector<unsigned long> DualRand::put () const {
00210 std::vector<unsigned long> v;
00211 v.push_back (engineIDulong<DualRand>());
00212 tausworthe.put(v);
00213 integerCong.put(v);
00214 return v;
00215 }
00216
00217 std::istream & DualRand::get(std::istream & is) {
00218 char beginMarker [MarkerLen];
00219 is >> std::ws;
00220 is.width(MarkerLen);
00221
00222
00223 is >> beginMarker;
00224 if (strcmp(beginMarker,"DualRand-begin")) {
00225 is.clear(std::ios::badbit | is.rdstate());
00226 std::cerr << "\nInput mispositioned or"
00227 << "\nDualRand state description missing or"
00228 << "\nwrong engine type found." << std::endl;
00229 return is;
00230 }
00231 return getState(is);
00232 }
00233
00234 std::string DualRand::beginTag ( ) {
00235 return "DualRand-begin";
00236 }
00237
00238 std::istream & DualRand::getState ( std::istream & is ) {
00239 if ( possibleKeywordInput ( is, "Uvec", theSeed ) ) {
00240 std::vector<unsigned long> v;
00241 unsigned long uu;
00242 for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) {
00243 is >> uu;
00244 if (!is) {
00245 is.clear(std::ios::badbit | is.rdstate());
00246 std::cerr << "\nDualRand state (vector) description improper."
00247 << "\ngetState() has failed."
00248 << "\nInput stream is probably mispositioned now." << std::endl;
00249 return is;
00250 }
00251 v.push_back(uu);
00252 }
00253 getState(v);
00254 return (is);
00255 }
00256
00257
00258
00259 char endMarker [MarkerLen];
00260 tausworthe.get(is);
00261 integerCong.get(is);
00262 is >> std::ws;
00263 is.width(MarkerLen);
00264 is >> endMarker;
00265 if (strcmp(endMarker,"DualRand-end")) {
00266 is.clear(std::ios::badbit | is.rdstate());
00267 std::cerr << "DualRand state description incomplete."
00268 << "\nInput stream is probably mispositioned now." << std::endl;
00269 return is;
00270 }
00271 return is;
00272 }
00273
00274 bool DualRand::get(const std::vector<unsigned long> & v) {
00275 if ((v[0] & 0xffffffffUL) != engineIDulong<DualRand>()) {
00276 std::cerr <<
00277 "\nDualRand get:state vector has wrong ID word - state unchanged\n";
00278 return false;
00279 }
00280 if (v.size() != VECTOR_STATE_SIZE) {
00281 std::cerr << "\nDualRand get:state vector has wrong size: "
00282 << v.size() << " - state unchanged\n";
00283 return false;
00284 }
00285 return getState(v);
00286 }
00287
00288 bool DualRand::getState (const std::vector<unsigned long> & v) {
00289 std::vector<unsigned long>::const_iterator iv = v.begin()+1;
00290 if (!tausworthe.get(iv)) return false;
00291 if (!integerCong.get(iv)) return false;
00292 if (iv != v.end()) {
00293 std::cerr <<
00294 "\nDualRand get:state vector has wrong size: " << v.size()
00295 << "\n Apparently " << iv-v.begin() << " words were consumed\n";
00296 return false;
00297 }
00298 return true;
00299 }
00300
00301 DualRand::Tausworthe::Tausworthe() {
00302 words[0] = 1234567;
00303 for (wordIndex = 1; wordIndex < 4; ++wordIndex) {
00304 words[wordIndex] = 69607 * words[wordIndex-1] + 54329;
00305 }
00306 }
00307
00308 DualRand::Tausworthe::Tausworthe(unsigned int seed) {
00309 words[0] = seed;
00310 for (wordIndex = 1; wordIndex < 4; ++wordIndex) {
00311 words[wordIndex] = 69607 * words[wordIndex-1] + 54329;
00312 }
00313 }
00314
00315 DualRand::Tausworthe::operator unsigned int() {
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361 if (wordIndex <= 0) {
00362 for (wordIndex = 0; wordIndex < 4; ++wordIndex) {
00363 words[wordIndex] = ( (words[(wordIndex+1) & 3] << 1 ) |
00364 (words[wordIndex] >> 31) )
00365 ^ ( (words[(wordIndex+1) & 3] << 31) |
00366 (words[wordIndex] >> 1) );
00367 }
00368 }
00369 return words[--wordIndex] & 0xffffffff;
00370 }
00371
00372 void DualRand::Tausworthe::put(std::ostream & os) const {
00373 char beginMarker[] = "Tausworthe-begin";
00374 char endMarker[] = "Tausworthe-end";
00375
00376 int pr=os.precision(20);
00377 os << " " << beginMarker << " ";
00378 for (int i = 0; i < 4; ++i) {
00379 os << words[i] << " ";
00380 }
00381 os << wordIndex;
00382 os << " " << endMarker << " ";
00383 os << std::endl;
00384 os.precision(pr);
00385 }
00386
00387 void DualRand::Tausworthe::put(std::vector<unsigned long> & v) const {
00388 for (int i = 0; i < 4; ++i) {
00389 v.push_back(static_cast<unsigned long>(words[i]));
00390 }
00391 v.push_back(static_cast<unsigned long>(wordIndex));
00392 }
00393
00394 void DualRand::Tausworthe::get(std::istream & is) {
00395 char beginMarker [MarkerLen];
00396 char endMarker [MarkerLen];
00397
00398 is >> std::ws;
00399 is.width(MarkerLen);
00400
00401
00402 is >> beginMarker;
00403 if (strcmp(beginMarker,"Tausworthe-begin")) {
00404 is.clear(std::ios::badbit | is.rdstate());
00405 std::cerr << "\nInput mispositioned or"
00406 << "\nTausworthe state description missing or"
00407 << "\nwrong engine type found." << std::endl;
00408 }
00409 for (int i = 0; i < 4; ++i) {
00410 is >> words[i];
00411 }
00412 is >> wordIndex;
00413 is >> std::ws;
00414 is.width(MarkerLen);
00415 is >> endMarker;
00416 if (strcmp(endMarker,"Tausworthe-end")) {
00417 is.clear(std::ios::badbit | is.rdstate());
00418 std::cerr << "\nTausworthe state description incomplete."
00419 << "\nInput stream is probably mispositioned now." << std::endl;
00420 }
00421 }
00422
00423 bool
00424 DualRand::Tausworthe::get(std::vector<unsigned long>::const_iterator & iv){
00425 for (int i = 0; i < 4; ++i) {
00426 words[i] = *iv++;
00427 }
00428 wordIndex = *iv++;
00429 return true;
00430 }
00431
00432 DualRand::IntegerCong::IntegerCong()
00433 : state((unsigned int)3758656018U),
00434 multiplier(66565),
00435 addend(12341)
00436 {
00437 }
00438
00439 DualRand::IntegerCong::IntegerCong(unsigned int seed, int streamNumber)
00440 : state(seed),
00441 multiplier(65536 + 1024 + 5 + (8 * 1017 * streamNumber)),
00442 addend(12341)
00443 {
00444
00445
00446
00447
00448
00449
00450
00451
00452
00453
00454
00455
00456
00457 }
00458
00459 DualRand::IntegerCong::operator unsigned int() {
00460 return state = (state * multiplier + addend) & 0xffffffff;
00461 }
00462
00463 void DualRand::IntegerCong::put(std::ostream & os) const {
00464 char beginMarker[] = "IntegerCong-begin";
00465 char endMarker[] = "IntegerCong-end";
00466
00467 int pr=os.precision(20);
00468 os << " " << beginMarker << " ";
00469 os << state << " " << multiplier << " " << addend;
00470 os << " " << endMarker << " ";
00471 os << std::endl;
00472 os.precision(pr);
00473 }
00474
00475 void DualRand::IntegerCong::put(std::vector<unsigned long> & v) const {
00476 v.push_back(static_cast<unsigned long>(state));
00477 v.push_back(static_cast<unsigned long>(multiplier));
00478 v.push_back(static_cast<unsigned long>(addend));
00479 }
00480
00481 void DualRand::IntegerCong::get(std::istream & is) {
00482 char beginMarker [MarkerLen];
00483 char endMarker [MarkerLen];
00484
00485 is >> std::ws;
00486 is.width(MarkerLen);
00487
00488
00489 is >> beginMarker;
00490 if (strcmp(beginMarker,"IntegerCong-begin")) {
00491 is.clear(std::ios::badbit | is.rdstate());
00492 std::cerr << "\nInput mispositioned or"
00493 << "\nIntegerCong state description missing or"
00494 << "\nwrong engine type found." << std::endl;
00495 }
00496 is >> state >> multiplier >> addend;
00497 is >> std::ws;
00498 is.width(MarkerLen);
00499 is >> endMarker;
00500 if (strcmp(endMarker,"IntegerCong-end")) {
00501 is.clear(std::ios::badbit | is.rdstate());
00502 std::cerr << "\nIntegerCong state description incomplete."
00503 << "\nInput stream is probably mispositioned now." << std::endl;
00504 }
00505 }
00506
00507 bool
00508 DualRand::IntegerCong::get(std::vector<unsigned long>::const_iterator & iv) {
00509 state = *iv++;
00510 multiplier = *iv++;
00511 addend = *iv++;
00512 return true;
00513 }
00514
00515 }