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 #include "CLHEP/Random/Random.h"
00040 #include "CLHEP/Random/JamesRandom.h"
00041 #include "CLHEP/Random/engineIDulong.h"
00042 #include "CLHEP/Random/DoubConv.h"
00043 #include <string.h>
00044 #include <cmath>
00045 #include <cstdlib>
00046
00047 namespace CLHEP {
00048
00049 static const int MarkerLen = 64;
00050
00051 std::string HepJamesRandom::name() const {return "HepJamesRandom";}
00052
00053
00054 int HepJamesRandom::numEngines = 0;
00055
00056
00057 int HepJamesRandom::maxIndex = 215;
00058
00059 HepJamesRandom::HepJamesRandom(long seed)
00060 : HepRandomEngine()
00061 {
00062 setSeed(seed,0);
00063 setSeeds(&theSeed,0);
00064 }
00065
00066 HepJamesRandom::HepJamesRandom()
00067 : HepRandomEngine()
00068 {
00069 long seeds[2];
00070 long seed;
00071
00072 int cycle = std::abs(int(numEngines/maxIndex));
00073 int curIndex = std::abs(int(numEngines%maxIndex));
00074 ++numEngines;
00075 long mask = ((cycle & 0x007fffff) << 8);
00076 HepRandom::getTheTableSeeds( seeds, curIndex );
00077 seed = seeds[0]^mask;
00078 setSeed(seed,0);
00079 setSeeds(&theSeed,0);
00080 }
00081
00082 HepJamesRandom::HepJamesRandom(int rowIndex, int colIndex)
00083 : HepRandomEngine()
00084 {
00085 long seed;
00086 long seeds[2];
00087
00088 int cycle = std::abs(int(rowIndex/maxIndex));
00089 int row = std::abs(int(rowIndex%maxIndex));
00090 int col = std::abs(int(colIndex%2));
00091 long mask = ((cycle & 0x000007ff) << 20);
00092 HepRandom::getTheTableSeeds( seeds, row );
00093 seed = (seeds[col])^mask;
00094 setSeed(seed,0);
00095 setSeeds(&theSeed,0);
00096 }
00097
00098 HepJamesRandom::HepJamesRandom(std::istream& is)
00099 : HepRandomEngine()
00100 {
00101 is >> *this;
00102 }
00103
00104 HepJamesRandom::~HepJamesRandom() {}
00105
00106 void HepJamesRandom::saveStatus( const char filename[] ) const
00107 {
00108 std::ofstream outFile( filename, std::ios::out ) ;
00109
00110 if (!outFile.bad()) {
00111 outFile << "Uvec\n";
00112 std::vector<unsigned long> v = put();
00113 for (unsigned int i=0; i<v.size(); ++i) {
00114 outFile << v[i] << "\n";
00115 }
00116 }
00117 }
00118
00119 void HepJamesRandom::restoreStatus( const char filename[] )
00120 {
00121 int ipos, jpos;
00122 std::ifstream inFile( filename, std::ios::in);
00123 if (!checkFile ( inFile, filename, engineName(), "restoreStatus" )) {
00124 std::cerr << " -- Engine state remains unchanged\n";
00125 return;
00126 }
00127 if ( possibleKeywordInput ( inFile, "Uvec", theSeed ) ) {
00128 std::vector<unsigned long> v;
00129 unsigned long xin;
00130 for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) {
00131 inFile >> xin;
00132 if (!inFile) {
00133 inFile.clear(std::ios::badbit | inFile.rdstate());
00134 std::cerr << "\nJamesRandom state (vector) description improper."
00135 << "\nrestoreStatus has failed."
00136 << "\nInput stream is probably mispositioned now." << std::endl;
00137 return;
00138 }
00139 v.push_back(xin);
00140 }
00141 getState(v);
00142 return;
00143 }
00144
00145 if (!inFile.bad() && !inFile.eof()) {
00146
00147 for (int i=0; i<97; ++i)
00148 inFile >> u[i];
00149 inFile >> c; inFile >> cd; inFile >> cm;
00150 inFile >> jpos;
00151 ipos = (64+jpos)%97;
00152 i97 = ipos;
00153 j97 = jpos;
00154 }
00155 }
00156
00157 void HepJamesRandom::showStatus() const
00158 {
00159 std::cout << std::endl;
00160 std::cout << "----- HepJamesRandom engine status -----" << std::endl;
00161 std::cout << " Initial seed = " << theSeed << std::endl;
00162 std::cout << " u[] = ";
00163 for (int i=0; i<97; ++i)
00164 std::cout << u[i] << " ";
00165 std::cout << std::endl;
00166 std::cout << " c = " << c << ", cd = " << cd << ", cm = " << cm
00167 << std::endl;
00168 std::cout << " i97 = " << i97 << ", u[i97] = " << u[i97] << std::endl;
00169 std::cout << " j97 = " << j97 << ", u[j97] = " << u[j97] << std::endl;
00170 std::cout << "----------------------------------------" << std::endl;
00171 }
00172
00173 void HepJamesRandom::setSeed(long seed, int)
00174 {
00175
00176
00177
00178
00179
00180
00181 int m, n;
00182 float s, t;
00183 long mm;
00184
00185 if (seed < 0) {
00186 std::cout << "Seed for HepJamesRandom must be non-negative\n"
00187 << "Seed value supplied was " << seed
00188 << "\nUsing its absolute value instead\n";
00189 seed = -seed;
00190 }
00191
00192 long ij = seed/30082;
00193 long kl = seed - 30082*ij;
00194 long i = (ij/177) % 177 + 2;
00195 long j = ij % 177 + 2;
00196 long k = (kl/169) % 178 + 1;
00197 long l = kl % 169;
00198
00199 theSeed = seed;
00200
00201 for ( n = 1 ; n < 98 ; n++ ) {
00202 s = 0.0;
00203 t = 0.5;
00204 for ( m = 1 ; m < 25 ; m++) {
00205 mm = ( ( (i*j) % 179 ) * k ) % 179;
00206 i = j;
00207 j = k;
00208 k = mm;
00209 l = ( 53 * l + 1 ) % 169;
00210 if ( (l*mm % 64 ) >= 32 )
00211 s += t;
00212 t *= 0.5;
00213 }
00214 u[n-1] = s;
00215 }
00216 c = 362436.0 / 16777216.0;
00217 cd = 7654321.0 / 16777216.0;
00218 cm = 16777213.0 / 16777216.0;
00219
00220 i97 = 96;
00221 j97 = 32;
00222
00223 }
00224
00225 void HepJamesRandom::setSeeds(const long* seeds, int)
00226 {
00227 setSeed(seeds ? *seeds : 19780503L, 0);
00228 theSeeds = seeds;
00229 }
00230
00231 double HepJamesRandom::flat()
00232 {
00233 double uni;
00234
00235 do {
00236 uni = u[i97] - u[j97];
00237 if ( uni < 0.0 ) uni++;
00238 u[i97] = uni;
00239
00240 if (i97 == 0) i97 = 96;
00241 else i97--;
00242
00243 if (j97 == 0) j97 = 96;
00244 else j97--;
00245
00246 c -= cd;
00247 if (c < 0.0) c += cm;
00248
00249 uni -= c;
00250 if (uni < 0.0) uni += 1.0;
00251 } while ( uni <= 0.0 || uni >= 1.0 );
00252
00253 return uni;
00254 }
00255
00256 void HepJamesRandom::flatArray(const int size, double* vect)
00257 {
00258
00259 int i;
00260
00261 for (i=0; i<size; ++i) {
00262 vect[i] = flat();
00263 }
00264 }
00265
00266 HepJamesRandom::operator unsigned int() {
00267 return ((unsigned int)(flat() * exponent_bit_32()) & 0xffffffff ) |
00268 (((unsigned int)( u[i97] * exponent_bit_32())>>16) & 0xff);
00269 }
00270
00271 std::ostream & HepJamesRandom::put ( std::ostream& os ) const {
00272 char beginMarker[] = "JamesRandom-begin";
00273 os << beginMarker << "\nUvec\n";
00274 std::vector<unsigned long> v = put();
00275 for (unsigned int i=0; i<v.size(); ++i) {
00276 os << v[i] << "\n";
00277 }
00278 return os;
00279 }
00280
00281 std::vector<unsigned long> HepJamesRandom::put () const {
00282 std::vector<unsigned long> v;
00283 v.push_back (engineIDulong<HepJamesRandom>());
00284 std::vector<unsigned long> t;
00285 for (int i=0; i<97; ++i) {
00286 t = DoubConv::dto2longs(u[i]);
00287 v.push_back(t[0]); v.push_back(t[1]);
00288 }
00289 t = DoubConv::dto2longs(c);
00290 v.push_back(t[0]); v.push_back(t[1]);
00291 t = DoubConv::dto2longs(cd);
00292 v.push_back(t[0]); v.push_back(t[1]);
00293 t = DoubConv::dto2longs(cm);
00294 v.push_back(t[0]); v.push_back(t[1]);
00295 v.push_back(static_cast<unsigned long>(j97));
00296 return v;
00297 }
00298
00299
00300 std::istream & HepJamesRandom::get ( std::istream& is) {
00301 char beginMarker [MarkerLen];
00302 is >> std::ws;
00303 is.width(MarkerLen);
00304
00305
00306 is >> beginMarker;
00307 if (strcmp(beginMarker,"JamesRandom-begin")) {
00308 is.clear(std::ios::badbit | is.rdstate());
00309 std::cerr << "\nInput stream mispositioned or"
00310 << "\nJamesRandom state description missing or"
00311 << "\nwrong engine type found." << std::endl;
00312 return is;
00313 }
00314 return getState(is);
00315 }
00316
00317 std::string HepJamesRandom::beginTag ( ) {
00318 return "JamesRandom-begin";
00319 }
00320
00321 std::istream & HepJamesRandom::getState ( std::istream& is) {
00322 if ( possibleKeywordInput ( is, "Uvec", theSeed ) ) {
00323 std::vector<unsigned long> v;
00324 unsigned long uu;
00325 for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) {
00326 is >> uu;
00327 if (!is) {
00328 is.clear(std::ios::badbit | is.rdstate());
00329 std::cerr << "\nJamesRandom state (vector) description improper."
00330 << "\ngetState() has failed."
00331 << "\nInput stream is probably mispositioned now." << std::endl;
00332 return is;
00333 }
00334 v.push_back(uu);
00335 }
00336 getState(v);
00337 return (is);
00338 }
00339
00340
00341
00342 int ipos, jpos;
00343 char endMarker [MarkerLen];
00344 for (int i=0; i<97; ++i) {
00345 is >> u[i];
00346 }
00347 is >> c; is >> cd; is >> cm;
00348 is >> jpos;
00349 is >> std::ws;
00350 is.width(MarkerLen);
00351 is >> endMarker;
00352 if(strcmp(endMarker,"JamesRandom-end")) {
00353 is.clear(std::ios::badbit | is.rdstate());
00354 std::cerr << "\nJamesRandom state description incomplete."
00355 << "\nInput stream is probably mispositioned now." << std::endl;
00356 return is;
00357 }
00358
00359 ipos = (64+jpos)%97;
00360 i97 = ipos;
00361 j97 = jpos;
00362 return is;
00363 }
00364
00365 bool HepJamesRandom::get (const std::vector<unsigned long> & v) {
00366 if ( (v[0] & 0xffffffffUL) != engineIDulong<HepJamesRandom>()) {
00367 std::cerr <<
00368 "\nHepJamesRandom get:state vector has wrong ID word - state unchanged\n";
00369 return false;
00370 }
00371 return getState(v);
00372 }
00373
00374 bool HepJamesRandom::getState (const std::vector<unsigned long> & v) {
00375 if (v.size() != VECTOR_STATE_SIZE ) {
00376 std::cerr <<
00377 "\nHepJamesRandom get:state vector has wrong length - state unchanged\n";
00378 return false;
00379 }
00380 std::vector<unsigned long> t(2);
00381 for (int i=0; i<97; ++i) {
00382 t[0] = v[2*i+1]; t[1] = v[2*i+2];
00383 u[i] = DoubConv::longs2double(t);
00384 }
00385 t[0] = v[195]; t[1] = v[196]; c = DoubConv::longs2double(t);
00386 t[0] = v[197]; t[1] = v[198]; cd = DoubConv::longs2double(t);
00387 t[0] = v[199]; t[1] = v[200]; cm = DoubConv::longs2double(t);
00388 j97 = v[201];
00389 i97 = (64+j97)%97;
00390 return true;
00391 }
00392
00393 }