| Geant4 Cross Reference |
1 //
1 //
2 // -*- C++ -*-
2 // -*- C++ -*-
3 //
3 //
4 // ------------------------------------------- 4 // -----------------------------------------------------------------------
5 // HEP Random
5 // HEP Random
6 // --- MixMaxRng ---
6 // --- MixMaxRng ---
7 // class implementation fi 7 // class implementation file
8 // ------------------------------------------- 8 // -----------------------------------------------------------------------
9 //
9 //
10 // This file interfaces the MixMax PseudoRando 10 // This file interfaces the MixMax PseudoRandom Number Generator
11 // proposed by:
11 // proposed by:
12 //
12 //
13 // G.K.Savvidy and N.G.Ter-Arutyunian,
13 // G.K.Savvidy and N.G.Ter-Arutyunian,
14 // On the Monte Carlo simulation of physical 14 // On the Monte Carlo simulation of physical systems,
15 // J.Comput.Phys. 97, 566 (1991);
15 // J.Comput.Phys. 97, 566 (1991);
16 // Preprint EPI-865-16-86, Yerevan, Jan. 198 16 // Preprint EPI-865-16-86, Yerevan, Jan. 1986
17 // http://dx.doi.org/10.1016/0021-9991(91)90 17 // http://dx.doi.org/10.1016/0021-9991(91)90015-D
18 //
18 //
19 // K.Savvidy
19 // K.Savvidy
20 // "The MIXMAX random number generator"
20 // "The MIXMAX random number generator"
21 // Comp. Phys. Commun. (2015)
21 // Comp. Phys. Commun. (2015)
22 // http://dx.doi.org/10.1016/j.cpc.2015.06.0 22 // http://dx.doi.org/10.1016/j.cpc.2015.06.003
23 //
23 //
24 // K.Savvidy and G.Savvidy
24 // K.Savvidy and G.Savvidy
25 // "Spectrum and Entropy of C-systems. MIXMA 25 // "Spectrum and Entropy of C-systems. MIXMAX random number generator"
26 // Chaos, Solitons & Fractals, Volume 91, (2 26 // Chaos, Solitons & Fractals, Volume 91, (2016) pp. 33-38
27 // http://dx.doi.org/10.1016/j.chaos.2016.05 27 // http://dx.doi.org/10.1016/j.chaos.2016.05.003
28 //
28 //
29 // =========================================== 29 // =======================================================================
30 // Implementation by Konstantin Savvidy - Copy << 30 // Implementation by Konstantin Savvidy - Copyright 2004-2017
31 // July 2023 - Updated class structure upon su <<
32 // September 2023 - fix (re-)initialization fr <<
33 // =========================================== 31 // =======================================================================
34
32
35 #include "CLHEP/Random/Random.h"
33 #include "CLHEP/Random/Random.h"
36 #include "CLHEP/Random/MixMaxRng.h"
34 #include "CLHEP/Random/MixMaxRng.h"
37 #include "CLHEP/Random/engineIDulong.h"
35 #include "CLHEP/Random/engineIDulong.h"
38 #include "CLHEP/Utility/atomic_int.h"
36 #include "CLHEP/Utility/atomic_int.h"
39
37
40 #include <atomic>
<<
41 #include <cmath>
<<
42 #include <algorithm>
<<
43 #include <iostream>
<<
44 #include <string.h> // for strcmp
38 #include <string.h> // for strcmp
45 #include <vector>
<< 39 #include <cmath>
46
40
47 const unsigned long MASK32=0xffffffff;
41 const unsigned long MASK32=0xffffffff;
48
42
49 namespace CLHEP {
43 namespace CLHEP {
50
44
51 namespace {
45 namespace {
52 // Number of instances with automatic seed s 46 // Number of instances with automatic seed selection
53 CLHEP_ATOMIC_INT_TYPE numberOfEngines(0);
47 CLHEP_ATOMIC_INT_TYPE numberOfEngines(0);
54 }
48 }
55
49
56 static const int MarkerLen = 64; // Enough roo 50 static const int MarkerLen = 64; // Enough room to hold a begin or end marker.
57
51
58 MixMaxRng::MixMaxRng()
52 MixMaxRng::MixMaxRng()
59 : HepRandomEngine()
53 : HepRandomEngine()
60 {
54 {
61 int numEngines = ++numberOfEngines;
55 int numEngines = ++numberOfEngines;
62 setSeed(static_cast<long>(numEngines));
56 setSeed(static_cast<long>(numEngines));
63 }
57 }
64
58
65 MixMaxRng::MixMaxRng(long seed)
59 MixMaxRng::MixMaxRng(long seed)
66 : HepRandomEngine()
60 : HepRandomEngine()
67 {
61 {
68 theSeed=seed;
62 theSeed=seed;
69 setSeed(seed);
63 setSeed(seed);
70 }
64 }
71
65
72 MixMaxRng::MixMaxRng(std::istream& is)
66 MixMaxRng::MixMaxRng(std::istream& is)
73 : HepRandomEngine()
67 : HepRandomEngine()
74 {
68 {
75 get(is);
69 get(is);
76 }
70 }
77
71
78 MixMaxRng::~MixMaxRng()
72 MixMaxRng::~MixMaxRng()
79 {
73 {
80 }
74 }
81
75
82 MixMaxRng::MixMaxRng(const MixMaxRng& rng)
76 MixMaxRng::MixMaxRng(const MixMaxRng& rng)
83 : HepRandomEngine(rng)
77 : HepRandomEngine(rng)
84 {
78 {
85 std::copy(rng.V, rng.V+N, V);
<< 79 S.V = rng.S.V;
86 sumtot= rng.sumtot;
<< 80 S.sumtot= rng.S.sumtot;
87 counter= rng.counter;
<< 81 S.counter= rng.S.counter;
88 }
82 }
89
83
90 MixMaxRng& MixMaxRng::operator=(const MixMaxRn 84 MixMaxRng& MixMaxRng::operator=(const MixMaxRng& rng)
91 {
85 {
92 // Check assignment to self
86 // Check assignment to self
93 //
87 //
94 if (this == &rng) { return *this; }
88 if (this == &rng) { return *this; }
95
89
96 // Copy base class data
90 // Copy base class data
97 //
91 //
98 HepRandomEngine::operator=(rng);
92 HepRandomEngine::operator=(rng);
99
93
100 std::copy(rng.V, rng.V+N, V);
<< 94 S.V = rng.S.V;
101 sumtot= rng.sumtot;
<< 95 S.sumtot= rng.S.sumtot;
102 counter= rng.counter;
<< 96 S.counter= rng.S.counter;
103
97
104 return *this;
98 return *this;
105 }
99 }
106
100
107 void MixMaxRng::saveStatus( const char filenam 101 void MixMaxRng::saveStatus( const char filename[] ) const
108 {
102 {
109 // Create a C file-handle or an object that 103 // Create a C file-handle or an object that can accept the same output
110 FILE *fh= fopen(filename, "w");
104 FILE *fh= fopen(filename, "w");
111 if( fh )
105 if( fh )
112 {
106 {
113 int j;
107 int j;
114 fprintf(fh, "mixmax state, file version 1 108 fprintf(fh, "mixmax state, file version 1.0\n" );
115 fprintf(fh, "N=%u; V[N]={", rng_get_N() ) 109 fprintf(fh, "N=%u; V[N]={", rng_get_N() );
116 for (j=0; (j< (rng_get_N()-1) ); ++j) {
<< 110 for (j=0; (j< (rng_get_N()-1) ); j++) {
117 fprintf(fh, "%llu, ", (unsigned long << 111 fprintf(fh, "%llu, ", S.V[j] );
118 }
112 }
119 fprintf(fh, "%llu", (unsigned long long)V << 113 fprintf(fh, "%llu", S.V[rng_get_N()-1] );
120 fprintf(fh, "}; " );
114 fprintf(fh, "}; " );
121 fprintf(fh, "counter=%u; ", counter );
<< 115 fprintf(fh, "counter=%u; ", S.counter );
122 fprintf(fh, "sumtot=%llu;\n", (unsigned l << 116 fprintf(fh, "sumtot=%llu;\n", S.sumtot );
123 fclose(fh);
117 fclose(fh);
124 }
118 }
125 }
119 }
126
120
127 void MixMaxRng::restoreStatus( const char file 121 void MixMaxRng::restoreStatus( const char filename[] )
128 {
122 {
129 // a function for reading the state from a 123 // a function for reading the state from a file
130 FILE* fin;
124 FILE* fin;
131 if( ( fin = fopen(filename, "r") ) )
125 if( ( fin = fopen(filename, "r") ) )
132 {
126 {
133 char l=0;
127 char l=0;
134 while ( l != '{' ) { // 0x7B = "{"
128 while ( l != '{' ) { // 0x7B = "{"
135 l=fgetc(fin); // proceed until hitting 129 l=fgetc(fin); // proceed until hitting opening bracket
136 }
130 }
137 ungetc(' ', fin);
131 ungetc(' ', fin);
138 }
132 }
139 else
133 else
140 {
134 {
141 fprintf(stderr, "mixmax -> read_state: e 135 fprintf(stderr, "mixmax -> read_state: error reading file %s\n", filename);
142 throw std::runtime_error("Error in readi 136 throw std::runtime_error("Error in reading state file");
143 }
137 }
144
138
145 myuint_t vecVal;
139 myuint_t vecVal;
146 //printf("mixmax -> read_state: starting to 140 //printf("mixmax -> read_state: starting to read state from file\n");
147 if (!fscanf(fin, "%llu", (unsigned long lon << 141 if (!fscanf(fin, "%llu", &S.V[0]) )
148 {
142 {
149 fprintf(stderr, "mixmax -> read_state: er 143 fprintf(stderr, "mixmax -> read_state: error reading file %s\n", filename);
150 throw std::runtime_error("Error in readin 144 throw std::runtime_error("Error in reading state file");
151 }
145 }
152
146
153 for (int i = 1; i < rng_get_N(); ++i)
<< 147 int i;
>> 148 for( i = 1; i < rng_get_N(); i++)
154 {
149 {
155 if (!fscanf(fin, ", %llu", (unsigned long << 150 if (!fscanf(fin, ", %llu", &vecVal) )
156 {
151 {
157 fprintf(stderr, "mixmax -> read_state: 152 fprintf(stderr, "mixmax -> read_state: error reading vector component i=%d from file %s\n", i, filename);
158 throw std::runtime_error("Error in read 153 throw std::runtime_error("Error in reading state file");
159 }
154 }
160 if( vecVal <= MixMaxRng::M61 )
<< 155 if( vecVal <= MixMaxRng::M61 )
161 {
156 {
162 V[i] = vecVal;
<< 157 S.V[i] = vecVal;
163 }
158 }
164 else
159 else
165 {
160 {
166 fprintf(stderr, "mixmax -> read_state: 161 fprintf(stderr, "mixmax -> read_state: Invalid state vector value= %llu"
167 " ( must be less than %llu ) "
162 " ( must be less than %llu ) "
168 " obtained from reading file %s 163 " obtained from reading file %s\n"
169 , (unsigned long long)vecVal, ( << 164 , vecVal, MixMaxRng::M61, filename);
170 }
165 }
171 }
166 }
172
167
173 int incounter;
<< 168 int counter;
174 if (!fscanf( fin, "}; counter=%i; ", &incou << 169 if (!fscanf( fin, "}; counter=%i; ", &counter))
175 {
170 {
176 fprintf(stderr, "mixmax -> read_state: er 171 fprintf(stderr, "mixmax -> read_state: error reading counter from file %s\n", filename);
177 throw std::runtime_error("Error in readin 172 throw std::runtime_error("Error in reading state file");
178 }
173 }
179 if( incounter <= rng_get_N() )
<< 174 if( counter <= rng_get_N() )
180 {
175 {
181 counter = incounter;
<< 176 S.counter= counter;
182 }
177 }
183 else
178 else
184 {
179 {
185 fprintf(stderr, "mixmax -> read_state: In 180 fprintf(stderr, "mixmax -> read_state: Invalid counter = %d"
186 " Must be 0 <= counter < %u\n" , 181 " Must be 0 <= counter < %u\n" , counter, rng_get_N());
187 print_state();
182 print_state();
188 throw std::runtime_error("Error in readin 183 throw std::runtime_error("Error in reading state counter");
189 }
184 }
190 precalc();
185 precalc();
191 myuint_t insumtot;
<< 186 myuint_t sumtot;
192 if (!fscanf( fin, "sumtot=%llu\n", (unsigne << 187 if (!fscanf( fin, "sumtot=%llu\n", &sumtot))
193 {
188 {
194 fprintf(stderr, "mixmax -> read_state: er 189 fprintf(stderr, "mixmax -> read_state: error reading checksum from file %s\n", filename);
195 throw std::runtime_error("Error in readin 190 throw std::runtime_error("Error in reading state file");
196 }
191 }
197
192
198 if (sumtot != insumtot)
<< 193 if (S.sumtot != sumtot)
199 {
194 {
200 fprintf(stderr, "mixmax -> checksum error 195 fprintf(stderr, "mixmax -> checksum error while reading state from file %s - corrupted?\n", filename);
201 throw std::runtime_error("Error in readin 196 throw std::runtime_error("Error in reading state checksum");
202 }
197 }
203 fclose(fin);
198 fclose(fin);
204 }
199 }
205
200
206 void MixMaxRng::showStatus() const
201 void MixMaxRng::showStatus() const
207 {
202 {
208 std::cout << std::endl;
203 std::cout << std::endl;
209 std::cout << "------- MixMaxRng engine stat 204 std::cout << "------- MixMaxRng engine status -------" << std::endl;
210
205
211 std::cout << " Current state vector is:" << 206 std::cout << " Current state vector is:" << std::endl;
212 print_state();
207 print_state();
213 std::cout << "----------------------------- 208 std::cout << "---------------------------------------" << std::endl;
214 }
209 }
215
210
>> 211 void MixMaxRng::setSeed(long longSeed, int /* extraSeed */)
>> 212 {
>> 213 //seed_uniquestream(0,0,0,longSeed);
>> 214 theSeed = longSeed;
>> 215 seed_spbox(longSeed);
>> 216 }
>> 217
216 // Preferred Seeding method
218 // Preferred Seeding method
217 // the values of 'Seeds' must be valid 32-bit 219 // the values of 'Seeds' must be valid 32-bit integers
218 // Higher order bits will be ignored!!
220 // Higher order bits will be ignored!!
219
221
220 void MixMaxRng::setSeeds(const long* Seeds, in 222 void MixMaxRng::setSeeds(const long* Seeds, int seedNum)
221 {
223 {
222 myID_t seed0, seed1= 0, seed2= 0, seed3= 0; << 224 unsigned long seed0, seed1= 0, seed2= 0, seed3= 0;
223
225
224 if( seedNum < 1 ) { // Assuming at least 2 226 if( seedNum < 1 ) { // Assuming at least 2 seeds in vector...
225 seed0= static_cast<myID_t>(Seeds[0]) & << 227 seed0= static_cast<unsigned long>(Seeds[0]) & MASK32;
226 seed1= static_cast<myID_t>(Seeds[1]) & << 228 seed1= static_cast<unsigned long>(Seeds[1]) & MASK32;
227 }
229 }
228 else
230 else
229 {
231 {
230 if( seedNum < 4 ) {
232 if( seedNum < 4 ) {
231 seed0= static_cast<myID_t>(Seeds[0]) & << 233 seed0= static_cast<unsigned long>(Seeds[0]) & MASK32;
232 if( seedNum > 1){ seed1= static_cast<my << 234 if( seedNum > 1){ seed1= static_cast<unsigned long>(Seeds[1]) & MASK32; }
233 if( seedNum > 2){ seed2= static_cast<my << 235 if( seedNum > 2){ seed2= static_cast<unsigned long>(Seeds[2]) & MASK32; }
234 }
236 }
235 if( seedNum >= 4 ) {
237 if( seedNum >= 4 ) {
236 seed0= static_cast<myID_t>(Seeds[0]) & << 238 seed0= static_cast<unsigned long>(Seeds[0]) & MASK32;
237 seed1= static_cast<myID_t>(Seeds[1]) & << 239 seed1= static_cast<unsigned long>(Seeds[1]) & MASK32;
238 seed2= static_cast<myID_t>(Seeds[2]) & << 240 seed2= static_cast<unsigned long>(Seeds[2]) & MASK32;
239 seed3= static_cast<myID_t>(Seeds[3]) & << 241 seed3= static_cast<unsigned long>(Seeds[3]) & MASK32;
240 }
242 }
241 }
243 }
242 theSeed = Seeds[0];
244 theSeed = Seeds[0];
243 theSeeds = Seeds;
245 theSeeds = Seeds;
244 seed_uniquestream(seed3, seed2, seed1, seed 246 seed_uniquestream(seed3, seed2, seed1, seed0);
245 }
247 }
246
248
247 std::string MixMaxRng::engineName()
249 std::string MixMaxRng::engineName()
248 {
250 {
249 return "MixMaxRng";
251 return "MixMaxRng";
250 }
252 }
251
253
252 constexpr int MixMaxRng::rng_get_N()
254 constexpr int MixMaxRng::rng_get_N()
253 {
255 {
254 return N;
256 return N;
255 }
257 }
256
258
>> 259 constexpr long long int MixMaxRng::rng_get_SPECIAL()
>> 260 {
>> 261 return SPECIAL;
>> 262 }
>> 263
>> 264 constexpr int MixMaxRng::rng_get_SPECIALMUL()
>> 265 {
>> 266 return SPECIALMUL;
>> 267 }
>> 268
>> 269 double MixMaxRng::generate(int i)
>> 270 {
>> 271 S.counter++;
>> 272 #if defined(__clang__) || defined(__llvm__)
>> 273 return INV_M61*static_cast<double>(S.V[i]);
>> 274 #elif defined(__GNUC__) && (__GNUC__ < 7) && (!defined(__ICC)) && defined(__x86_64__) && defined(__SSE2_MATH__)
>> 275 int64_t Z=S.V[i];
>> 276 double F=0.0;
>> 277 //#warning Using the inline assembler
>> 278 /* using SSE inline assemly to zero the xmm register, just before int64 -> double conversion,
>> 279 not necessary in GCC-5 or better, but huge penalty on earlier compilers
>> 280 */
>> 281 __asm__ __volatile__( "pxor %0, %0;"
>> 282 "cvtsi2sdq %1, %0;"
>> 283 :"=x"(F)
>> 284 :"r"(Z)
>> 285 );
>> 286 return F*INV_M61;
>> 287 #else
>> 288 //#warning other method
>> 289 return convert1double(S.V[i]); //get_next_float_packbits();
>> 290 #endif
>> 291 }
>> 292
>> 293 double MixMaxRng::iterate()
>> 294 {
>> 295 myuint_t* Y=S.V.data();
>> 296 myuint_t tempP, tempV;
>> 297 Y[0] = ( tempV = S.sumtot);
>> 298 myuint_t sumtot = Y[0], ovflow = 0; // will keep a running sum of all new elements
>> 299 tempP = 0; // will keep a partial sum of all old elements
>> 300 myuint_t tempPO;
>> 301 tempPO = MULWU(tempP); tempP = modadd(tempP, Y[1] ); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); Y[1] = tempV; sumtot += tempV; if (sumtot < tempV) {ovflow++;};
>> 302 tempPO = MULWU(tempP); tempP = modadd(tempP, Y[2] ); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); Y[2] = tempV; sumtot += tempV; if (sumtot < tempV) {ovflow++;};
>> 303 tempPO = MULWU(tempP); tempP = modadd(tempP, Y[3] ); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); Y[3] = tempV; sumtot += tempV; if (sumtot < tempV) {ovflow++;};
>> 304 tempPO = MULWU(tempP); tempP = modadd(tempP, Y[4] ); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); Y[4] = tempV; sumtot += tempV; if (sumtot < tempV) {ovflow++;};
>> 305 tempPO = MULWU(tempP); tempP = modadd(tempP, Y[5] ); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); Y[5] = tempV; sumtot += tempV; if (sumtot < tempV) {ovflow++;};
>> 306 tempPO = MULWU(tempP); tempP = modadd(tempP, Y[6] ); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); Y[6] = tempV; sumtot += tempV; if (sumtot < tempV) {ovflow++;};
>> 307 tempPO = MULWU(tempP); tempP = modadd(tempP, Y[7] ); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); Y[7] = tempV; sumtot += tempV; if (sumtot < tempV) {ovflow++;};
>> 308 tempPO = MULWU(tempP); tempP = modadd(tempP, Y[8] ); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); Y[8] = tempV; sumtot += tempV; if (sumtot < tempV) {ovflow++;};
>> 309 tempPO = MULWU(tempP); tempP = modadd(tempP, Y[9] ); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); Y[9] = tempV; sumtot += tempV; if (sumtot < tempV) {ovflow++;};
>> 310 tempPO = MULWU(tempP); tempP = modadd(tempP, Y[10]); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); Y[10] = tempV; sumtot += tempV; if (sumtot < tempV) {ovflow++;};
>> 311 tempPO = MULWU(tempP); tempP = modadd(tempP, Y[11]); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); Y[11] = tempV; sumtot += tempV; if (sumtot < tempV) {ovflow++;};
>> 312 tempPO = MULWU(tempP); tempP = modadd(tempP, Y[12]); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); Y[12] = tempV; sumtot += tempV; if (sumtot < tempV) {ovflow++;};
>> 313 tempPO = MULWU(tempP); tempP = modadd(tempP, Y[13]); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); Y[13] = tempV; sumtot += tempV; if (sumtot < tempV) {ovflow++;};
>> 314 tempPO = MULWU(tempP); tempP = modadd(tempP, Y[14]); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); Y[14] = tempV; sumtot += tempV; if (sumtot < tempV) {ovflow++;};
>> 315 tempPO = MULWU(tempP); tempP = modadd(tempP, Y[15]); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); Y[15] = tempV; sumtot += tempV; if (sumtot < tempV) {ovflow++;};
>> 316 tempPO = MULWU(tempP); tempP = modadd(tempP, Y[16]); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); Y[16] = tempV; sumtot += tempV; if (sumtot < tempV) {ovflow++;};
>> 317 S.sumtot = MIXMAX_MOD_MERSENNE(MIXMAX_MOD_MERSENNE(sumtot) + (ovflow <<3 ));
>> 318
>> 319 S.counter=2;
>> 320 return double(S.V[1])*INV_M61;
>> 321 }
>> 322
257 void MixMaxRng::flatArray(const int size, doub 323 void MixMaxRng::flatArray(const int size, double* vect )
258 {
324 {
>> 325 // fill_array( S, size, arrayDbl );
259 for (int i=0; i<size; ++i) { vect[i] = flat 326 for (int i=0; i<size; ++i) { vect[i] = flat(); }
260 }
327 }
261
328
>> 329 MixMaxRng::operator double()
>> 330 {
>> 331 return flat();
>> 332 }
>> 333
>> 334 MixMaxRng::operator float()
>> 335 {
>> 336 return float( flat() );
>> 337 }
>> 338
>> 339 MixMaxRng::operator unsigned int()
>> 340 {
>> 341 return static_cast<unsigned int>(get_next());
>> 342 // clhep_get_next returns a 64-bit integer, of which the lower 61 bits
>> 343 // are random and upper 3 bits are zero
>> 344 }
>> 345
262 std::ostream & MixMaxRng::put ( std::ostream& 346 std::ostream & MixMaxRng::put ( std::ostream& os ) const
263 {
347 {
264 char beginMarker[] = "MixMaxRng-begin";
348 char beginMarker[] = "MixMaxRng-begin";
265 char endMarker[] = "MixMaxRng-end";
349 char endMarker[] = "MixMaxRng-end";
266
350
267 long pr = os.precision(24);
<< 351 int pr = os.precision(24);
268 os << beginMarker << " ";
352 os << beginMarker << " ";
269 os << theSeed << "\n";
353 os << theSeed << "\n";
270 for (int i=0; i<rng_get_N(); ++i) {
354 for (int i=0; i<rng_get_N(); ++i) {
271 os << V[i] << "\n";
<< 355 os << S.V[i] << "\n";
272 }
356 }
273 os << counter << "\n";
<< 357 os << S.counter << "\n";
274 os << sumtot << "\n";
<< 358 os << S.sumtot << "\n";
275 os << endMarker << "\n";
359 os << endMarker << "\n";
276 os.precision(pr);
360 os.precision(pr);
277 return os;
361 return os;
278 }
362 }
279
363
280 std::vector<unsigned long> MixMaxRng::put () c 364 std::vector<unsigned long> MixMaxRng::put () const
281 {
365 {
282 std::vector<unsigned long> vec;
<< 366 std::vector<unsigned long> v;
283 vec.push_back (engineIDulong<MixMaxRng>()); << 367 v.push_back (engineIDulong<MixMaxRng>());
284 for (int i=0; i<rng_get_N(); ++i)
368 for (int i=0; i<rng_get_N(); ++i)
285 {
369 {
286 vec.push_back(static_cast<unsigned long>( << 370 v.push_back(static_cast<unsigned long>(S.V[i] & MASK32));
287 // little-ended order on all platforms
371 // little-ended order on all platforms
288 vec.push_back(static_cast<unsigned long>( << 372 v.push_back(static_cast<unsigned long>(S.V[i] >> 32 ));
289 // pack uint64 into a data structure wh 373 // pack uint64 into a data structure which is 32-bit on some platforms
290 }
374 }
291 vec.push_back(static_cast<unsigned long>(co << 375 v.push_back(static_cast<unsigned long>(S.counter));
292 vec.push_back(static_cast<unsigned long>(su << 376 v.push_back(static_cast<unsigned long>(S.sumtot & MASK32));
293 vec.push_back(static_cast<unsigned long>(su << 377 v.push_back(static_cast<unsigned long>(S.sumtot >> 32));
294 return vec;
<< 378 return v;
295 }
379 }
296
380
297 std::istream & MixMaxRng::get ( std::istream& 381 std::istream & MixMaxRng::get ( std::istream& is)
298 {
382 {
299 char beginMarker [MarkerLen];
383 char beginMarker [MarkerLen];
300 is >> std::ws;
384 is >> std::ws;
301 is.width(MarkerLen); // causes the next re 385 is.width(MarkerLen); // causes the next read to the char* to be <=
302 // that many bytes, I 386 // that many bytes, INCLUDING A TERMINATION \0
303 // (Stroustrup, secti 387 // (Stroustrup, section 21.3.2)
304 is >> beginMarker;
388 is >> beginMarker;
305 if (strcmp(beginMarker,"MixMaxRng-begin")) 389 if (strcmp(beginMarker,"MixMaxRng-begin")) {
306 is.clear(std::ios::badbit | is.rdstate() 390 is.clear(std::ios::badbit | is.rdstate());
307 std::cerr << "\nInput stream misposition 391 std::cerr << "\nInput stream mispositioned or"
308 << "\nMixMaxRng state descript 392 << "\nMixMaxRng state description missing or"
309 << "\nwrong engine type found. 393 << "\nwrong engine type found." << std::endl;
310 return is;
394 return is;
311 }
395 }
312 return getState(is);
396 return getState(is);
313 }
397 }
314
398
315 std::string MixMaxRng::beginTag ()
399 std::string MixMaxRng::beginTag ()
316 {
400 {
317 return "MixMaxRng-begin";
401 return "MixMaxRng-begin";
318 }
402 }
319
403
320 std::istream & MixMaxRng::getState ( std::ist 404 std::istream & MixMaxRng::getState ( std::istream& is )
321 {
405 {
322 char endMarker[MarkerLen];
406 char endMarker[MarkerLen];
323 is >> theSeed;
407 is >> theSeed;
324 for (int i=0; i<rng_get_N(); ++i) is >> V[ << 408 for (int i=0; i<rng_get_N(); ++i) is >> S.V[i];
325 is >> counter;
<< 409 is >> S.counter;
326 myuint_t checksum;
410 myuint_t checksum;
327 is >> checksum;
411 is >> checksum;
328 is >> std::ws;
412 is >> std::ws;
329 is.width(MarkerLen);
413 is.width(MarkerLen);
330 is >> endMarker;
414 is >> endMarker;
331 if (strcmp(endMarker,"MixMaxRng-end")) {
415 if (strcmp(endMarker,"MixMaxRng-end")) {
332 is.clear(std::ios::badbit | is.rdstate( 416 is.clear(std::ios::badbit | is.rdstate());
333 std::cerr << "\nMixMaxRng state descrip 417 std::cerr << "\nMixMaxRng state description incomplete."
334 << "\nInput stream is probabl 418 << "\nInput stream is probably mispositioned now.\n";
335 return is;
419 return is;
336 }
420 }
337 if ( counter < 0 || counter > rng_get_N() ) << 421 if ( S.counter < 0 || S.counter > rng_get_N() ) {
338 std::cerr << "\nMixMaxRng::getState(): 422 std::cerr << "\nMixMaxRng::getState(): "
339 << "vector read wrong value o 423 << "vector read wrong value of counter from file!"
340 << "\nInput stream is probabl 424 << "\nInput stream is probably mispositioned now.\n";
341 return is;
425 return is;
342 }
426 }
343 precalc();
427 precalc();
344 if ( checksum != sumtot) {
<< 428 if ( checksum != S.sumtot) {
345 std::cerr << "\nMixMaxRng::getState(): 429 std::cerr << "\nMixMaxRng::getState(): "
346 << "checksum disagrees with v 430 << "checksum disagrees with value stored in file!"
347 << "\nInput stream is probabl 431 << "\nInput stream is probably mispositioned now.\n";
348 return is;
432 return is;
349 }
433 }
350 return is;
434 return is;
351 }
435 }
352
436
353 bool MixMaxRng::get (const std::vector<unsigne << 437 bool MixMaxRng::get (const std::vector<unsigned long> & v)
354 {
438 {
355 if ((vec[0] & 0xffffffffUL) != engineIDulon << 439 if ((v[0] & 0xffffffffUL) != engineIDulong<MixMaxRng>()) {
356 {
<<
357 std::cerr <<
440 std::cerr <<
358 "\nMixMaxRng::get(): vector has wrong 441 "\nMixMaxRng::get(): vector has wrong ID word - state unchanged\n";
359 return false;
442 return false;
360 }
443 }
361 return getState(vec);
<< 444 return getState(v);
362 }
445 }
363
446
364 bool MixMaxRng::getState (const std::vector<un << 447 bool MixMaxRng::getState (const std::vector<unsigned long> & v)
365 {
448 {
366 if (vec.size() != VECTOR_STATE_SIZE ) {
<< 449 if (v.size() != VECTOR_STATE_SIZE ) {
367 std::cerr <<
450 std::cerr <<
368 "\nMixMaxRng::getState(): vector has w 451 "\nMixMaxRng::getState(): vector has wrong length - state unchanged\n";
369 return false;
452 return false;
370 }
453 }
371 for (int i=1; i<2*rng_get_N() ; i=i+2) {
454 for (int i=1; i<2*rng_get_N() ; i=i+2) {
372 V[i/2]= ( (vec[i] & MASK32) | ( (myuint_t << 455 S.V[i/2]= ( (v[i] & MASK32) | ( (myuint_t)(v[i+1]) << 32 ) );
373 // unpack from a data structure which is 456 // unpack from a data structure which is 32-bit on some platforms
374 }
457 }
375 counter = (int)vec[2*rng_get_N()+1];
<< 458 S.counter = v[2*rng_get_N()+1];
376 precalc();
459 precalc();
377 if ( ( (vec[2*rng_get_N()+2] & MASK32)
<< 460 if ( ( (v[2*rng_get_N()+2] & MASK32)
378 | ( (myuint_t)(vec[2*rng_get_N()+3]) < << 461 | ( (myuint_t)(v[2*rng_get_N()+3]) << 32 ) ) != S.sumtot) {
379 std::cerr << "\nMixMaxRng::getState(): ve 462 std::cerr << "\nMixMaxRng::getState(): vector has wrong checksum!"
380 << "\nInput vector is probably 463 << "\nInput vector is probably mispositioned now.\n";
381 return false;
464 return false;
382 }
465 }
383 return true;
466 return true;
384 }
467 }
385
468
>> 469 myuint_t MixMaxRng ::MOD_MULSPEC(myuint_t k)
>> 470 {
>> 471 switch (N)
>> 472 {
>> 473 case 17:
>> 474 return 0;
>> 475 break;
>> 476 case 8:
>> 477 return 0;
>> 478 break;
>> 479 case 240:
>> 480 return fmodmulM61( 0, SPECIAL , (k) );
>> 481 break;
>> 482 default:
>> 483 std::cerr << "MIXMAX ERROR: " << "Disallowed value of parameter N\n";
>> 484 std::terminate();
>> 485 break;
>> 486 }
>> 487 }
>> 488
>> 489 myuint_t MixMaxRng::MULWU (myuint_t k)
>> 490 {
>> 491 return (( (k)<<(SPECIALMUL) & M61) ^ ( (k) >> (BITS-SPECIALMUL)) );
>> 492 }
>> 493
386 myuint_t MixMaxRng::iterate_raw_vec(myuint_t* 494 myuint_t MixMaxRng::iterate_raw_vec(myuint_t* Y, myuint_t sumtotOld)
387 {
495 {
388 // operates with a raw vector, uses known s 496 // operates with a raw vector, uses known sum of elements of Y
>> 497 int i;
389
498
390 myuint_t tempP, tempV;
499 myuint_t tempP, tempV;
391 Y[0] = ( tempV = sumtotOld);
500 Y[0] = ( tempV = sumtotOld);
392 myuint_t insumtot = Y[0], ovflow = 0; // wi << 501 myuint_t sumtot = Y[0], ovflow = 0; // will keep a running sum of all new elements
393 tempP = 0; // will keep a part 502 tempP = 0; // will keep a partial sum of all old elements
394 for (int i=1; (i<N); ++i)
<< 503 for (i=1; (i<N); i++)
395 {
504 {
396 myuint_t tempPO = MULWU(tempP);
505 myuint_t tempPO = MULWU(tempP);
397 tempP = modadd(tempP, Y[i]);
506 tempP = modadd(tempP, Y[i]);
398 tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+t 507 tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); // new Y[i] = old Y[i] + old partial * m
399 Y[i] = tempV;
508 Y[i] = tempV;
400 insumtot += tempV; if (insumtot < tempV) << 509 sumtot += tempV; if (sumtot < tempV) {ovflow++;}
401 }
510 }
402 return MIXMAX_MOD_MERSENNE(MIXMAX_MOD_MERSE << 511 return MIXMAX_MOD_MERSENNE(MIXMAX_MOD_MERSENNE(sumtot) + (ovflow <<3 ));
403 }
512 }
404
513
405 myuint_t MixMaxRng::get_next()
514 myuint_t MixMaxRng::get_next()
406 {
515 {
407 int i = counter;
<< 516 int i;
>> 517 i=S.counter;
408
518
409 if ((i<=(N-1)) )
519 if ((i<=(N-1)) )
410 {
520 {
411 ++counter;
<< 521 S.counter++;
412 return V[i];
<< 522 return S.V[i];
413 }
523 }
414 else
524 else
415 {
525 {
416 sumtot = iterate_raw_vec(V, sumtot);
<< 526 S.sumtot = iterate_raw_vec(S.V.data(), S.sumtot);
417 counter=2;
<< 527 S.counter=2;
418 return V[1];
<< 528 return S.V[1];
419 }
529 }
420 }
530 }
421
531
422 myuint_t MixMaxRng::precalc()
532 myuint_t MixMaxRng::precalc()
423 {
533 {
424 myuint_t temp = 0;
<< 534 int i;
425 for (int i=0; i < N; ++i) { temp = MIXMAX_M << 535 myuint_t temp;
426 sumtot = temp;
<< 536 temp = 0;
>> 537 for (i=0; i < N; i++){
>> 538 temp = MIXMAX_MOD_MERSENNE(temp + S.V[i]);
>> 539 }
>> 540 S.sumtot = temp;
427 return temp;
541 return temp;
428 }
542 }
>> 543
>> 544 double MixMaxRng::get_next_float_packbits()
>> 545 {
>> 546 myuint_t Z=get_next();
>> 547 return convert1double(Z);
>> 548 }
429
549
430 void MixMaxRng::state_init()
<< 550 void MixMaxRng::seed_vielbein(unsigned int index)
431 {
551 {
432 for (int i=1; i < N; ++i) { V[i] = 0; }
<< 552 int i;
433 V[0] = 1;
<< 553 if (index<N)
434 counter = N; // set the counter to N if it << 554 {
435 sumtot = 1;
<< 555 for (i=0; i < N; i++){
>> 556 S.V[i] = 0;
>> 557 }
>> 558 S.V[index] = 1;
>> 559 }
>> 560 else
>> 561 {
>> 562 std::terminate();
>> 563 }
>> 564 S.counter = N; // set the counter to N if iteration should happen right away
>> 565 S.sumtot = 1;
436 }
566 }
437
567
438 void MixMaxRng::seed_spbox(myuint_t seed)
568 void MixMaxRng::seed_spbox(myuint_t seed)
439 {
569 {
440 // a 64-bit LCG from Knuth line 26, in comb 570 // a 64-bit LCG from Knuth line 26, in combination with a bit swap is used to seed
441
571
442 if (seed == 0)
<<
443 throw std::runtime_error("try seeding wit <<
444
<<
445 const myuint_t MULT64=6364136223846793005UL 572 const myuint_t MULT64=6364136223846793005ULL;
446 sumtot = 0;
<< 573 int i;
>> 574
>> 575 myuint_t sumtot=0,ovflow=0;
>> 576 if (seed == 0) throw std::runtime_error("try seeding with nonzero seed next time");
447
577
448 myuint_t l = seed;
578 myuint_t l = seed;
449
579
450 for (int i=0; i < N; ++i)
<< 580 for (i=0; i < N; i++){
451 {
<<
452 l*=MULT64; l = (l << 32) ^ (l>>32);
581 l*=MULT64; l = (l << 32) ^ (l>>32);
453 V[i] = l & M61;
<< 582 S.V[i] = l & M61;
454 sumtot = MIXMAX_MOD_MERSENNE(sumtot + V[( << 583 sumtot += S.V[(i)]; if (sumtot < S.V[(i)]) {ovflow++;}
455 }
584 }
456 counter = N; // set the counter to N if it << 585 S.counter = N; // set the counter to N if iteration should happen right away
>> 586 S.sumtot = MIXMAX_MOD_MERSENNE(MIXMAX_MOD_MERSENNE(sumtot) + (ovflow <<3 ));
457 }
587 }
458
588
459 void MixMaxRng::seed_uniquestream( myID_t clus 589 void MixMaxRng::seed_uniquestream( myID_t clusterID, myID_t machineID, myID_t runID, myID_t streamID )
460 {
590 {
461 state_init();
<< 591 seed_vielbein(0);
462 sumtot = apply_bigskip(V, V, clusterID, mac << 592 S.sumtot = apply_bigskip(S.V.data(), S.V.data(), clusterID, machineID, runID, streamID );
463 counter = 1;
<< 593 S.counter = 1;
464 }
594 }
465
595
466 myuint_t MixMaxRng::apply_bigskip( myuint_t* V 596 myuint_t MixMaxRng::apply_bigskip( myuint_t* Vout, myuint_t* Vin, myID_t clusterID, myID_t machineID, myID_t runID, myID_t streamID )
467 {
597 {
468 /*
598 /*
469 makes a derived state vector, Vout, from t 599 makes a derived state vector, Vout, from the mother state vector Vin
470 by skipping a large number of steps, deter 600 by skipping a large number of steps, determined by the given seeding ID's
471
601
472 it is mathematically guaranteed that the s 602 it is mathematically guaranteed that the substreams derived in this way from the SAME (!!!) Vin will not collide provided
473 1) at least one bit of ID is different
603 1) at least one bit of ID is different
474 2) less than 10^100 numbers are drawn from 604 2) less than 10^100 numbers are drawn from the stream
475 (this is good enough : a single CPU will n 605 (this is good enough : a single CPU will not exceed this in the lifetime of the universe, 10^19 sec,
476 even if it had a clock cycle of Planch tim 606 even if it had a clock cycle of Planch time, 10^44 Hz )
477
607
478 Caution: never apply this to a derived vec << 608 Caution: never apply this to a derived vector, just choose some mother vector Vin, for example the unit vector by seed_vielbein(X,0),
479 and use it in all your runs, just change r 609 and use it in all your runs, just change runID to get completely nonoverlapping streams of random numbers on a different day.
480
610
481 clusterID and machineID are provided for t 611 clusterID and machineID are provided for the benefit of large organizations who wish to ensure that a simulation
482 which is running in parallel on a large nu 612 which is running in parallel on a large number of clusters and machines will have non-colliding source of random numbers.
483
613
484 did i repeat it enough times? the non-coll 614 did i repeat it enough times? the non-collision guarantee is absolute, not probabilistic
485
615
486 */
616 */
487
617
488 const myuint_t skipMat17[128][17] =
618 const myuint_t skipMat17[128][17] =
489 #include "CLHEP/Random/mixmax_skip_N17.icc" 619 #include "CLHEP/Random/mixmax_skip_N17.icc"
490 ;
620 ;
491
621
492 const myuint_t* skipMat[128];
622 const myuint_t* skipMat[128];
493 for (int i=0; i<128; ++i) { skipMat[i] = sk << 623 for (int i=0; i<128; i++) { skipMat[i] = skipMat17[i];}
494
624
495 myID_t IDvec[4] = {streamID, runID, machine 625 myID_t IDvec[4] = {streamID, runID, machineID, clusterID};
496 int r,i,j, IDindex;
626 int r,i,j, IDindex;
497 myID_t id;
627 myID_t id;
498 myuint_t Y[N], cum[N];
628 myuint_t Y[N], cum[N];
499 myuint_t coeff;
629 myuint_t coeff;
500 myuint_t* rowPtr;
630 myuint_t* rowPtr;
501 myuint_t insumtot=0;
<< 631 myuint_t sumtot=0;
502
632
503 for (i=0; i<N; ++i) { Y[i] = Vin[i]; insumt << 633 for (i=0; i<N; i++) { Y[i] = Vin[i]; sumtot = modadd( sumtot, Vin[i]); } ;
504 for (IDindex=0; IDindex<4; ++IDindex)
<< 634 for (IDindex=0; IDindex<4; IDindex++)
505 { // go from lower order to higher order ID 635 { // go from lower order to higher order ID
506 id=IDvec[IDindex];
636 id=IDvec[IDindex];
507 //printf("now doing ID at level %d, with 637 //printf("now doing ID at level %d, with ID = %d\n", IDindex, id);
508 r = 0;
638 r = 0;
509 while (id)
639 while (id)
510 {
640 {
511 if (id & 1)
641 if (id & 1)
512 {
642 {
513 rowPtr = (myuint_t*)skipMat[r + IDind 643 rowPtr = (myuint_t*)skipMat[r + IDindex*8*sizeof(myID_t)];
514 for (i=0; i<N; ++i){ cum[i] = 0; }
<< 644 for (i=0; i<N; i++){ cum[i] = 0; }
515 for (j=0; j<N; ++j)
<< 645 for (j=0; j<N; j++)
516 { // j is lag, enumerates terms of th 646 { // j is lag, enumerates terms of the poly
517 // for zero lag Y is already given
647 // for zero lag Y is already given
518 coeff = rowPtr[j]; // same coeff fo 648 coeff = rowPtr[j]; // same coeff for all i
519 for (i =0; i<N; ++i){
<< 649 for (i =0; i<N; i++){
520 cum[i] = fmodmulM61( cum[i], coe 650 cum[i] = fmodmulM61( cum[i], coeff , Y[i] ) ;
521 }
651 }
522 insumtot = iterate_raw_vec(Y, insum << 652 sumtot = iterate_raw_vec(Y, sumtot);
523 }
653 }
524 insumtot=0;
<< 654 sumtot=0;
525 for (i=0; i<N; ++i){ Y[i] = cum[i]; i << 655 for (i=0; i<N; i++){ Y[i] = cum[i]; sumtot = modadd( sumtot, cum[i]); } ;
526 }
656 }
527 id = (id >> 1); ++r; // bring up the r- << 657 id = (id >> 1); r++; // bring up the r-th bit in the ID
528 }
658 }
529 }
659 }
530 insumtot=0;
<< 660 sumtot=0;
531 for (i=0; i<N; ++i){ Vout[i] = Y[i]; insumt << 661 for (i=0; i<N; i++){ Vout[i] = Y[i]; sumtot = modadd( sumtot, Y[i]); }
532 // returns sumtot, and copy the vector over 662 // returns sumtot, and copy the vector over to Vout
533 return insumtot;
<< 663 return (sumtot) ;
534 }
664 }
535
665
536 #if defined(__x86_64__)
666 #if defined(__x86_64__)
537 myuint_t MixMaxRng::mod128(__uint128_t s)
667 myuint_t MixMaxRng::mod128(__uint128_t s)
538 {
668 {
539 myuint_t s1;
669 myuint_t s1;
540 s1 = ( ( ((myuint_t)s)&M61 ) + ( ((myuin 670 s1 = ( ( ((myuint_t)s)&M61 ) + ( ((myuint_t)(s>>64)) * 8 ) + ( ((myuint_t)s) >>BITS) );
541 return MIXMAX_MOD_MERSENNE(s1);
671 return MIXMAX_MOD_MERSENNE(s1);
542 }
672 }
543 myuint_t MixMaxRng::fmodmulM61(myuint_t cum, 673 myuint_t MixMaxRng::fmodmulM61(myuint_t cum, myuint_t a, myuint_t b)
544 {
674 {
545 __uint128_t temp;
675 __uint128_t temp;
546 temp = (__uint128_t)a*(__uint128_t)b + cum 676 temp = (__uint128_t)a*(__uint128_t)b + cum;
547 return mod128(temp);
677 return mod128(temp);
548 }
678 }
549 #else // on all other platforms, including 32- 679 #else // on all other platforms, including 32-bit linux, PPC and PPC64, ARM and all Windows
550 myuint_t MixMaxRng::fmodmulM61(myuint_t cum, 680 myuint_t MixMaxRng::fmodmulM61(myuint_t cum, myuint_t s, myuint_t a)
551 {
681 {
552 const myuint_t MASK32=0xFFFFFFFFULL;
682 const myuint_t MASK32=0xFFFFFFFFULL;
553 myuint_t o,ph,pl,ah,al;
683 myuint_t o,ph,pl,ah,al;
554 o=(s)*a;
684 o=(s)*a;
555 ph = ((s)>>32);
685 ph = ((s)>>32);
556 pl = (s) & MASK32;
686 pl = (s) & MASK32;
557 ah = a>>32;
687 ah = a>>32;
558 al = a & MASK32;
688 al = a & MASK32;
559 o = (o & M61) + ((ph*ah)<<3) + ((ah*pl+al* 689 o = (o & M61) + ((ph*ah)<<3) + ((ah*pl+al*ph + ((al*pl)>>32))>>29) ;
560 o += cum;
690 o += cum;
561 o = (o & M61) + ((o>>61));
691 o = (o & M61) + ((o>>61));
562 return o;
692 return o;
563 }
693 }
564 #endif
694 #endif
565
695
>> 696 myuint_t MixMaxRng::modadd(myuint_t foo, myuint_t bar)
>> 697 {
>> 698 #if defined(__x86_64__) && defined(__GNUC__) && (!defined(__ICC))
>> 699 //#warning Using assembler routine in modadd
>> 700 myuint_t out;
>> 701 /* Assembler trick suggested by Andrzej Gòˆrlich */
>> 702 __asm__ ("addq %2, %0; "
>> 703 "btrq $61, %0; "
>> 704 "adcq $0, %0; "
>> 705 :"=r"(out)
>> 706 :"0"(foo), "r"(bar)
>> 707 );
>> 708 return out;
>> 709 #else
>> 710 return MIXMAX_MOD_MERSENNE(foo+bar);
>> 711 #endif
>> 712 }
>> 713
566 void MixMaxRng::print_state() const
714 void MixMaxRng::print_state() const
567 {
715 {
>> 716 int j;
568 std::cout << "mixmax state, file version 1. 717 std::cout << "mixmax state, file version 1.0\n";
569 std::cout << "N=" << rng_get_N() << "; V[N] 718 std::cout << "N=" << rng_get_N() << "; V[N]={";
570 for (int j=0; (j< (rng_get_N()-1) ); ++j) { << 719 for (j=0; (j< (rng_get_N()-1) ); j++) {
571 std::cout << V[rng_get_N()-1];
<< 720 std::cout << S.V[j] << ", ";
>> 721 }
>> 722 std::cout << S.V[rng_get_N()-1];
572 std::cout << "}; ";
723 std::cout << "}; ";
573 std::cout << "counter= " << counter;
<< 724 std::cout << "counter= " << S.counter;
574 std::cout << "sumtot= " << sumtot << "\n";
<< 725 std::cout << "sumtot= " << S.sumtot << "\n";
575 }
726 }
576
727
577 MixMaxRng MixMaxRng::Branch()
728 MixMaxRng MixMaxRng::Branch()
578 {
729 {
579 sumtot = iterate_raw_vec(V, sumtot); counte << 730 S.sumtot = iterate_raw_vec(S.V.data(), S.sumtot); S.counter = 1;
580 MixMaxRng tmp=*this;
731 MixMaxRng tmp=*this;
581 tmp.BranchInplace(0); // daughter id
732 tmp.BranchInplace(0); // daughter id
582 return tmp;
733 return tmp;
583 }
734 }
584
735
585 void MixMaxRng::BranchInplace(int id)
736 void MixMaxRng::BranchInplace(int id)
586 {
737 {
587 // Dont forget to iterate the mother, when 738 // Dont forget to iterate the mother, when branching the daughter, or else will have collisions!
588 // a 64-bit LCG from Knuth line 26, is used 739 // a 64-bit LCG from Knuth line 26, is used to mangle a vector component
589 constexpr myuint_t MULT64=63641362238467930 740 constexpr myuint_t MULT64=6364136223846793005ULL;
590 myuint_t tmp=V[id];
<< 741 myuint_t tmp=S.V[id];
591 V[1] *= MULT64; V[id] &= M61;
<< 742 S.V[1] *= MULT64; S.V[id] &= M61;
592 sumtot = MIXMAX_MOD_MERSENNE( sumtot + V[id << 743 S.sumtot = MIXMAX_MOD_MERSENNE( S.sumtot + S.V[id] - tmp + M61);
593 sumtot = iterate_raw_vec(V, sumtot);
<< 744 S.sumtot = iterate_raw_vec(S.V.data(), S.sumtot);// printf("iterating!\n");
594 counter = 1;
<< 745 S.counter = 1;
595 }
746 }
596
747
597 } // namespace CLHEP
748 } // namespace CLHEP
598 749