CLHEP VERSION Reference Documentation
   
CLHEP Home Page     CLHEP Documentation     CLHEP Bug Reports

JamesRandom.cc
Go to the documentation of this file.
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 <string.h> // for strcmp
45#include <cmath>
46#include <cstdlib>
47
48//#define TRACE_IO
49
50namespace CLHEP {
51
52static const int MarkerLen = 64; // Enough room to hold a begin or end marker.
53
54std::string HepJamesRandom::name() const {return "HepJamesRandom";}
55
56// Number of instances with automatic seed selection
57int HepJamesRandom::numEngines = 0;
58
59// Maximum index into the seed table
60int HepJamesRandom::maxIndex = 215;
61
64{
65 setSeed(seed,0);
66 setSeeds(&theSeed,0);
67}
68
69HepJamesRandom::HepJamesRandom() // 15 Feb. 1998 JMM
71{
72 long seeds[2];
73 long seed;
74
75 int cycle = std::abs(int(numEngines/maxIndex));
76 int curIndex = std::abs(int(numEngines%maxIndex));
77 ++numEngines;
78 long mask = ((cycle & 0x007fffff) << 8);
79 HepRandom::getTheTableSeeds( seeds, curIndex );
80 seed = seeds[0]^mask;
81 setSeed(seed,0);
82 setSeeds(&theSeed,0);
83}
84
85HepJamesRandom::HepJamesRandom(int rowIndex, int colIndex) // 15 Feb. 1998 JMM
87{
88 long seed;
89 long seeds[2];
90
91 int cycle = std::abs(int(rowIndex/maxIndex));
92 int row = std::abs(int(rowIndex%maxIndex));
93 int col = std::abs(int(colIndex%2));
94 long mask = ((cycle & 0x000007ff) << 20);
95 HepRandom::getTheTableSeeds( seeds, row );
96 seed = (seeds[col])^mask;
97 setSeed(seed,0);
98 setSeeds(&theSeed,0);
99}
100
103{
104 is >> *this;
105}
106
108
109void HepJamesRandom::saveStatus( const char filename[] ) const
110{
111 std::ofstream outFile( filename, std::ios::out ) ;
112
113 if (!outFile.bad()) {
114 outFile << "Uvec\n";
115 std::vector<unsigned long> v = put();
116 #ifdef TRACE_IO
117 std::cout << "Result of v = put() is:\n";
118 #endif
119 for (unsigned int i=0; i<v.size(); ++i) {
120 outFile << v[i] << "\n";
121 #ifdef TRACE_IO
122 std::cout << v[i] << " ";
123 if (i%6==0) std::cout << "\n";
124 #endif
125 }
126 #ifdef TRACE_IO
127 std::cout << "\n";
128 #endif
129 }
130#ifdef REMOVED
131 int pos = j97;
132 outFile << theSeed << std::endl;
133 for (int i=0; i<97; ++i)
134 outFile << std::setprecision(20) << u[i] << " ";
135 outFile << std::endl;
136 outFile << std::setprecision(20) << c << " ";
137 outFile << std::setprecision(20) << cd << " ";
138 outFile << std::setprecision(20) << cm << std::endl;
139 outFile << pos << std::endl;
140#endif
141}
142
143void HepJamesRandom::restoreStatus( const char filename[] )
144{
145 int ipos, jpos;
146 std::ifstream inFile( filename, std::ios::in);
147 if (!checkFile ( inFile, filename, engineName(), "restoreStatus" )) {
148 std::cerr << " -- Engine state remains unchanged\n";
149 return;
150 }
151 if ( possibleKeywordInput ( inFile, "Uvec", theSeed ) ) {
152 std::vector<unsigned long> v;
153 unsigned long xin;
154 for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) {
155 inFile >> xin;
156 #ifdef TRACE_IO
157 std::cout << "ivec = " << ivec << " xin = " << xin << " ";
158 if (ivec%3 == 0) std::cout << "\n";
159 #endif
160 if (!inFile) {
161 inFile.clear(std::ios::badbit | inFile.rdstate());
162 std::cerr << "\nJamesRandom state (vector) description improper."
163 << "\nrestoreStatus has failed."
164 << "\nInput stream is probably mispositioned now." << std::endl;
165 return;
166 }
167 v.push_back(xin);
168 }
169 getState(v);
170 return;
171 }
172
173 if (!inFile.bad() && !inFile.eof()) {
174// inFile >> theSeed; removed -- encompased by possibleKeywordInput
175 for (int i=0; i<97; ++i)
176 inFile >> u[i];
177 inFile >> c; inFile >> cd; inFile >> cm;
178 inFile >> jpos;
179 ipos = (64+jpos)%97;
180 i97 = ipos;
181 j97 = jpos;
182 }
183}
184
186{
187 std::cout << std::endl;
188 std::cout << "----- HepJamesRandom engine status -----" << std::endl;
189 std::cout << " Initial seed = " << theSeed << std::endl;
190 std::cout << " u[] = ";
191 for (int i=0; i<97; ++i)
192 std::cout << u[i] << " ";
193 std::cout << std::endl;
194 std::cout << " c = " << c << ", cd = " << cd << ", cm = " << cm
195 << std::endl;
196 std::cout << " i97 = " << i97 << ", u[i97] = " << u[i97] << std::endl;
197 std::cout << " j97 = " << j97 << ", u[j97] = " << u[j97] << std::endl;
198 std::cout << "----------------------------------------" << std::endl;
199}
200
201void HepJamesRandom::setSeed(long seed, int)
202{
203 // The input value for "seed" should be within the range [0,900000000]
204 //
205 // Negative seeds result in serious flaws in the randomness;
206 // seeds above 900000000 are OK because of the %177 in the expression for i,
207 // but may have the same effect as other seeds below 900000000.
208
209 int m, n;
210 float s, t;
211 long mm;
212
213 if (seed < 0) {
214 std::cout << "Seed for HepJamesRandom must be non-negative\n"
215 << "Seed value supplied was " << seed
216 << "\nUsing its absolute value instead\n";
217 seed = -seed;
218 }
219
220 long ij = seed/30082;
221 long kl = seed - 30082*ij;
222 long i = (ij/177) % 177 + 2;
223 long j = ij % 177 + 2;
224 long k = (kl/169) % 178 + 1;
225 long l = kl % 169;
226
227 theSeed = seed;
228
229 for ( n = 1 ; n < 98 ; n++ ) {
230 s = 0.0;
231 t = 0.5;
232 for ( m = 1 ; m < 25 ; m++) {
233 mm = ( ( (i*j) % 179 ) * k ) % 179;
234 i = j;
235 j = k;
236 k = mm;
237 l = ( 53 * l + 1 ) % 169;
238 if ( (l*mm % 64 ) >= 32 )
239 s += t;
240 t *= 0.5;
241 }
242 u[n-1] = s;
243 }
244 c = 362436.0 / 16777216.0;
245 cd = 7654321.0 / 16777216.0;
246 cm = 16777213.0 / 16777216.0;
247
248 i97 = 96;
249 j97 = 32;
250
251}
252
253void HepJamesRandom::setSeeds(const long* seeds, int)
254{
255 setSeed(seeds ? *seeds : 19780503L, 0);
256 theSeeds = seeds;
257}
258
260{
261 double uni;
262
263 do {
264 uni = u[i97] - u[j97];
265 if ( uni < 0.0 ) uni++;
266 u[i97] = uni;
267
268 if (i97 == 0) i97 = 96;
269 else i97--;
270
271 if (j97 == 0) j97 = 96;
272 else j97--;
273
274 c -= cd;
275 if (c < 0.0) c += cm;
276
277 uni -= c;
278 if (uni < 0.0) uni += 1.0;
279 } while ( uni <= 0.0 || uni >= 1.0 );
280
281 return uni;
282}
283
284void HepJamesRandom::flatArray(const int size, double* vect)
285{
286// double uni;
287 int i;
288
289 for (i=0; i<size; ++i) {
290 vect[i] = flat();
291 }
292}
293
294HepJamesRandom::operator unsigned int() {
295 return ((unsigned int)(flat() * exponent_bit_32()) & 0xffffffff ) |
296 (((unsigned int)( u[i97] * exponent_bit_32())>>16) & 0xff);
297}
298
299std::ostream & HepJamesRandom::put ( std::ostream& os ) const {
300 char beginMarker[] = "JamesRandom-begin";
301 os << beginMarker << "\nUvec\n";
302 std::vector<unsigned long> v = put();
303 for (unsigned int i=0; i<v.size(); ++i) {
304 os << v[i] << "\n";
305 }
306 return os;
307#ifdef REMOVED
308 char endMarker[] = "JamesRandom-end";
309 int pos = j97;
310 int pr = os.precision(20);
311 os << " " << beginMarker << " ";
312 os << theSeed << " ";
313 for (int i=0; i<97; ++i) {
314 os << std::setprecision(20) << u[i] << "\n";
315 }
316 os << std::setprecision(20) << c << " ";
317 os << std::setprecision(20) << cd << " ";
318 os << std::setprecision(20) << cm << " ";
319 os << pos << "\n";
320 os << endMarker << "\n";
321 os.precision(pr);
322 return os;
323#endif
324}
325
326std::vector<unsigned long> HepJamesRandom::put () const {
327 std::vector<unsigned long> v;
328 v.push_back (engineIDulong<HepJamesRandom>());
329 std::vector<unsigned long> t;
330 for (int i=0; i<97; ++i) {
331 t = DoubConv::dto2longs(u[i]);
332 v.push_back(t[0]); v.push_back(t[1]);
333 }
334 t = DoubConv::dto2longs(c);
335 v.push_back(t[0]); v.push_back(t[1]);
336 t = DoubConv::dto2longs(cd);
337 v.push_back(t[0]); v.push_back(t[1]);
338 t = DoubConv::dto2longs(cm);
339 v.push_back(t[0]); v.push_back(t[1]);
340 v.push_back(static_cast<unsigned long>(j97));
341 return v;
342}
343
344
345std::istream & HepJamesRandom::get ( std::istream& is) {
346 char beginMarker [MarkerLen];
347 is >> std::ws;
348 is.width(MarkerLen); // causes the next read to the char* to be <=
349 // that many bytes, INCLUDING A TERMINATION \0
350 // (Stroustrup, section 21.3.2)
351 is >> beginMarker;
352 if (strcmp(beginMarker,"JamesRandom-begin")) {
353 is.clear(std::ios::badbit | is.rdstate());
354 std::cerr << "\nInput stream mispositioned or"
355 << "\nJamesRandom state description missing or"
356 << "\nwrong engine type found." << std::endl;
357 return is;
358 }
359 return getState(is);
360}
361
362std::string HepJamesRandom::beginTag ( ) {
363 return "JamesRandom-begin";
364}
365
366std::istream & HepJamesRandom::getState ( std::istream& is) {
367 if ( possibleKeywordInput ( is, "Uvec", theSeed ) ) {
368 std::vector<unsigned long> v;
369 unsigned long uu;
370 for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) {
371 is >> uu;
372 if (!is) {
373 is.clear(std::ios::badbit | is.rdstate());
374 std::cerr << "\nJamesRandom state (vector) description improper."
375 << "\ngetState() has failed."
376 << "\nInput stream is probably mispositioned now." << std::endl;
377 return is;
378 }
379 v.push_back(uu);
380 }
381 getState(v);
382 return (is);
383 }
384
385// is >> theSeed; Removed, encompassed by possibleKeywordInput()
386
387 int ipos, jpos;
388 char endMarker [MarkerLen];
389 for (int i=0; i<97; ++i) {
390 is >> u[i];
391 }
392 is >> c; is >> cd; is >> cm;
393 is >> jpos;
394 is >> std::ws;
395 is.width(MarkerLen);
396 is >> endMarker;
397 if(strcmp(endMarker,"JamesRandom-end")) {
398 is.clear(std::ios::badbit | is.rdstate());
399 std::cerr << "\nJamesRandom state description incomplete."
400 << "\nInput stream is probably mispositioned now." << std::endl;
401 return is;
402 }
403
404 ipos = (64+jpos)%97;
405 i97 = ipos;
406 j97 = jpos;
407 return is;
408}
409
410bool HepJamesRandom::get (const std::vector<unsigned long> & v) {
411 if ( (v[0] & 0xffffffffUL) != engineIDulong<HepJamesRandom>()) {
412 std::cerr <<
413 "\nHepJamesRandom get:state vector has wrong ID word - state unchanged\n";
414 return false;
415 }
416 return getState(v);
417}
418
419bool HepJamesRandom::getState (const std::vector<unsigned long> & v) {
420 if (v.size() != VECTOR_STATE_SIZE ) {
421 std::cerr <<
422 "\nHepJamesRandom get:state vector has wrong length - state unchanged\n";
423 return false;
424 }
425 std::vector<unsigned long> t(2);
426 for (int i=0; i<97; ++i) {
427 t[0] = v[2*i+1]; t[1] = v[2*i+2];
428 u[i] = DoubConv::longs2double(t);
429 }
430 t[0] = v[195]; t[1] = v[196]; c = DoubConv::longs2double(t);
431 t[0] = v[197]; t[1] = v[198]; cd = DoubConv::longs2double(t);
432 t[0] = v[199]; t[1] = v[200]; cm = DoubConv::longs2double(t);
433 j97 = v[201];
434 i97 = (64+j97)%97;
435 return true;
436}
437
438} // namespace CLHEP
static double longs2double(const std::vector< unsigned long > &v)
Definition DoubConv.cc:106
static std::vector< unsigned long > dto2longs(double d)
Definition DoubConv.cc:90
void flatArray(const int size, double *vect)
static std::string beginTag()
void setSeeds(const long *seeds, int dum=0)
static const unsigned int VECTOR_STATE_SIZE
virtual std::istream & getState(std::istream &is)
std::string name() const
std::vector< unsigned long > put() const
void setSeed(long seed, int dum=0)
virtual std::istream & get(std::istream &is)
void saveStatus(const char filename[]="JamesRand.conf") const
void restoreStatus(const char filename[]="JamesRand.conf")
void showStatus() const
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)
Definition Random.cc:152
bool possibleKeywordInput(IS &is, const std::string &key, T &t)