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

RanecuEngine.cc
Go to the documentation of this file.
1// $Id: RanecuEngine.cc,v 1.7 2010/07/20 18:06:02 garren Exp $
2// -*- C++ -*-
3//
4// -----------------------------------------------------------------------
5// HEP Random
6// --- RanecuEngine ---
7// class implementation file
8// -----------------------------------------------------------------------
9// This file is part of Geant4 (simulation toolkit for HEP).
10//
11// RANECU Random Engine - algorithm originally written in FORTRAN77
12// as part of the MATHLIB HEP library.
13
14// =======================================================================
15// Gabriele Cosmo - Created - 2nd February 1996
16// - Minor corrections: 31st October 1996
17// - Added methods for engine status: 19th November 1996
18// - Added abs for setting seed index: 11th July 1997
19// - Modified setSeeds() to handle default index: 16th Oct 1997
20// - setSeed() now resets the engine status to the original
21// values in the static table of HepRandom: 19th Mar 1998
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// M. Fischler - Add endl to the end of saveStatus 10 Apr 2001
28// M. Fischler - In restore, checkFile for file not found 03 Dec 2004
29// M. Fischler - Methods for distrib. instance 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 - put/get for vectors of ulongs 3/14/05
33// M. Fischler - State-saving using only ints, for portability 4/12/05
34// M. Fischler - Modify ctor and setSeed to utilize all info provided
35// and avoid coincidence of same state from different
36// seeds 6/22/10
37//
38// =======================================================================
39
40#include "CLHEP/Random/defs.h"
41#include "CLHEP/Random/Random.h"
42#include "CLHEP/Random/RanecuEngine.h"
43#include "CLHEP/Random/engineIDulong.h"
44#include <string.h> // for strcmp
45#include <cmath>
46#include <cstdlib>
47
48namespace CLHEP {
49
50static const int MarkerLen = 64; // Enough room to hold a begin or end marker.
51
52static const double prec = 4.6566128E-10;
53
54std::string RanecuEngine::name() const {return "RanecuEngine";}
55
56void RanecuEngine::further_randomize (int seq1, int col, int index, int modulus)
57{
58 table[seq1][col] -= (index&0x3FFFFFFF);
59 while (table[seq1][col] <= 0) table[seq1][col] += (modulus-1);
60} // mf 6/22/10
61
62// Number of instances with automatic seed selection
63int RanecuEngine::numEngines = 0;
64
67{
68 int cycle = std::abs(int(numEngines/maxSeq));
69 seq = std::abs(int(numEngines%maxSeq));
70 numEngines += 1;
71 theSeed = seq;
72 long mask = ((cycle & 0x007fffff) << 8);
73 for (int i=0; i<2; ++i) {
74 for (int j=0; j<maxSeq; ++j) {
76 table[j][i] ^= mask;
77 }
78 }
79 theSeeds = &table[seq][0];
80}
81
84{
85 int cycle = std::abs(int(index/maxSeq));
86 seq = std::abs(int(index%maxSeq));
87 theSeed = seq;
88 long mask = ((cycle & 0x000007ff) << 20);
89 for (int j=0; j<maxSeq; ++j) {
91 table[j][0] ^= mask;
92 table[j][1] ^= mask;
93 }
94 theSeeds = &table[seq][0];
95 further_randomize (seq, 0, index, shift1); // mf 6/22/10
96}
97
100{
101 is >> *this;
102}
103
105
106void RanecuEngine::setSeed(long index, int dum)
107{
108 seq = std::abs(int(index%maxSeq));
109 theSeed = seq;
110 HepRandom::getTheTableSeeds(table[seq],seq);
111 theSeeds = &table[seq][0];
112 further_randomize (seq, 0, index, shift1); // mf 6/22/10
113 further_randomize (seq, 1, dum, shift2); // mf 6/22/10
114}
115
116void RanecuEngine::setSeeds(const long* seeds, int pos)
117{
118 if (pos != -1) {
119 seq = std::abs(int(pos%maxSeq));
120 theSeed = seq;
121 }
122 // only positive seeds are allowed
123 table[seq][0] = std::abs(seeds[0])%shift1;
124 table[seq][1] = std::abs(seeds[1])%shift2;
125 theSeeds = &table[seq][0];
126}
127
129{
130 seq = std::abs(int(index%maxSeq));
131 theSeed = seq;
132 theSeeds = &table[seq][0];
133}
134
135void RanecuEngine::saveStatus( const char filename[] ) const
136{
137 std::ofstream outFile( filename, std::ios::out ) ;
138
139 if (!outFile.bad()) {
140 outFile << "Uvec\n";
141 std::vector<unsigned long> v = put();
142 #ifdef TRACE_IO
143 std::cout << "Result of v = put() is:\n";
144 #endif
145 for (unsigned int i=0; i<v.size(); ++i) {
146 outFile << v[i] << "\n";
147 #ifdef TRACE_IO
148 std::cout << v[i] << " ";
149 if (i%6==0) std::cout << "\n";
150 #endif
151 }
152 #ifdef TRACE_IO
153 std::cout << "\n";
154 #endif
155 }
156#ifdef REMOVED
157 if (!outFile.bad()) {
158 outFile << theSeed << std::endl;
159 for (int i=0; i<2; ++i)
160 outFile << table[theSeed][i] << " ";
161 outFile << std::endl;
162 }
163#endif
164}
165
166void RanecuEngine::restoreStatus( const char filename[] )
167{
168 std::ifstream inFile( filename, std::ios::in);
169 if (!checkFile ( inFile, filename, engineName(), "restoreStatus" )) {
170 std::cerr << " -- Engine state remains unchanged\n";
171 return;
172 }
173 if ( possibleKeywordInput ( inFile, "Uvec", theSeed ) ) {
174 std::vector<unsigned long> v;
175 unsigned long xin;
176 for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) {
177 inFile >> xin;
178 #ifdef TRACE_IO
179 std::cout << "ivec = " << ivec << " xin = " << xin << " ";
180 if (ivec%3 == 0) std::cout << "\n";
181 #endif
182 if (!inFile) {
183 inFile.clear(std::ios::badbit | inFile.rdstate());
184 std::cerr << "\nJamesRandom state (vector) description improper."
185 << "\nrestoreStatus has failed."
186 << "\nInput stream is probably mispositioned now." << std::endl;
187 return;
188 }
189 v.push_back(xin);
190 }
191 getState(v);
192 return;
193 }
194
195 if (!inFile.bad() && !inFile.eof()) {
196// inFile >> theSeed; removed -- encompased by possibleKeywordInput
197 for (int i=0; i<2; ++i)
198 inFile >> table[theSeed][i];
199 seq = int(theSeed);
200 }
201}
202
204{
205 std::cout << std::endl;
206 std::cout << "--------- Ranecu engine status ---------" << std::endl;
207 std::cout << " Initial seed (index) = " << theSeed << std::endl;
208 std::cout << " Current couple of seeds = "
209 << table[theSeed][0] << ", "
210 << table[theSeed][1] << std::endl;
211 std::cout << "----------------------------------------" << std::endl;
212}
213
215{
216 const int index = seq;
217 long seed1 = table[index][0];
218 long seed2 = table[index][1];
219
220 int k1 = (int)(seed1/ecuyer_b);
221 int k2 = (int)(seed2/ecuyer_e);
222
223 seed1 = ecuyer_a*(seed1-k1*ecuyer_b)-k1*ecuyer_c;
224 if (seed1 < 0) seed1 += shift1;
225 seed2 = ecuyer_d*(seed2-k2*ecuyer_e)-k2*ecuyer_f;
226 if (seed2 < 0) seed2 += shift2;
227
228 table[index][0] = seed1;
229 table[index][1] = seed2;
230
231 long diff = seed1-seed2;
232
233 if (diff <= 0) diff += (shift1-1);
234 return (double)(diff*prec);
235}
236
237void RanecuEngine::flatArray(const int size, double* vect)
238{
239 const int index = seq;
240 long seed1 = table[index][0];
241 long seed2 = table[index][1];
242 int k1, k2;
243 register int i;
244
245 for (i=0; i<size; ++i)
246 {
247 k1 = (int)(seed1/ecuyer_b);
248 k2 = (int)(seed2/ecuyer_e);
249
250 seed1 = ecuyer_a*(seed1-k1*ecuyer_b)-k1*ecuyer_c;
251 if (seed1 < 0) seed1 += shift1;
252 seed2 = ecuyer_d*(seed2-k2*ecuyer_e)-k2*ecuyer_f;
253 if (seed2 < 0) seed2 += shift2;
254
255 long diff = seed1-seed2;
256 if (diff <= 0) diff += (shift1-1);
257
258 vect[i] = (double)(diff*prec);
259 }
260 table[index][0] = seed1;
261 table[index][1] = seed2;
262}
263
264RanecuEngine::operator unsigned int() {
265 const int index = seq;
266 long seed1 = table[index][0];
267 long seed2 = table[index][1];
268
269 int k1 = (int)(seed1/ecuyer_b);
270 int k2 = (int)(seed2/ecuyer_e);
271
272 seed1 = ecuyer_a*(seed1-k1*ecuyer_b)-k1*ecuyer_c;
273 if (seed1 < 0) seed1 += shift1;
274 seed2 = ecuyer_d*(seed2-k2*ecuyer_e)-k2*ecuyer_f;
275 if (seed2 < 0) seed2 += shift2;
276
277 table[index][0] = seed1;
278 table[index][1] = seed2;
279 long diff = seed1-seed2;
280 if( diff <= 0 ) diff += (shift1-1);
281
282 return ((diff << 1) | (seed1&1))& 0xffffffff;
283}
284
285std::ostream & RanecuEngine::put( std::ostream& os ) const
286{
287 char beginMarker[] = "RanecuEngine-begin";
288 os << beginMarker << "\nUvec\n";
289 std::vector<unsigned long> v = put();
290 for (unsigned int i=0; i<v.size(); ++i) {
291 os << v[i] << "\n";
292 }
293 return os;
294#ifdef REMOVED
295 char endMarker[] = "RanecuEngine-end";
296 os << " " << beginMarker << "\n";
297 os << theSeed << " ";
298 for (int i=0; i<2; ++i) {
299 os << table[theSeed][i] << "\n";
300 }
301 os << endMarker << "\n";
302 return os;
303#endif
304}
305
306std::vector<unsigned long> RanecuEngine::put () const {
307 std::vector<unsigned long> v;
308 v.push_back (engineIDulong<RanecuEngine>());
309 v.push_back(static_cast<unsigned long>(theSeed));
310 v.push_back(static_cast<unsigned long>(table[theSeed][0]));
311 v.push_back(static_cast<unsigned long>(table[theSeed][1]));
312 return v;
313}
314
315std::istream & RanecuEngine::get ( std::istream& is )
316{
317 char beginMarker [MarkerLen];
318
319 is >> std::ws;
320 is.width(MarkerLen); // causes the next read to the char* to be <=
321 // that many bytes, INCLUDING A TERMINATION \0
322 // (Stroustrup, section 21.3.2)
323 is >> beginMarker;
324 if (strcmp(beginMarker,"RanecuEngine-begin")) {
325 is.clear(std::ios::badbit | is.rdstate());
326 std::cerr << "\nInput stream mispositioned or"
327 << "\nRanecuEngine state description missing or"
328 << "\nwrong engine type found." << std::endl;
329 return is;
330 }
331 return getState(is);
332}
333
334std::string RanecuEngine::beginTag ( ) {
335 return "RanecuEngine-begin";
336}
337
338std::istream & RanecuEngine::getState ( std::istream& is )
339{
340 if ( possibleKeywordInput ( is, "Uvec", theSeed ) ) {
341 std::vector<unsigned long> v;
342 unsigned long uu;
343 for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) {
344 is >> uu;
345 if (!is) {
346 is.clear(std::ios::badbit | is.rdstate());
347 std::cerr << "\nRanecuEngine state (vector) description improper."
348 << "\ngetState() has failed."
349 << "\nInput stream is probably mispositioned now." << std::endl;
350 return is;
351 }
352 v.push_back(uu);
353 }
354 getState(v);
355 return (is);
356 }
357
358// is >> theSeed; Removed, encompassed by possibleKeywordInput()
359 char endMarker [MarkerLen];
360 for (int i=0; i<2; ++i) {
361 is >> table[theSeed][i];
362 }
363 is >> std::ws;
364 is.width(MarkerLen);
365 is >> endMarker;
366 if (strcmp(endMarker,"RanecuEngine-end")) {
367 is.clear(std::ios::badbit | is.rdstate());
368 std::cerr << "\nRanecuEngine state description incomplete."
369 << "\nInput stream is probably mispositioned now." << std::endl;
370 return is;
371 }
372
373 seq = int(theSeed);
374 return is;
375}
376
377bool RanecuEngine::get (const std::vector<unsigned long> & v) {
378 if ((v[0] & 0xffffffffUL) != engineIDulong<RanecuEngine>()) {
379 std::cerr <<
380 "\nRanecuEngine get:state vector has wrong ID word - state unchanged\n";
381 return false;
382 }
383 return getState(v);
384}
385
386bool RanecuEngine::getState (const std::vector<unsigned long> & v) {
387 if (v.size() != VECTOR_STATE_SIZE ) {
388 std::cerr <<
389 "\nRanecuEngine get:state vector has wrong length - state unchanged\n";
390 return false;
391 }
392 theSeed = v[1];
393 table[theSeed][0] = v[2];
394 table[theSeed][1] = v[3];
395 seq = int(theSeed);
396 return true;
397}
398
399
400} // namespace CLHEP
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
void restoreStatus(const char filename[]="Ranecu.conf")
void saveStatus(const char filename[]="Ranecu.conf") const
void showStatus() const
std::vector< unsigned long > put() const
void setSeed(long index, int dum=0)
virtual std::istream & get(std::istream &is)
std::string name() const
void flatArray(const int size, double *vect)
virtual std::istream & getState(std::istream &is)
void setIndex(long index)
void setSeeds(const long *seeds, int index=-1)
static const unsigned int VECTOR_STATE_SIZE
static std::string beginTag()
#define double(obj)
bool possibleKeywordInput(IS &is, const std::string &key, T &t)