| Geant4 Cross Reference |
1 // -*- C++ -*- 1
2 //
3 // -------------------------------------------
4 // HEP Random
5 // --- RanecuEngine ---
6 // class implementation f
7 // -------------------------------------------
8 // This file is part of Geant4 (simulation too
9 //
10 // RANECU Random Engine - algorithm originally
11 // as part of the MATHL
12
13 // ===========================================
14 // Gabriele Cosmo - Created - 2nd February 199
15 // - Minor corrections: 31st Oc
16 // - Added methods for engine s
17 // - Added abs for setting seed
18 // - Modified setSeeds() to han
19 // - setSeed() now resets the e
20 // values in the static table
21 // J.Marraffino - Added stream operators and
22 // Added automatic seed selec
23 // engine counter: 16th Feb 1
24 // Ken Smith - Added conversion operators
25 // J. Marraffino - Remove dependence on hepSt
26 // M. Fischler - Add endl to the end of sav
27 // M. Fischler - In restore, checkFile for
28 // M. Fischler - Methods for distrib. insta
29 // M. Fischler - split get() into tag valid
30 // getState() for anonymous r
31 // M. Fischler - put/get for vectors of ulo
32 // M. Fischler - State-saving using only in
33 // M. Fischler - Modify ctor and setSeed to
34 // and avoid coincidence of s
35 // seeds
36 //
37 // ===========================================
38
39 #include "CLHEP/Random/Random.h"
40 #include "CLHEP/Random/RanecuEngine.h"
41 #include "CLHEP/Random/engineIDulong.h"
42 #include "CLHEP/Utility/atomic_int.h"
43
44 #include <atomic>
45 #include <cstdlib>
46 #include <cmath>
47 #include <iostream>
48 #include <string>
49 #include <string.h> // for strcmp
50 #include <vector>
51
52 namespace CLHEP {
53
54 namespace {
55 // Number of instances with automatic seed s
56 CLHEP_ATOMIC_INT_TYPE numberOfEngines(0);
57 }
58
59 static const int MarkerLen = 64; // Enough roo
60
61 static const double prec = 4.6566128E-10;
62
63 std::string RanecuEngine::name() const {return
64
65 void RanecuEngine::further_randomize (int seq1
66 {
67 table[seq1][col] -= (index&0x3FFFFFFF);
68 while (table[seq1][col] <= 0) table[seq1][co
69 } // mf 6/22/10
70
71 RanecuEngine::RanecuEngine()
72 : HepRandomEngine()
73 {
74 int numEngines = numberOfEngines++;
75 int cycle = std::abs(int(numEngines/maxSeq))
76 seq = std::abs(int(numEngines%maxSeq));
77
78 theSeed = seq;
79 long mask = ((cycle & 0x007fffff) << 8);
80 for (int i=0; i<2; ++i) {
81 for (int j=0; j<maxSeq; ++j) {
82 HepRandom::getTheTableSeeds(table[j],j);
83 table[j][i] ^= mask;
84 }
85 }
86 theSeeds = &table[seq][0];
87 }
88
89 RanecuEngine::RanecuEngine(int index)
90 : HepRandomEngine()
91 {
92 int cycle = std::abs(int(index/maxSeq));
93 seq = std::abs(int(index%maxSeq));
94 theSeed = seq;
95 long mask = ((cycle & 0x000007ff) << 20);
96 for (int j=0; j<maxSeq; ++j) {
97 HepRandom::getTheTableSeeds(table[j],j);
98 table[j][0] ^= mask;
99 table[j][1] ^= mask;
100 }
101 theSeeds = &table[seq][0];
102 further_randomize (seq, 0, index, shift1);
103 }
104
105 RanecuEngine::RanecuEngine(std::istream& is)
106 : HepRandomEngine()
107 {
108 is >> *this;
109 }
110
111 RanecuEngine::~RanecuEngine() {}
112
113 void RanecuEngine::setSeed(long index, int dum
114 {
115 seq = std::abs(int(index%maxSeq));
116 theSeed = seq;
117 HepRandom::getTheTableSeeds(table[seq],seq);
118 theSeeds = &table[seq][0];
119 further_randomize (seq, 0, (int)index, shift
120 further_randomize (seq, 1, dum, shift2);
121 }
122
123 void RanecuEngine::setSeeds(const long* seeds,
124 {
125 if (pos != -1) {
126 seq = std::abs(int(pos%maxSeq));
127 theSeed = seq;
128 }
129 // only positive seeds are allowed
130 table[seq][0] = std::abs(seeds[0])%shift1;
131 table[seq][1] = std::abs(seeds[1])%shift2;
132 theSeeds = &table[seq][0];
133 }
134
135 void RanecuEngine::setIndex(long index)
136 {
137 seq = std::abs(int(index%maxSeq));
138 theSeed = seq;
139 theSeeds = &table[seq][0];
140 }
141
142 void RanecuEngine::saveStatus( const char file
143 {
144 std::ofstream outFile( filename, std::ios::
145
146 if (!outFile.bad()) {
147 outFile << "Uvec\n";
148 std::vector<unsigned long> v = put();
149 for (unsigned int i=0; i<v.size(); ++i) {
150 outFile << v[i] << "\n";
151 }
152 }
153 }
154
155 void RanecuEngine::restoreStatus( const char f
156 {
157 std::ifstream inFile( filename, std::ios::in
158 if (!checkFile ( inFile, filename, engineNam
159 std::cerr << " -- Engine state remains un
160 return;
161 }
162 if ( possibleKeywordInput ( inFile, "Uvec",
163 std::vector<unsigned long> v;
164 unsigned long xin;
165 for (unsigned int ivec=0; ivec < VECTOR_ST
166 inFile >> xin;
167 if (!inFile) {
168 inFile.clear(std::ios::badbit | inFile
169 std::cerr << "\nJamesRandom state (vec
170 << "\nrestoreStatus has failed."
171 << "\nInput stream is probably mispos
172 return;
173 }
174 v.push_back(xin);
175 }
176 getState(v);
177 return;
178 }
179
180 if (!inFile.bad() && !inFile.eof()) {
181 // inFile >> theSeed; removed -- encompas
182 for (int i=0; i<2; ++i)
183 inFile >> table[theSeed][i];
184 seq = int(theSeed);
185 }
186 }
187
188 void RanecuEngine::showStatus() const
189 {
190 std::cout << std::endl;
191 std::cout << "--------- Ranecu engine statu
192 std::cout << " Initial seed (index) = " <<
193 std::cout << " Current couple of seeds = "
194 << table[theSeed][0] << ", "
195 << table[theSeed][1] << std::endl;
196 std::cout << "-----------------------------
197 }
198
199 double RanecuEngine::flat()
200 {
201 const int index = seq;
202 long seed1 = table[index][0];
203 long seed2 = table[index][1];
204
205 int k1 = (int)(seed1/ecuyer_b);
206 int k2 = (int)(seed2/ecuyer_e);
207
208 seed1 = ecuyer_a*(seed1-k1*ecuyer_b)-k1*ecu
209 if (seed1 < 0) seed1 += shift1;
210 seed2 = ecuyer_d*(seed2-k2*ecuyer_e)-k2*ecu
211 if (seed2 < 0) seed2 += shift2;
212
213 table[index][0] = seed1;
214 table[index][1] = seed2;
215
216 long diff = seed1-seed2;
217
218 if (diff <= 0) diff += (shift1-1);
219 return (double)(diff*prec);
220 }
221
222 void RanecuEngine::flatArray(const int size, d
223 {
224 const int index = seq;
225 long seed1 = table[index][0];
226 long seed2 = table[index][1];
227 int k1, k2;
228 int i;
229
230 for (i=0; i<size; ++i)
231 {
232 k1 = (int)(seed1/ecuyer_b);
233 k2 = (int)(seed2/ecuyer_e);
234
235 seed1 = ecuyer_a*(seed1-k1*ecuyer_b)-k1*e
236 if (seed1 < 0) seed1 += shift1;
237 seed2 = ecuyer_d*(seed2-k2*ecuyer_e)-k2*e
238 if (seed2 < 0) seed2 += shift2;
239
240 long diff = seed1-seed2;
241 if (diff <= 0) diff += (shift1-1);
242
243 vect[i] = (double)(diff*prec);
244 }
245 table[index][0] = seed1;
246 table[index][1] = seed2;
247 }
248
249 RanecuEngine::operator double() {
250 return flat();
251 }
252
253 RanecuEngine::operator float() {
254 return float( flat() );
255 }
256
257 RanecuEngine::operator unsigned int() {
258 const int index = seq;
259 long seed1 = table[index][0];
260 long seed2 = table[index][1];
261
262 int k1 = (int)(seed1/ecuyer_b);
263 int k2 = (int)(seed2/ecuyer_e);
264
265 seed1 = ecuyer_a*(seed1-k1*ecuyer_b)-k1*ecu
266 if (seed1 < 0) seed1 += shift1;
267 seed2 = ecuyer_d*(seed2-k2*ecuyer_e)-k2*ecu
268 if (seed2 < 0) seed2 += shift2;
269
270 table[index][0] = seed1;
271 table[index][1] = seed2;
272 long diff = seed1-seed2;
273 if( diff <= 0 ) diff += (shift1-1);
274
275 return ((diff << 1) | (seed1&1))& 0xfffffff
276 }
277
278 std::ostream & RanecuEngine::put( std::ostream
279 {
280 char beginMarker[] = "RanecuEngine-begin";
281 os << beginMarker << "\nUvec\n";
282 std::vector<unsigned long> v = put();
283 for (unsigned int i=0; i<v.size(); ++i) {
284 os << v[i] << "\n";
285 }
286 return os;
287 }
288
289 std::vector<unsigned long> RanecuEngine::put (
290 std::vector<unsigned long> v;
291 v.push_back (engineIDulong<RanecuEngine>());
292 v.push_back(static_cast<unsigned long>(theSe
293 v.push_back(static_cast<unsigned long>(table
294 v.push_back(static_cast<unsigned long>(table
295 return v;
296 }
297
298 std::istream & RanecuEngine::get ( std::istrea
299 {
300 char beginMarker [MarkerLen];
301
302 is >> std::ws;
303 is.width(MarkerLen); // causes the next rea
304 // that many bytes, INCLUDING A TERMINAT
305 // (Stroustrup, section 21.3.2)
306 is >> beginMarker;
307 if (strcmp(beginMarker,"RanecuEngine-begin")
308 is.clear(std::ios::badbit | is.rdstate())
309 std::cerr << "\nInput stream mispositione
310 << "\nRanecuEngine state description
311 << "\nwrong engine type found." << st
312 return is;
313 }
314 return getState(is);
315 }
316
317 std::string RanecuEngine::beginTag ( ) {
318 return "RanecuEngine-begin";
319 }
320
321 std::istream & RanecuEngine::getState ( std::i
322 {
323 if ( possibleKeywordInput ( is, "Uvec", theS
324 std::vector<unsigned long> v;
325 unsigned long uu;
326 for (unsigned int ivec=0; ivec < VECTOR_ST
327 is >> uu;
328 if (!is) {
329 is.clear(std::ios::badbit | is.rdstate
330 std::cerr << "\nRanecuEngine state (ve
331 << "\ngetState() has failed."
332 << "\nInput stream is probably mispos
333 return is;
334 }
335 v.push_back(uu);
336 }
337 getState(v);
338 return (is);
339 }
340
341 // is >> theSeed; Removed, encompassed by po
342 char endMarker [MarkerLen];
343 for (int i=0; i<2; ++i) {
344 is >> table[theSeed][i];
345 }
346 is >> std::ws;
347 is.width(MarkerLen);
348 is >> endMarker;
349 if (strcmp(endMarker,"RanecuEngine-end")) {
350 is.clear(std::ios::badbit | is.rdstate())
351 std::cerr << "\nRanecuEngine state descri
352 << "\nInput stream is probably mispos
353 return is;
354 }
355
356 seq = int(theSeed);
357 return is;
358 }
359
360 bool RanecuEngine::get (const std::vector<unsi
361 if ((v[0] & 0xffffffffUL) != engineIDulong<R
362 std::cerr <<
363 "\nRanecuEngine get:state vector has wro
364 return false;
365 }
366 return getState(v);
367 }
368
369 bool RanecuEngine::getState (const std::vector
370 if (v.size() != VECTOR_STATE_SIZE ) {
371 std::cerr <<
372 "\nRanecuEngine get:state vector has wro
373 return false;
374 }
375 theSeed = v[1];
376 table[theSeed][0] = v[2];
377 table[theSeed][1] = v[3];
378 seq = int(theSeed);
379 return true;
380 }
381
382
383 } // namespace CLHEP
384