| Geant4 Cross Reference |
1 // -*- C++ -*- 1 // -*- C++ -*-
2 // 2 //
3 // ------------------------------------------- 3 // -----------------------------------------------------------------------
4 // HEP Random 4 // HEP Random
5 // --- Ranlux64Engine -- 5 // --- Ranlux64Engine ---
6 // class implementation f 6 // class implementation file
7 // ------------------------------------------- 7 // -----------------------------------------------------------------------
8 // A double-precision implementation of the Ra 8 // A double-precision implementation of the RanluxEngine generator as
9 // decsribed by the notes of the original ranl 9 // decsribed by the notes of the original ranlux author (Martin Luscher)
10 // 10 //
11 // See the note by Martin Luscher, December 19 11 // See the note by Martin Luscher, December 1997, entitiled
12 // Double-precision implementation of the rand 12 // Double-precision implementation of the random number generator ranlux
13 // 13 //
14 // =========================================== 14 // =======================================================================
15 // Ken Smith - Initial draft: 14th Jul 19 15 // Ken Smith - Initial draft: 14th Jul 1998
16 // - Removed pow() from flat me 16 // - Removed pow() from flat method 14th Jul 1998
17 // - Added conversion operators 17 // - Added conversion operators: 6th Aug 1998
18 // 18 //
19 // Mark Fischler The following were modified 19 // Mark Fischler The following were modified mostly to make the routine
20 // exactly match the Luscher algorithm in 20 // exactly match the Luscher algorithm in generating 48-bit
21 // randoms: 21 // randoms:
22 // 9/9/98 - Substantial changes in what used 22 // 9/9/98 - Substantial changes in what used to be flat() to match
23 // algorithm in Luscher's ranlxd.c 23 // algorithm in Luscher's ranlxd.c
24 // - Added update() method for 12 numbers 24 // - Added update() method for 12 numbers, making flat() trivial
25 // - Added advance() method to hold the u 25 // - Added advance() method to hold the unrolled loop for update
26 // - Distinction between three forms of s 26 // - Distinction between three forms of seeding such that it
27 // is impossible to get same sequence f 27 // is impossible to get same sequence from different forms -
28 // done by discarding some fraction of 28 // done by discarding some fraction of one macro cycle which
29 // is different for the three cases 29 // is different for the three cases
30 // - Change the misnomer "seed_table" to 30 // - Change the misnomer "seed_table" to the more accurate
31 // "randoms" 31 // "randoms"
32 // - Removed the no longer needed count12 32 // - Removed the no longer needed count12, i_lag, j_lag, etc.
33 // - Corrected seed procedure which had b 33 // - Corrected seed procedure which had been filling bits past
34 // 2^-48. This actually was very bad, 34 // 2^-48. This actually was very bad, invalidating the
35 // number theory behind the proof that 35 // number theory behind the proof that ranlxd is good.
36 // - Addition of 2**(-49) to generated nu 36 // - Addition of 2**(-49) to generated number to prevent zero
37 // from being returned; this does not a 37 // from being returned; this does not affect the sequence
38 // itself. 38 // itself.
39 // - Corrected ecu seeding, which had bee 39 // - Corrected ecu seeding, which had been supplying only
40 // numbers less than 1/2. This is prob 40 // numbers less than 1/2. This is probably moot.
41 // 9/15/98 - Modified use of the various ex 41 // 9/15/98 - Modified use of the various exponents of 2
42 // to avoid per-instance spac 42 // to avoid per-instance space overhead. Note that these
43 // are initialized in setSeed, which EV 43 // are initialized in setSeed, which EVERY constructor
44 // must invoke. 44 // must invoke.
45 // J. Marraffino - Remove dependence on hepSt 45 // J. Marraffino - Remove dependence on hepString class 13 May 1999
46 // M. Fischler - In restore, checkFile for 46 // M. Fischler - In restore, checkFile for file not found 03 Dec 2004
47 // M. Fischler - put get Methods for distri 47 // M. Fischler - put get Methods for distrib instance save/restore 12/8/04
48 // M. Fischler - split get() into tag valid 48 // M. Fischler - split get() into tag validation and
49 // getState() for anonymous r 49 // getState() for anonymous restores 12/27/04
50 // M. Fischler - put/get for vectors of ulo 50 // M. Fischler - put/get for vectors of ulongs 3/14/05
51 // M. Fischler - State-saving using only in 51 // M. Fischler - State-saving using only ints, for portability 4/12/05
52 // 52 //
53 // =========================================== 53 // =======================================================================
54 54
55 #include "CLHEP/Random/Random.h" 55 #include "CLHEP/Random/Random.h"
56 #include "CLHEP/Random/Ranlux64Engine.h" 56 #include "CLHEP/Random/Ranlux64Engine.h"
57 #include "CLHEP/Random/engineIDulong.h" 57 #include "CLHEP/Random/engineIDulong.h"
58 #include "CLHEP/Random/DoubConv.h" 58 #include "CLHEP/Random/DoubConv.h"
59 #include "CLHEP/Utility/atomic_int.h" 59 #include "CLHEP/Utility/atomic_int.h"
60 60
61 #include <atomic> << 61 #include <string.h> // for strcmp
62 #include <cstdlib> // for std::abs(int) 62 #include <cstdlib> // for std::abs(int)
63 #include <iostream> <<
64 #include <limits> // for numeric_limits 63 #include <limits> // for numeric_limits
65 #include <string.h> // for strcmp <<
66 #include <vector> <<
67 64
68 namespace CLHEP { 65 namespace CLHEP {
69 66
70 namespace { 67 namespace {
71 // Number of instances with automatic seed s 68 // Number of instances with automatic seed selection
72 CLHEP_ATOMIC_INT_TYPE numberOfEngines(0); 69 CLHEP_ATOMIC_INT_TYPE numberOfEngines(0);
73 70
74 // Maximum index into the seed table 71 // Maximum index into the seed table
75 const int maxIndex = 215; 72 const int maxIndex = 215;
76 } 73 }
77 74
78 static const int MarkerLen = 64; // Enough roo 75 static const int MarkerLen = 64; // Enough room to hold a begin or end marker.
79 76
80 77
81 #ifndef WIN32 78 #ifndef WIN32
82 namespace detail { 79 namespace detail {
83 80
84 template< std::size_t n, 81 template< std::size_t n,
85 bool = n < std::size_t(std::numeric_ 82 bool = n < std::size_t(std::numeric_limits<unsigned long>::digits) >
86 struct do_right_shift; 83 struct do_right_shift;
87 template< std::size_t n > 84 template< std::size_t n >
88 struct do_right_shift<n,true> 85 struct do_right_shift<n,true>
89 { 86 {
90 unsigned long operator()(unsigned long value 87 unsigned long operator()(unsigned long value) { return value >> n; }
91 }; 88 };
92 template< std::size_t n > 89 template< std::size_t n >
93 struct do_right_shift<n,false> 90 struct do_right_shift<n,false>
94 { 91 {
95 unsigned long operator()(unsigned long) { re 92 unsigned long operator()(unsigned long) { return 0ul; }
96 }; 93 };
97 94
98 template< std::size_t nbits > 95 template< std::size_t nbits >
99 unsigned long rshift( unsigned long value ) 96 unsigned long rshift( unsigned long value )
100 { return do_right_shift<nbits>()(value); } 97 { return do_right_shift<nbits>()(value); }
101 98
102 } // namespace detail 99 } // namespace detail
103 #endif 100 #endif
104 101
105 std::string Ranlux64Engine::name() const {retu 102 std::string Ranlux64Engine::name() const {return "Ranlux64Engine";}
106 103
107 Ranlux64Engine::Ranlux64Engine() 104 Ranlux64Engine::Ranlux64Engine()
108 : HepRandomEngine() 105 : HepRandomEngine()
109 { 106 {
110 luxury = 1; 107 luxury = 1;
111 int numEngines = numberOfEngines++; 108 int numEngines = numberOfEngines++;
112 int cycle = std::abs(int(numEngines/maxI 109 int cycle = std::abs(int(numEngines/maxIndex));
113 int curIndex = std::abs(int(numEngines%maxI 110 int curIndex = std::abs(int(numEngines%maxIndex));
114 111
115 long mask = ((cycle & 0x007fffff) << 8); 112 long mask = ((cycle & 0x007fffff) << 8);
116 long seedlist[2]; 113 long seedlist[2];
117 HepRandom::getTheTableSeeds( seedlist, curI 114 HepRandom::getTheTableSeeds( seedlist, curIndex );
118 seedlist[0] ^= mask; 115 seedlist[0] ^= mask;
119 seedlist[1] = 0; 116 seedlist[1] = 0;
120 117
121 setSeeds(seedlist, luxury); 118 setSeeds(seedlist, luxury);
122 advance ( 8 ); // Discard some iteratio 119 advance ( 8 ); // Discard some iterations and ensure that
123 // this sequence won't match one where 120 // this sequence won't match one where seeds
124 // were provided. 121 // were provided.
125 } 122 }
126 123
127 Ranlux64Engine::Ranlux64Engine(long seed, int 124 Ranlux64Engine::Ranlux64Engine(long seed, int lux)
128 : HepRandomEngine() 125 : HepRandomEngine()
129 { 126 {
130 luxury = lux; 127 luxury = lux;
131 long seedlist[2]={seed,0}; 128 long seedlist[2]={seed,0};
132 setSeeds(seedlist, lux); 129 setSeeds(seedlist, lux);
133 advance ( 2*lux + 1 ); // Discard some it 130 advance ( 2*lux + 1 ); // Discard some iterations to use a different
134 // point in the sequence. 131 // point in the sequence.
135 } 132 }
136 133
137 Ranlux64Engine::Ranlux64Engine(int rowIndex, i 134 Ranlux64Engine::Ranlux64Engine(int rowIndex, int, int lux)
138 : HepRandomEngine() 135 : HepRandomEngine()
139 { 136 {
140 luxury = lux; 137 luxury = lux;
141 int cycle = std::abs(int(rowIndex/maxIndex) 138 int cycle = std::abs(int(rowIndex/maxIndex));
142 int row = std::abs(int(rowIndex%maxIndex) 139 int row = std::abs(int(rowIndex%maxIndex));
143 long mask = (( cycle & 0x000007ff ) << 20 ) 140 long mask = (( cycle & 0x000007ff ) << 20 );
144 long seedlist[2]; 141 long seedlist[2];
145 HepRandom::getTheTableSeeds( seedlist, row 142 HepRandom::getTheTableSeeds( seedlist, row );
146 seedlist[0] ^= mask; 143 seedlist[0] ^= mask;
147 seedlist[1]= 0; 144 seedlist[1]= 0;
148 setSeeds(seedlist, lux); 145 setSeeds(seedlist, lux);
149 } 146 }
150 147
151 Ranlux64Engine::Ranlux64Engine( std::istream& 148 Ranlux64Engine::Ranlux64Engine( std::istream& is )
152 : HepRandomEngine() 149 : HepRandomEngine()
153 { 150 {
154 is >> *this; 151 is >> *this;
155 } 152 }
156 153
157 Ranlux64Engine::~Ranlux64Engine() {} 154 Ranlux64Engine::~Ranlux64Engine() {}
158 155
159 double Ranlux64Engine::flat() { 156 double Ranlux64Engine::flat() {
160 // Luscher improves the speed by computing s 157 // Luscher improves the speed by computing several numbers in a shot,
161 // in a manner similar to that of the Tauswo 158 // in a manner similar to that of the Tausworth in DualRand or the Hurd
162 // engines. Thus, the real work is done in 159 // engines. Thus, the real work is done in update(). Here we merely ensure
163 // that zero, which the algorithm can produc 160 // that zero, which the algorithm can produce, is never returned by flat().
164 161
165 if (index <= 0) update(); 162 if (index <= 0) update();
166 return randoms[--index] + twoToMinus_49(); 163 return randoms[--index] + twoToMinus_49();
167 } 164 }
168 165
169 void Ranlux64Engine::update() { 166 void Ranlux64Engine::update() {
170 // Update the stash of twelve random numbers 167 // Update the stash of twelve random numbers.
171 // When this routione is entered, index is a 168 // When this routione is entered, index is always 0. The randoms
172 // contains the last 12 numbers in the seque 169 // contains the last 12 numbers in the sequents: s[0] is x[a+11],
173 // s[1] is x[a+10] ... and s[11] is x[a] for 170 // s[1] is x[a+10] ... and s[11] is x[a] for some a. Carry contains
174 // the last carry value (c[a+11]). 171 // the last carry value (c[a+11]).
175 // 172 //
176 // The recursion relation (3) in Luscher's n 173 // The recursion relation (3) in Luscher's note says
177 // delta[n] = x[n-s] = x[n-r] -c[n-1] or f 174 // delta[n] = x[n-s] = x[n-r] -c[n-1] or for n=a+12,
178 // delta[a+12] = x[a+7] - x[a] -c[a+11] wh 175 // delta[a+12] = x[a+7] - x[a] -c[a+11] where we use r=12, s=5 per eqn. (7)
179 // This reduces to 176 // This reduces to
180 // s[11] = s[4] - s[11] - carry. 177 // s[11] = s[4] - s[11] - carry.
181 // The next number similarly will be given b 178 // The next number similarly will be given by s[10] = s[3] - s[10] - carry,
182 // and so forth until s[0] is filled. 179 // and so forth until s[0] is filled.
183 // 180 //
184 // However, we need to skip 397, 202 or 109 181 // However, we need to skip 397, 202 or 109 numbers - these are not divisible
185 // by 12 - to "fare well in the spectral tes 182 // by 12 - to "fare well in the spectral test".
186 183
187 advance(pDozens); 184 advance(pDozens);
188 185
189 // Since we wish at the end to have the 12 l 186 // Since we wish at the end to have the 12 last numbers in the order of
190 // s[11] first, till s[0] last, we will have 187 // s[11] first, till s[0] last, we will have to do 1, 10, or 1 iterations
191 // and then re-arrange to place to get the o 188 // and then re-arrange to place to get the oldest one in s[11].
192 // Generically, this will imply re-arranging 189 // Generically, this will imply re-arranging the s array at the end,
193 // but we can treat the special case of endI 190 // but we can treat the special case of endIters = 1 separately for superior
194 // efficiency in the cases of levels 0 and 2 191 // efficiency in the cases of levels 0 and 2.
195 192
196 double y1; 193 double y1;
197 194
198 if ( endIters == 1 ) { // Luxury levels 0 195 if ( endIters == 1 ) { // Luxury levels 0 and 2 will go here
199 y1 = randoms[ 4] - randoms[11] - carry; 196 y1 = randoms[ 4] - randoms[11] - carry;
200 if ( y1 < 0.0 ) { 197 if ( y1 < 0.0 ) {
201 y1 += 1.0; 198 y1 += 1.0;
202 carry = twoToMinus_48(); 199 carry = twoToMinus_48();
203 } else { 200 } else {
204 carry = 0.0; 201 carry = 0.0;
205 } 202 }
206 randoms[11] = randoms[10]; 203 randoms[11] = randoms[10];
207 randoms[10] = randoms[ 9]; 204 randoms[10] = randoms[ 9];
208 randoms[ 9] = randoms[ 8]; 205 randoms[ 9] = randoms[ 8];
209 randoms[ 8] = randoms[ 7]; 206 randoms[ 8] = randoms[ 7];
210 randoms[ 7] = randoms[ 6]; 207 randoms[ 7] = randoms[ 6];
211 randoms[ 6] = randoms[ 5]; 208 randoms[ 6] = randoms[ 5];
212 randoms[ 5] = randoms[ 4]; 209 randoms[ 5] = randoms[ 4];
213 randoms[ 4] = randoms[ 3]; 210 randoms[ 4] = randoms[ 3];
214 randoms[ 3] = randoms[ 2]; 211 randoms[ 3] = randoms[ 2];
215 randoms[ 2] = randoms[ 1]; 212 randoms[ 2] = randoms[ 1];
216 randoms[ 1] = randoms[ 0]; 213 randoms[ 1] = randoms[ 0];
217 randoms[ 0] = y1; 214 randoms[ 0] = y1;
218 215
219 } else { 216 } else {
220 217
221 int m, nr, ns; 218 int m, nr, ns;
222 for ( m = 0, nr = 11, ns = 4; m < endIters 219 for ( m = 0, nr = 11, ns = 4; m < endIters; ++m, --nr ) {
223 y1 = randoms [ns] - randoms[nr] - carry; 220 y1 = randoms [ns] - randoms[nr] - carry;
224 if ( y1 < 0.0 ) { 221 if ( y1 < 0.0 ) {
225 y1 += 1.0; 222 y1 += 1.0;
226 carry = twoToMinus_48(); 223 carry = twoToMinus_48();
227 } else { 224 } else {
228 carry = 0.0; 225 carry = 0.0;
229 } 226 }
230 randoms[nr] = y1; 227 randoms[nr] = y1;
231 --ns; 228 --ns;
232 if ( ns < 0 ) { 229 if ( ns < 0 ) {
233 ns = 11; 230 ns = 11;
234 } 231 }
235 } // loop on m 232 } // loop on m
236 233
237 double temp[12]; 234 double temp[12];
238 for (m=0; m<12; m++) { 235 for (m=0; m<12; m++) {
239 temp[m]=randoms[m]; 236 temp[m]=randoms[m];
240 } 237 }
241 238
242 ns = 11 - endIters; 239 ns = 11 - endIters;
243 for (m=11; m>=0; --m) { 240 for (m=11; m>=0; --m) {
244 randoms[m] = temp[ns]; 241 randoms[m] = temp[ns];
245 --ns; 242 --ns;
246 if ( ns < 0 ) { 243 if ( ns < 0 ) {
247 ns = 11; 244 ns = 11;
248 } 245 }
249 } 246 }
250 247
251 } 248 }
252 249
253 // Now when we return, there are 12 fresh us 250 // Now when we return, there are 12 fresh usable numbers in s[11] ... s[0]
254 251
255 index = 12; << 252 index = 11;
256 253
257 } // update() 254 } // update()
258 255
259 void Ranlux64Engine::advance(int dozens) { 256 void Ranlux64Engine::advance(int dozens) {
260 257
261 double y1, y2, y3; 258 double y1, y2, y3;
262 double cValue = twoToMinus_48(); 259 double cValue = twoToMinus_48();
263 double zero = 0.0; 260 double zero = 0.0;
264 double one = 1.0; 261 double one = 1.0;
265 262
266 // Technical note: We use Luscher's trick 263 // Technical note: We use Luscher's trick to only do the
267 // carry subtraction when we really have t 264 // carry subtraction when we really have to. Like him, we use
268 // three registers instead of two so that 265 // three registers instead of two so that we avoid sequences
269 // like storing y1 then immediately replac 266 // like storing y1 then immediately replacing its value:
270 // some architectures lose time when this 267 // some architectures lose time when this is done.
271 268
272 // Luscher's ranlxd.c fills the stash go 269 // Luscher's ranlxd.c fills the stash going
273 // upward. We fill it downward to save a 270 // upward. We fill it downward to save a bit of time in the
274 // flat() routine at no cost later. This 271 // flat() routine at no cost later. This means that while
275 // Luscher's ir is jr+5, our n-r is (n-s)- 272 // Luscher's ir is jr+5, our n-r is (n-s)-5. (Note that
276 // though ranlxd.c initializes ir and jr t 273 // though ranlxd.c initializes ir and jr to 11 and 7, ir as
277 // used is 5 more than jr because update i 274 // used is 5 more than jr because update is entered after
278 // incrementing ir.) 275 // incrementing ir.)
279 // 276 //
280 277
281 // I have CAREFULLY checked that the algor 278 // I have CAREFULLY checked that the algorithms do match
282 // in all details. 279 // in all details.
283 280
284 int k; 281 int k;
285 for ( k = dozens; k > 0; --k ) { 282 for ( k = dozens; k > 0; --k ) {
286 283
287 y1 = randoms[ 4] - randoms[11] - carry; 284 y1 = randoms[ 4] - randoms[11] - carry;
288 285
289 y2 = randoms[ 3] - randoms[10]; 286 y2 = randoms[ 3] - randoms[10];
290 if ( y1 < zero ) { 287 if ( y1 < zero ) {
291 y1 += one; 288 y1 += one;
292 y2 -= cValue; 289 y2 -= cValue;
293 } 290 }
294 randoms[11] = y1; 291 randoms[11] = y1;
295 292
296 y3 = randoms[ 2] - randoms[ 9]; 293 y3 = randoms[ 2] - randoms[ 9];
297 if ( y2 < zero ) { 294 if ( y2 < zero ) {
298 y2 += one; 295 y2 += one;
299 y3 -= cValue; 296 y3 -= cValue;
300 } 297 }
301 randoms[10] = y2; 298 randoms[10] = y2;
302 299
303 y1 = randoms[ 1] - randoms[ 8]; 300 y1 = randoms[ 1] - randoms[ 8];
304 if ( y3 < zero ) { 301 if ( y3 < zero ) {
305 y3 += one; 302 y3 += one;
306 y1 -= cValue; 303 y1 -= cValue;
307 } 304 }
308 randoms[ 9] = y3; 305 randoms[ 9] = y3;
309 306
310 y2 = randoms[ 0] - randoms[ 7]; 307 y2 = randoms[ 0] - randoms[ 7];
311 if ( y1 < zero ) { 308 if ( y1 < zero ) {
312 y1 += one; 309 y1 += one;
313 y2 -= cValue; 310 y2 -= cValue;
314 } 311 }
315 randoms[ 8] = y1; 312 randoms[ 8] = y1;
316 313
317 y3 = randoms[11] - randoms[ 6]; 314 y3 = randoms[11] - randoms[ 6];
318 if ( y2 < zero ) { 315 if ( y2 < zero ) {
319 y2 += one; 316 y2 += one;
320 y3 -= cValue; 317 y3 -= cValue;
321 } 318 }
322 randoms[ 7] = y2; 319 randoms[ 7] = y2;
323 320
324 y1 = randoms[10] - randoms[ 5]; 321 y1 = randoms[10] - randoms[ 5];
325 if ( y3 < zero ) { 322 if ( y3 < zero ) {
326 y3 += one; 323 y3 += one;
327 y1 -= cValue; 324 y1 -= cValue;
328 } 325 }
329 randoms[ 6] = y3; 326 randoms[ 6] = y3;
330 327
331 y2 = randoms[ 9] - randoms[ 4]; 328 y2 = randoms[ 9] - randoms[ 4];
332 if ( y1 < zero ) { 329 if ( y1 < zero ) {
333 y1 += one; 330 y1 += one;
334 y2 -= cValue; 331 y2 -= cValue;
335 } 332 }
336 randoms[ 5] = y1; 333 randoms[ 5] = y1;
337 334
338 y3 = randoms[ 8] - randoms[ 3]; 335 y3 = randoms[ 8] - randoms[ 3];
339 if ( y2 < zero ) { 336 if ( y2 < zero ) {
340 y2 += one; 337 y2 += one;
341 y3 -= cValue; 338 y3 -= cValue;
342 } 339 }
343 randoms[ 4] = y2; 340 randoms[ 4] = y2;
344 341
345 y1 = randoms[ 7] - randoms[ 2]; 342 y1 = randoms[ 7] - randoms[ 2];
346 if ( y3 < zero ) { 343 if ( y3 < zero ) {
347 y3 += one; 344 y3 += one;
348 y1 -= cValue; 345 y1 -= cValue;
349 } 346 }
350 randoms[ 3] = y3; 347 randoms[ 3] = y3;
351 348
352 y2 = randoms[ 6] - randoms[ 1]; 349 y2 = randoms[ 6] - randoms[ 1];
353 if ( y1 < zero ) { 350 if ( y1 < zero ) {
354 y1 += one; 351 y1 += one;
355 y2 -= cValue; 352 y2 -= cValue;
356 } 353 }
357 randoms[ 2] = y1; 354 randoms[ 2] = y1;
358 355
359 y3 = randoms[ 5] - randoms[ 0]; 356 y3 = randoms[ 5] - randoms[ 0];
360 if ( y2 < zero ) { 357 if ( y2 < zero ) {
361 y2 += one; 358 y2 += one;
362 y3 -= cValue; 359 y3 -= cValue;
363 } 360 }
364 randoms[ 1] = y2; 361 randoms[ 1] = y2;
365 362
366 if ( y3 < zero ) { 363 if ( y3 < zero ) {
367 y3 += one; 364 y3 += one;
368 carry = cValue; 365 carry = cValue;
369 } 366 }
370 randoms[ 0] = y3; 367 randoms[ 0] = y3;
371 368
372 } // End of major k loop doing 12 numbers at 369 } // End of major k loop doing 12 numbers at each cycle
373 370
374 } // advance(dozens) 371 } // advance(dozens)
375 372
376 void Ranlux64Engine::flatArray(const int size, 373 void Ranlux64Engine::flatArray(const int size, double* vect) {
377 for( int i=0; i < size; ++i ) { 374 for( int i=0; i < size; ++i ) {
378 vect[i] = flat(); 375 vect[i] = flat();
379 } 376 }
380 } 377 }
381 378
382 void Ranlux64Engine::setSeed(long seed, int lu 379 void Ranlux64Engine::setSeed(long seed, int lux) {
383 380
384 // The initialization is carried out using a M 381 // The initialization is carried out using a Multiplicative
385 // Congruential generator using formula consta 382 // Congruential generator using formula constants of L'Ecuyer
386 // as described in "A review of pseudorandom n 383 // as described in "A review of pseudorandom number generators"
387 // (Fred James) published in Computer Physics 384 // (Fred James) published in Computer Physics Communications 60 (1990)
388 // pages 329-344 385 // pages 329-344
389 386
390 const int ecuyer_a(53668); 387 const int ecuyer_a(53668);
391 const int ecuyer_b(40014); 388 const int ecuyer_b(40014);
392 const int ecuyer_c(12211); 389 const int ecuyer_c(12211);
393 const int ecuyer_d(2147483563); 390 const int ecuyer_d(2147483563);
394 391
395 const int lux_levels[3] = {109, 202, 397}; 392 const int lux_levels[3] = {109, 202, 397};
396 theSeed = seed; 393 theSeed = seed;
397 394
398 if( (lux > 2)||(lux < 0) ){ 395 if( (lux > 2)||(lux < 0) ){
399 pDiscard = (lux >= 12) ? (lux-12) : lux_l 396 pDiscard = (lux >= 12) ? (lux-12) : lux_levels[1];
400 }else{ 397 }else{
401 pDiscard = lux_levels[luxury]; 398 pDiscard = lux_levels[luxury];
402 } 399 }
403 pDozens = pDiscard / 12; 400 pDozens = pDiscard / 12;
404 endIters = pDiscard % 12; 401 endIters = pDiscard % 12;
405 402
406 long init_table[24]; 403 long init_table[24];
407 long next_seed = seed; 404 long next_seed = seed;
408 long k_multiple; 405 long k_multiple;
409 int i; 406 int i;
410 next_seed &= 0xffffffff; 407 next_seed &= 0xffffffff;
411 while( next_seed >= ecuyer_d ) { 408 while( next_seed >= ecuyer_d ) {
412 next_seed -= ecuyer_d; 409 next_seed -= ecuyer_d;
413 } 410 }
414 411
415 for(i = 0;i != 24;i++){ 412 for(i = 0;i != 24;i++){
416 k_multiple = next_seed / ecuyer_a; 413 k_multiple = next_seed / ecuyer_a;
417 next_seed = ecuyer_b * (next_seed - k_mul 414 next_seed = ecuyer_b * (next_seed - k_multiple * ecuyer_a)
418 - k_mul 415 - k_multiple * ecuyer_c;
419 if(next_seed < 0) { 416 if(next_seed < 0) {
420 next_seed += ecuyer_d; 417 next_seed += ecuyer_d;
421 } 418 }
422 next_seed &= 0xffffffff; 419 next_seed &= 0xffffffff;
423 init_table[i] = next_seed; 420 init_table[i] = next_seed;
424 } 421 }
425 // are we on a 64bit machine? 422 // are we on a 64bit machine?
426 if( sizeof(long) >= 8 ) { 423 if( sizeof(long) >= 8 ) {
427 int64_t topbits1, topbits2; 424 int64_t topbits1, topbits2;
428 #ifdef WIN32 425 #ifdef WIN32
429 topbits1 = ( (int64_t) seed >> 32) & 0xff 426 topbits1 = ( (int64_t) seed >> 32) & 0xffff ;
430 topbits2 = ( (int64_t) seed >> 48) & 0xff 427 topbits2 = ( (int64_t) seed >> 48) & 0xffff ;
431 #else 428 #else
432 topbits1 = detail::rshift<32>(seed) & 0xf 429 topbits1 = detail::rshift<32>(seed) & 0xffff ;
433 topbits2 = detail::rshift<48>(seed) & 0xf 430 topbits2 = detail::rshift<48>(seed) & 0xffff ;
434 #endif 431 #endif
435 init_table[0] ^= topbits1; 432 init_table[0] ^= topbits1;
436 init_table[2] ^= topbits2; 433 init_table[2] ^= topbits2;
437 //std::cout << " init_table[0] " << init_ 434 //std::cout << " init_table[0] " << init_table[0] << " from " << topbits1 << std::endl;
438 //std::cout << " init_table[2] " << init_ 435 //std::cout << " init_table[2] " << init_table[2] << " from " << topbits2 << std::endl;
439 } 436 }
440 437
441 for(i = 0;i < 12; i++){ 438 for(i = 0;i < 12; i++){
442 randoms[i] = (init_table[2*i ] ) * 439 randoms[i] = (init_table[2*i ] ) * 2.0 * twoToMinus_32() +
443 (init_table[2*i+1] >> 15) * 440 (init_table[2*i+1] >> 15) * twoToMinus_48();
444 //if( randoms[i] < 0. || randoms[i] > 1. 441 //if( randoms[i] < 0. || randoms[i] > 1. ) {
445 //std::cout << "setSeed: init_table " << 442 //std::cout << "setSeed: init_table " << init_table[2*i ] << std::endl;
446 //std::cout << "setSeed: init_table " << 443 //std::cout << "setSeed: init_table " << init_table[2*i+1] << std::endl;
447 //std::cout << "setSeed: random " << i < 444 //std::cout << "setSeed: random " << i << " is " << randoms[i] << std::endl;
448 //} 445 //}
449 } 446 }
450 447
451 carry = 0.0; 448 carry = 0.0;
452 if ( randoms[11] == 0. ) carry = twoToMinus_ 449 if ( randoms[11] == 0. ) carry = twoToMinus_48();
453 // Perform an update before returning the fi << 450 index = 11;
454 index = -1; <<
455 451
456 } // setSeed() 452 } // setSeed()
457 453
458 void Ranlux64Engine::setSeeds(const long * see 454 void Ranlux64Engine::setSeeds(const long * seeds, int lux) {
459 // old code only uses the first long in seeds 455 // old code only uses the first long in seeds
460 // setSeed( *seeds ? *seeds : 32767, lux ); 456 // setSeed( *seeds ? *seeds : 32767, lux );
461 // theSeeds = seeds; 457 // theSeeds = seeds;
462 458
463 // using code from Ranlux - even those are 32b 459 // using code from Ranlux - even those are 32bit seeds,
464 // that is good enough to completely different 460 // that is good enough to completely differentiate the sequences
465 461
466 const int ecuyer_a = 53668; 462 const int ecuyer_a = 53668;
467 const int ecuyer_b = 40014; 463 const int ecuyer_b = 40014;
468 const int ecuyer_c = 12211; 464 const int ecuyer_c = 12211;
469 const int ecuyer_d = 2147483563; 465 const int ecuyer_d = 2147483563;
470 466
471 const int lux_levels[3] = {109, 202, 397}; 467 const int lux_levels[3] = {109, 202, 397};
472 const long *seedptr; 468 const long *seedptr;
473 469
474 theSeeds = seeds; 470 theSeeds = seeds;
475 seedptr = seeds; 471 seedptr = seeds;
476 472
477 if(seeds == 0){ 473 if(seeds == 0){
478 setSeed(theSeed,lux); 474 setSeed(theSeed,lux);
479 theSeeds = &theSeed; 475 theSeeds = &theSeed;
480 return; 476 return;
481 } 477 }
482 478
483 theSeed = *seeds; 479 theSeed = *seeds;
484 480
485 // number of additional random numbers that ne 481 // number of additional random numbers that need to be 'thrown away'
486 // every 24 numbers is set using luxury level 482 // every 24 numbers is set using luxury level variable.
487 483
488 if( (lux > 2)||(lux < 0) ){ 484 if( (lux > 2)||(lux < 0) ){
489 pDiscard = (lux >= 12) ? (lux-12) : lux_l 485 pDiscard = (lux >= 12) ? (lux-12) : lux_levels[1];
490 }else{ 486 }else{
491 pDiscard = lux_levels[luxury]; 487 pDiscard = lux_levels[luxury];
492 } 488 }
493 pDozens = pDiscard / 12; 489 pDozens = pDiscard / 12;
494 endIters = pDiscard % 12; 490 endIters = pDiscard % 12;
495 491
496 long init_table[24]; 492 long init_table[24];
497 long next_seed = *seeds; 493 long next_seed = *seeds;
498 long k_multiple; 494 long k_multiple;
499 int i; 495 int i;
500 496
501 for( i = 0;(i != 24)&&(*seedptr != 0);i++){ 497 for( i = 0;(i != 24)&&(*seedptr != 0);i++){
502 init_table[i] = *seedptr & 0xffffffff; 498 init_table[i] = *seedptr & 0xffffffff;
503 seedptr++; 499 seedptr++;
504 } 500 }
505 501
506 if(i != 24){ 502 if(i != 24){
507 next_seed = init_table[i-1]; 503 next_seed = init_table[i-1];
508 for(;i != 24;i++){ 504 for(;i != 24;i++){
509 k_multiple = next_seed / ecuyer_a; 505 k_multiple = next_seed / ecuyer_a;
510 next_seed = ecuyer_b * (next_seed - k_multip 506 next_seed = ecuyer_b * (next_seed - k_multiple * ecuyer_a)
511 - k_multip 507 - k_multiple * ecuyer_c;
512 if(next_seed < 0) { 508 if(next_seed < 0) {
513 next_seed += ecuyer_d; 509 next_seed += ecuyer_d;
514 } 510 }
515 next_seed &= 0xffffffff; 511 next_seed &= 0xffffffff;
516 init_table[i] = next_seed; 512 init_table[i] = next_seed;
517 } 513 }
518 } 514 }
519 515
520 for(i = 0;i < 12; i++){ 516 for(i = 0;i < 12; i++){
521 randoms[i] = (init_table[2*i ] ) * 517 randoms[i] = (init_table[2*i ] ) * 2.0 * twoToMinus_32() +
522 (init_table[2*i+1] >> 15) * 518 (init_table[2*i+1] >> 15) * twoToMinus_48();
523 } 519 }
524 520
525 carry = 0.0; 521 carry = 0.0;
526 if ( randoms[11] == 0. ) carry = twoToMinus_ 522 if ( randoms[11] == 0. ) carry = twoToMinus_48();
527 // Perform an update before returning the fi << 523 index = 11;
528 index = -1; <<
529 524
530 } 525 }
531 526
532 void Ranlux64Engine::saveStatus( const char fi 527 void Ranlux64Engine::saveStatus( const char filename[] ) const
533 { 528 {
534 std::ofstream outFile( filename, std::ios:: 529 std::ofstream outFile( filename, std::ios::out ) ;
535 if (!outFile.bad()) { 530 if (!outFile.bad()) {
536 outFile << "Uvec\n"; 531 outFile << "Uvec\n";
537 std::vector<unsigned long> v = put(); 532 std::vector<unsigned long> v = put();
538 for (unsigned int i=0; i<v.size(); ++i) { 533 for (unsigned int i=0; i<v.size(); ++i) {
539 outFile << v[i] << "\n"; 534 outFile << v[i] << "\n";
540 } 535 }
541 } 536 }
542 } 537 }
543 538
544 void Ranlux64Engine::restoreStatus( const char 539 void Ranlux64Engine::restoreStatus( const char filename[] )
545 { 540 {
546 std::ifstream inFile( filename, std::ios::i 541 std::ifstream inFile( filename, std::ios::in);
547 if (!checkFile ( inFile, filename, engineNa 542 if (!checkFile ( inFile, filename, engineName(), "restoreStatus" )) {
548 std::cerr << " -- Engine state remains u 543 std::cerr << " -- Engine state remains unchanged\n";
549 return; 544 return;
550 } 545 }
551 if ( possibleKeywordInput ( inFile, "Uvec", 546 if ( possibleKeywordInput ( inFile, "Uvec", theSeed ) ) {
552 std::vector<unsigned long> v; 547 std::vector<unsigned long> v;
553 unsigned long xin; 548 unsigned long xin;
554 for (unsigned int ivec=0; ivec < VECTOR_ST 549 for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) {
555 inFile >> xin; 550 inFile >> xin;
556 if (!inFile) { 551 if (!inFile) {
557 inFile.clear(std::ios::badbit | inFile 552 inFile.clear(std::ios::badbit | inFile.rdstate());
558 std::cerr << "\nJamesRandom state (vec 553 std::cerr << "\nJamesRandom state (vector) description improper."
559 << "\nrestoreStatus has failed." 554 << "\nrestoreStatus has failed."
560 << "\nInput stream is probably mispos 555 << "\nInput stream is probably mispositioned now." << std::endl;
561 return; 556 return;
562 } 557 }
563 v.push_back(xin); 558 v.push_back(xin);
564 } 559 }
565 getState(v); 560 getState(v);
566 return; 561 return;
567 } 562 }
568 563
569 if (!inFile.bad() && !inFile.eof()) { 564 if (!inFile.bad() && !inFile.eof()) {
570 // inFile >> theSeed; removed -- encompas 565 // inFile >> theSeed; removed -- encompased by possibleKeywordInput
571 for (int i=0; i<12; ++i) { 566 for (int i=0; i<12; ++i) {
572 inFile >> randoms[i]; 567 inFile >> randoms[i];
573 } 568 }
574 inFile >> carry; inFile >> index; 569 inFile >> carry; inFile >> index;
575 inFile >> luxury; inFile >> pDiscard; 570 inFile >> luxury; inFile >> pDiscard;
576 pDozens = pDiscard / 12; 571 pDozens = pDiscard / 12;
577 endIters = pDiscard % 12; 572 endIters = pDiscard % 12;
578 } 573 }
579 } 574 }
580 575
581 void Ranlux64Engine::showStatus() const 576 void Ranlux64Engine::showStatus() const
582 { 577 {
583 std::cout << std::endl; 578 std::cout << std::endl;
584 std::cout << "--------- Ranlux engine statu 579 std::cout << "--------- Ranlux engine status ---------" << std::endl;
585 std::cout << " Initial seed = " << theSeed 580 std::cout << " Initial seed = " << theSeed << std::endl;
586 std::cout << " randoms[] = "; 581 std::cout << " randoms[] = ";
587 for (int i=0; i<12; ++i) { 582 for (int i=0; i<12; ++i) {
588 std::cout << randoms[i] << std::endl; 583 std::cout << randoms[i] << std::endl;
589 } 584 }
590 std::cout << std::endl; 585 std::cout << std::endl;
591 std::cout << " carry = " << carry << ", ind 586 std::cout << " carry = " << carry << ", index = " << index << std::endl;
592 std::cout << " luxury = " << luxury << " pD 587 std::cout << " luxury = " << luxury << " pDiscard = "
593 << pDiscard << std::endl; 588 << pDiscard << std::endl;
594 std::cout << "----------------------------- 589 std::cout << "----------------------------------------" << std::endl;
595 } 590 }
596 591
597 std::ostream & Ranlux64Engine::put( std::ostre 592 std::ostream & Ranlux64Engine::put( std::ostream& os ) const
598 { 593 {
599 char beginMarker[] = "Ranlux64Engine-begin" 594 char beginMarker[] = "Ranlux64Engine-begin";
600 os << beginMarker << "\nUvec\n"; 595 os << beginMarker << "\nUvec\n";
601 std::vector<unsigned long> v = put(); 596 std::vector<unsigned long> v = put();
602 for (unsigned int i=0; i<v.size(); ++i) { 597 for (unsigned int i=0; i<v.size(); ++i) {
603 os << v[i] << "\n"; 598 os << v[i] << "\n";
604 } 599 }
605 return os; 600 return os;
606 } 601 }
607 602
608 std::vector<unsigned long> Ranlux64Engine::put 603 std::vector<unsigned long> Ranlux64Engine::put () const {
609 std::vector<unsigned long> v; 604 std::vector<unsigned long> v;
610 v.push_back (engineIDulong<Ranlux64Engine>() 605 v.push_back (engineIDulong<Ranlux64Engine>());
611 std::vector<unsigned long> t; 606 std::vector<unsigned long> t;
612 for (int i=0; i<12; ++i) { 607 for (int i=0; i<12; ++i) {
613 t = DoubConv::dto2longs(randoms[i]); 608 t = DoubConv::dto2longs(randoms[i]);
614 v.push_back(t[0]); v.push_back(t[1]); 609 v.push_back(t[0]); v.push_back(t[1]);
615 } 610 }
616 t = DoubConv::dto2longs(carry); 611 t = DoubConv::dto2longs(carry);
617 v.push_back(t[0]); v.push_back(t[1]); 612 v.push_back(t[0]); v.push_back(t[1]);
618 v.push_back(static_cast<unsigned long>(index 613 v.push_back(static_cast<unsigned long>(index));
619 v.push_back(static_cast<unsigned long>(luxur 614 v.push_back(static_cast<unsigned long>(luxury));
620 v.push_back(static_cast<unsigned long>(pDisc 615 v.push_back(static_cast<unsigned long>(pDiscard));
621 return v; 616 return v;
622 } 617 }
623 618
624 std::istream & Ranlux64Engine::get ( std::istr 619 std::istream & Ranlux64Engine::get ( std::istream& is )
625 { 620 {
626 char beginMarker [MarkerLen]; 621 char beginMarker [MarkerLen];
627 is >> std::ws; 622 is >> std::ws;
628 is.width(MarkerLen); // causes the next rea 623 is.width(MarkerLen); // causes the next read to the char* to be <=
629 // that many bytes, INCLUDING A TERMINAT 624 // that many bytes, INCLUDING A TERMINATION \0
630 // (Stroustrup, section 21.3.2) 625 // (Stroustrup, section 21.3.2)
631 is >> beginMarker; 626 is >> beginMarker;
632 if (strcmp(beginMarker,"Ranlux64Engine-begin 627 if (strcmp(beginMarker,"Ranlux64Engine-begin")) {
633 is.clear(std::ios::badbit | is.rdstate()) 628 is.clear(std::ios::badbit | is.rdstate());
634 std::cerr << "\nInput stream mispositione 629 std::cerr << "\nInput stream mispositioned or"
635 << "\nRanlux64Engine state descriptio 630 << "\nRanlux64Engine state description missing or"
636 << "\nwrong engine type found." << st 631 << "\nwrong engine type found." << std::endl;
637 return is; 632 return is;
638 } 633 }
639 return getState(is); 634 return getState(is);
640 } 635 }
641 636
642 std::string Ranlux64Engine::beginTag ( ) { 637 std::string Ranlux64Engine::beginTag ( ) {
643 return "Ranlux64Engine-begin"; 638 return "Ranlux64Engine-begin";
644 } 639 }
645 640
646 std::istream & Ranlux64Engine::getState ( std: 641 std::istream & Ranlux64Engine::getState ( std::istream& is )
647 { 642 {
648 if ( possibleKeywordInput ( is, "Uvec", theS 643 if ( possibleKeywordInput ( is, "Uvec", theSeed ) ) {
649 std::vector<unsigned long> v; 644 std::vector<unsigned long> v;
650 unsigned long uu; 645 unsigned long uu;
651 for (unsigned int ivec=0; ivec < VECTOR_ST 646 for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) {
652 is >> uu; 647 is >> uu;
653 if (!is) { 648 if (!is) {
654 is.clear(std::ios::badbit | is.rdstate 649 is.clear(std::ios::badbit | is.rdstate());
655 std::cerr << "\nRanlux64Engine state ( 650 std::cerr << "\nRanlux64Engine state (vector) description improper."
656 << "\ngetState() has failed." 651 << "\ngetState() has failed."
657 << "\nInput stream is probably mispos 652 << "\nInput stream is probably mispositioned now." << std::endl;
658 return is; 653 return is;
659 } 654 }
660 v.push_back(uu); 655 v.push_back(uu);
661 } 656 }
662 getState(v); 657 getState(v);
663 return (is); 658 return (is);
664 } 659 }
665 660
666 // is >> theSeed; Removed, encompassed by po 661 // is >> theSeed; Removed, encompassed by possibleKeywordInput()
667 662
668 char endMarker [MarkerLen]; 663 char endMarker [MarkerLen];
669 for (int i=0; i<12; ++i) { 664 for (int i=0; i<12; ++i) {
670 is >> randoms[i]; 665 is >> randoms[i];
671 } 666 }
672 is >> carry; is >> index; 667 is >> carry; is >> index;
673 is >> luxury; is >> pDiscard; 668 is >> luxury; is >> pDiscard;
674 pDozens = pDiscard / 12; 669 pDozens = pDiscard / 12;
675 endIters = pDiscard % 12; 670 endIters = pDiscard % 12;
676 is >> std::ws; 671 is >> std::ws;
677 is.width(MarkerLen); 672 is.width(MarkerLen);
678 is >> endMarker; 673 is >> endMarker;
679 if (strcmp(endMarker,"Ranlux64Engine-end")) 674 if (strcmp(endMarker,"Ranlux64Engine-end")) {
680 is.clear(std::ios::badbit | is.rdstate()) 675 is.clear(std::ios::badbit | is.rdstate());
681 std::cerr << "\nRanlux64Engine state desc 676 std::cerr << "\nRanlux64Engine state description incomplete."
682 << "\nInput stream is probably mispos 677 << "\nInput stream is probably mispositioned now." << std::endl;
683 return is; 678 return is;
684 } 679 }
685 return is; 680 return is;
686 } 681 }
687 682
688 bool Ranlux64Engine::get (const std::vector<un 683 bool Ranlux64Engine::get (const std::vector<unsigned long> & v) {
689 if ((v[0] & 0xffffffffUL) != engineIDulong<R 684 if ((v[0] & 0xffffffffUL) != engineIDulong<Ranlux64Engine>()) {
690 std::cerr << 685 std::cerr <<
691 "\nRanlux64Engine get:state vector has w 686 "\nRanlux64Engine get:state vector has wrong ID word - state unchanged\n";
692 return false; 687 return false;
693 } 688 }
694 return getState(v); 689 return getState(v);
695 } 690 }
696 691
697 bool Ranlux64Engine::getState (const std::vect 692 bool Ranlux64Engine::getState (const std::vector<unsigned long> & v) {
698 if (v.size() != VECTOR_STATE_SIZE ) { 693 if (v.size() != VECTOR_STATE_SIZE ) {
699 std::cerr << 694 std::cerr <<
700 "\nRanlux64Engine get:state vector has w 695 "\nRanlux64Engine get:state vector has wrong length - state unchanged\n";
701 return false; 696 return false;
702 } 697 }
703 std::vector<unsigned long> t(2); 698 std::vector<unsigned long> t(2);
704 for (int i=0; i<12; ++i) { 699 for (int i=0; i<12; ++i) {
705 t[0] = v[2*i+1]; t[1] = v[2*i+2]; 700 t[0] = v[2*i+1]; t[1] = v[2*i+2];
706 randoms[i] = DoubConv::longs2double(t); 701 randoms[i] = DoubConv::longs2double(t);
707 } 702 }
708 t[0] = v[25]; t[1] = v[26]; 703 t[0] = v[25]; t[1] = v[26];
709 carry = DoubConv::longs2double(t); 704 carry = DoubConv::longs2double(t);
710 index = (int)v[27]; << 705 index = v[27];
711 luxury = (int)v[28]; << 706 luxury = v[28];
712 pDiscard = (int)v[29]; << 707 pDiscard = v[29];
713 return true; 708 return true;
714 } 709 }
715 710
716 } // namespace CLHEP 711 } // namespace CLHEP
717 712