1 // $Id: JamesRandom.cc,v 1.6 2010/06/16 17:24:53 garren Exp $
2 // -*- C++ -*-
3 //
4 // -----------------------------------------------------------------------
5 // HEP Random
6 // --- HepJamesRandom ---
7 // class implementation file
8 // -----------------------------------------------------------------------
9 // This file is part of Geant4 (simulation toolkit for HEP).
10 //
11 // This algorithm implements the original universal random number generator
12 // as proposed by Marsaglia & Zaman in report FSU-SCRI-87-50 and coded
13 // in FORTRAN77 by Fred James as the RANMAR generator, part of the MATHLIB
14 // HEP library.
15
16 // =======================================================================
17 // Gabriele Cosmo - Created: 5th September 1995
18 // - Fixed a bug in setSeed(): 26th February 1996
19 // - Minor corrections: 31st October 1996
20 // - Added methods for engine status: 19th November 1996
21 // - Fixed bug in setSeeds(): 15th September 1997
22 // J.Marraffino - Added stream operators and related constructor.
23 // Added automatic seed selection from seed table and
24 // engine counter: 16th Feb 1998
25 // Ken Smith - Added conversion operators: 6th Aug 1998
26 // J. Marraffino - Remove dependence on hepString class 13 May 1999
27 // V. Innocente - changed pointers to indices 3 may 2000
28 // M. Fischler - In restore, checkFile for file not found 03 Dec 2004
29 // M. Fischler - Methods for distrib. instacne save/restore 12/8/04
30 // M. Fischler - split get() into tag validation and
31 // getState() for anonymous restores 12/27/04
32 // M. Fischler - Enforcement that seeds be non-negative
33 // (lest the sequence be non-random) 2/14/05
34 // M. Fischler - put/get for vectors of ulongs 3/14/05
35 // M. Fischler - State-saving using only ints, for portability 4/12/05
36 //
37 // =======================================================================
38
39 #include "CLHEP/Random/defs.h"
40 #include "CLHEP/Random/Random.h"
41 #include "CLHEP/Random/JamesRandom.h"
42 #include "CLHEP/Random/engineIDulong.h"
43 #include "CLHEP/Random/DoubConv.hh"
44 #include "CLHEP/Utility/atomic_int.h"
45
46 #include <string.h> // for strcmp
47 #include <cmath>
48 #include <cstdlib>
49
50 //#define TRACE_IO
51
52 namespace CLHEP {
53
54 namespace {
55 // Number of instances with automatic seed selection
56 CLHEP_ATOMIC_INT_TYPE numberOfEngines(0);
57
58 // Maximum index into the seed table
59 const int maxIndex = 215;
60 }
61
62 static const int MarkerLen = 64; // Enough room to hold a begin or end marker.
63
name() const64 std::string HepJamesRandom::name() const {return "HepJamesRandom";}
65
HepJamesRandom(long seed)66 HepJamesRandom::HepJamesRandom(long seed)
67 : HepRandomEngine()
68 {
69 setSeed(seed,0);
70 setSeeds(&theSeed,0);
71 }
72
HepJamesRandom()73 HepJamesRandom::HepJamesRandom() // 15 Feb. 1998 JMM
74 : HepRandomEngine()
75 {
76 long seeds[2];
77 long seed;
78
79 int numEngines = numberOfEngines++;
80 int cycle = std::abs(int(numEngines/maxIndex));
81 int curIndex = std::abs(int(numEngines%maxIndex));
82
83 long mask = ((cycle & 0x007fffff) << 8);
84 HepRandom::getTheTableSeeds( seeds, curIndex );
85 seed = seeds[0]^mask;
86 setSeed(seed,0);
87 setSeeds(&theSeed,0);
88 }
89
HepJamesRandom(int rowIndex,int colIndex)90 HepJamesRandom::HepJamesRandom(int rowIndex, int colIndex) // 15 Feb. 1998 JMM
91 : HepRandomEngine()
92 {
93 long seed;
94 long seeds[2];
95
96 int cycle = std::abs(int(rowIndex/maxIndex));
97 int row = std::abs(int(rowIndex%maxIndex));
98 int col = std::abs(int(colIndex%2));
99 long mask = ((cycle & 0x000007ff) << 20);
100 HepRandom::getTheTableSeeds( seeds, row );
101 seed = (seeds[col])^mask;
102 setSeed(seed,0);
103 setSeeds(&theSeed,0);
104 }
105
HepJamesRandom(std::istream & is)106 HepJamesRandom::HepJamesRandom(std::istream& is)
107 : HepRandomEngine()
108 {
109 is >> *this;
110 }
111
~HepJamesRandom()112 HepJamesRandom::~HepJamesRandom() {}
113
saveStatus(const char filename[]) const114 void HepJamesRandom::saveStatus( const char filename[] ) const
115 {
116 std::ofstream outFile( filename, std::ios::out ) ;
117
118 if (!outFile.bad()) {
119 outFile << "Uvec\n";
120 std::vector<unsigned long> v = put();
121 #ifdef TRACE_IO
122 std::cout << "Result of v = put() is:\n";
123 #endif
124 for (unsigned int i=0; i<v.size(); ++i) {
125 outFile << v[i] << "\n";
126 #ifdef TRACE_IO
127 std::cout << v[i] << " ";
128 if (i%6==0) std::cout << "\n";
129 #endif
130 }
131 #ifdef TRACE_IO
132 std::cout << "\n";
133 #endif
134 }
135 #ifdef REMOVED
136 int pos = j97;
137 outFile << theSeed << std::endl;
138 for (int i=0; i<97; ++i)
139 outFile << std::setprecision(20) << u[i] << " ";
140 outFile << std::endl;
141 outFile << std::setprecision(20) << c << " ";
142 outFile << std::setprecision(20) << cd << " ";
143 outFile << std::setprecision(20) << cm << std::endl;
144 outFile << pos << std::endl;
145 #endif
146 }
147
restoreStatus(const char filename[])148 void HepJamesRandom::restoreStatus( const char filename[] )
149 {
150 int ipos, jpos;
151 std::ifstream inFile( filename, std::ios::in);
152 if (!checkFile ( inFile, filename, engineName(), "restoreStatus" )) {
153 std::cerr << " -- Engine state remains unchanged\n";
154 return;
155 }
156 if ( possibleKeywordInput ( inFile, "Uvec", theSeed ) ) {
157 std::vector<unsigned long> v;
158 unsigned long xin;
159 for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) {
160 inFile >> xin;
161 #ifdef TRACE_IO
162 std::cout << "ivec = " << ivec << " xin = " << xin << " ";
163 if (ivec%3 == 0) std::cout << "\n";
164 #endif
165 if (!inFile) {
166 inFile.clear(std::ios::badbit | inFile.rdstate());
167 std::cerr << "\nJamesRandom state (vector) description improper."
168 << "\nrestoreStatus has failed."
169 << "\nInput stream is probably mispositioned now." << std::endl;
170 return;
171 }
172 v.push_back(xin);
173 }
174 getState(v);
175 return;
176 }
177
178 if (!inFile.bad() && !inFile.eof()) {
179 // inFile >> theSeed; removed -- encompased by possibleKeywordInput
180 for (int i=0; i<97; ++i)
181 inFile >> u[i];
182 inFile >> c; inFile >> cd; inFile >> cm;
183 inFile >> jpos;
184 ipos = (64+jpos)%97;
185 i97 = ipos;
186 j97 = jpos;
187 }
188 }
189
showStatus() const190 void HepJamesRandom::showStatus() const
191 {
192 std::cout << std::endl;
193 std::cout << "----- HepJamesRandom engine status -----" << std::endl;
194 std::cout << " Initial seed = " << theSeed << std::endl;
195 std::cout << " u[] = ";
196 for (int i=0; i<97; ++i)
197 std::cout << u[i] << " ";
198 std::cout << std::endl;
199 std::cout << " c = " << c << ", cd = " << cd << ", cm = " << cm
200 << std::endl;
201 std::cout << " i97 = " << i97 << ", u[i97] = " << u[i97] << std::endl;
202 std::cout << " j97 = " << j97 << ", u[j97] = " << u[j97] << std::endl;
203 std::cout << "----------------------------------------" << std::endl;
204 }
205
setSeed(long seed,int)206 void HepJamesRandom::setSeed(long seed, int)
207 {
208 // The input value for "seed" should be within the range [0,900000000]
209 //
210 // Negative seeds result in serious flaws in the randomness;
211 // seeds above 900000000 are OK because of the %177 in the expression for i,
212 // but may have the same effect as other seeds below 900000000.
213
214 int m, n;
215 float s, t;
216 long mm;
217
218 if (seed < 0) {
219 std::cout << "Seed for HepJamesRandom must be non-negative\n"
220 << "Seed value supplied was " << seed
221 << "\nUsing its absolute value instead\n";
222 seed = -seed;
223 }
224
225 long ij = seed/30082;
226 long kl = seed - 30082*ij;
227 long i = (ij/177) % 177 + 2;
228 long j = ij % 177 + 2;
229 long k = (kl/169) % 178 + 1;
230 long l = kl % 169;
231
232 theSeed = seed;
233
234 for ( n = 1 ; n < 98 ; n++ ) {
235 s = 0.0;
236 t = 0.5;
237 for ( m = 1 ; m < 25 ; m++) {
238 mm = ( ( (i*j) % 179 ) * k ) % 179;
239 i = j;
240 j = k;
241 k = mm;
242 l = ( 53 * l + 1 ) % 169;
243 if ( (l*mm % 64 ) >= 32 )
244 s += t;
245 t *= 0.5;
246 }
247 u[n-1] = s;
248 }
249 c = 362436.0 / 16777216.0;
250 cd = 7654321.0 / 16777216.0;
251 cm = 16777213.0 / 16777216.0;
252
253 i97 = 96;
254 j97 = 32;
255
256 }
257
setSeeds(const long * seeds,int)258 void HepJamesRandom::setSeeds(const long* seeds, int)
259 {
260 setSeed(seeds ? *seeds : 19780503L, 0);
261 theSeeds = seeds;
262 }
263
flat()264 double HepJamesRandom::flat()
265 {
266 double uni;
267
268 do {
269 uni = u[i97] - u[j97];
270 if ( uni < 0.0 ) uni++;
271 u[i97] = uni;
272
273 if (i97 == 0) i97 = 96;
274 else i97--;
275
276 if (j97 == 0) j97 = 96;
277 else j97--;
278
279 c -= cd;
280 if (c < 0.0) c += cm;
281
282 uni -= c;
283 if (uni < 0.0) uni += 1.0;
284 } while ( uni <= 0.0 || uni >= 1.0 );
285
286 return uni;
287 }
288
flatArray(const int size,double * vect)289 void HepJamesRandom::flatArray(const int size, double* vect)
290 {
291 // double uni;
292 int i;
293
294 for (i=0; i<size; ++i) {
295 vect[i] = flat();
296 }
297 }
298
operator double()299 HepJamesRandom::operator double() {
300 return flat();
301 }
302
operator float()303 HepJamesRandom::operator float() {
304 return float( flat() );
305 }
306
operator unsigned int()307 HepJamesRandom::operator unsigned int() {
308 return ((unsigned int)(flat() * exponent_bit_32()) & 0xffffffff ) |
309 (((unsigned int)( u[i97] * exponent_bit_32())>>16) & 0xff);
310 }
311
put(std::ostream & os) const312 std::ostream & HepJamesRandom::put ( std::ostream& os ) const {
313 char beginMarker[] = "JamesRandom-begin";
314 os << beginMarker << "\nUvec\n";
315 std::vector<unsigned long> v = put();
316 for (unsigned int i=0; i<v.size(); ++i) {
317 os << v[i] << "\n";
318 }
319 return os;
320 #ifdef REMOVED
321 char endMarker[] = "JamesRandom-end";
322 int pos = j97;
323 int pr = os.precision(20);
324 os << " " << beginMarker << " ";
325 os << theSeed << " ";
326 for (int i=0; i<97; ++i) {
327 os << std::setprecision(20) << u[i] << "\n";
328 }
329 os << std::setprecision(20) << c << " ";
330 os << std::setprecision(20) << cd << " ";
331 os << std::setprecision(20) << cm << " ";
332 os << pos << "\n";
333 os << endMarker << "\n";
334 os.precision(pr);
335 return os;
336 #endif
337 }
338
put() const339 std::vector<unsigned long> HepJamesRandom::put () const {
340 std::vector<unsigned long> v;
341 v.push_back (engineIDulong<HepJamesRandom>());
342 std::vector<unsigned long> t;
343 for (int i=0; i<97; ++i) {
344 t = DoubConv::dto2longs(u[i]);
345 v.push_back(t[0]); v.push_back(t[1]);
346 }
347 t = DoubConv::dto2longs(c);
348 v.push_back(t[0]); v.push_back(t[1]);
349 t = DoubConv::dto2longs(cd);
350 v.push_back(t[0]); v.push_back(t[1]);
351 t = DoubConv::dto2longs(cm);
352 v.push_back(t[0]); v.push_back(t[1]);
353 v.push_back(static_cast<unsigned long>(j97));
354 return v;
355 }
356
357
get(std::istream & is)358 std::istream & HepJamesRandom::get ( std::istream& is) {
359 char beginMarker [MarkerLen];
360 is >> std::ws;
361 is.width(MarkerLen); // causes the next read to the char* to be <=
362 // that many bytes, INCLUDING A TERMINATION \0
363 // (Stroustrup, section 21.3.2)
364 is >> beginMarker;
365 if (strcmp(beginMarker,"JamesRandom-begin")) {
366 is.clear(std::ios::badbit | is.rdstate());
367 std::cerr << "\nInput stream mispositioned or"
368 << "\nJamesRandom state description missing or"
369 << "\nwrong engine type found." << std::endl;
370 return is;
371 }
372 return getState(is);
373 }
374
beginTag()375 std::string HepJamesRandom::beginTag ( ) {
376 return "JamesRandom-begin";
377 }
378
getState(std::istream & is)379 std::istream & HepJamesRandom::getState ( std::istream& is) {
380 if ( possibleKeywordInput ( is, "Uvec", theSeed ) ) {
381 std::vector<unsigned long> v;
382 unsigned long uu;
383 for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) {
384 is >> uu;
385 if (!is) {
386 is.clear(std::ios::badbit | is.rdstate());
387 std::cerr << "\nJamesRandom state (vector) description improper."
388 << "\ngetState() has failed."
389 << "\nInput stream is probably mispositioned now." << std::endl;
390 return is;
391 }
392 v.push_back(uu);
393 }
394 getState(v);
395 return (is);
396 }
397
398 // is >> theSeed; Removed, encompassed by possibleKeywordInput()
399
400 int ipos, jpos;
401 char endMarker [MarkerLen];
402 for (int i=0; i<97; ++i) {
403 is >> u[i];
404 }
405 is >> c; is >> cd; is >> cm;
406 is >> jpos;
407 is >> std::ws;
408 is.width(MarkerLen);
409 is >> endMarker;
410 if(strcmp(endMarker,"JamesRandom-end")) {
411 is.clear(std::ios::badbit | is.rdstate());
412 std::cerr << "\nJamesRandom state description incomplete."
413 << "\nInput stream is probably mispositioned now." << std::endl;
414 return is;
415 }
416
417 ipos = (64+jpos)%97;
418 i97 = ipos;
419 j97 = jpos;
420 return is;
421 }
422
get(const std::vector<unsigned long> & v)423 bool HepJamesRandom::get (const std::vector<unsigned long> & v) {
424 if ( (v[0] & 0xffffffffUL) != engineIDulong<HepJamesRandom>()) {
425 std::cerr <<
426 "\nHepJamesRandom get:state vector has wrong ID word - state unchanged\n";
427 return false;
428 }
429 return getState(v);
430 }
431
getState(const std::vector<unsigned long> & v)432 bool HepJamesRandom::getState (const std::vector<unsigned long> & v) {
433 if (v.size() != VECTOR_STATE_SIZE ) {
434 std::cerr <<
435 "\nHepJamesRandom get:state vector has wrong length - state unchanged\n";
436 return false;
437 }
438 std::vector<unsigned long> t(2);
439 for (int i=0; i<97; ++i) {
440 t[0] = v[2*i+1]; t[1] = v[2*i+2];
441 u[i] = DoubConv::longs2double(t);
442 }
443 t[0] = v[195]; t[1] = v[196]; c = DoubConv::longs2double(t);
444 t[0] = v[197]; t[1] = v[198]; cd = DoubConv::longs2double(t);
445 t[0] = v[199]; t[1] = v[200]; cm = DoubConv::longs2double(t);
446 j97 = v[201];
447 i97 = (64+j97)%97;
448 return true;
449 }
450
451 } // namespace CLHEP
452