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