| Geant4 Cross Reference |
1 // 1 //
2 // -*- C++ -*- 2 // -*- C++ -*-
3 // 3 //
4 // ------------------------------------------- 4 // -----------------------------------------------------------------------
5 // HEP Random 5 // HEP Random
6 // --- RanluxppEngine -- 6 // --- RanluxppEngine ---
7 // class implementation fi 7 // class implementation file
8 // ------------------------------------------- 8 // -----------------------------------------------------------------------
9 // Implementation of the RANLUX++ generator 9 // Implementation of the RANLUX++ generator
10 // 10 //
11 // RANLUX++ is an LCG equivalent of RANLUX usi 11 // RANLUX++ is an LCG equivalent of RANLUX using 576 bit numbers.
12 // 12 //
13 // The idea of the generator (such as the init 13 // The idea of the generator (such as the initialization method) and the algorithm
14 // for the modulo operation are described in 14 // for the modulo operation are described in
15 // A. Sibidanov, *A revision of the subtract-w 15 // A. Sibidanov, *A revision of the subtract-with-borrow random numbergenerators*,
16 // *Computer Physics Communications*, 221(2017 16 // *Computer Physics Communications*, 221(2017), 299-303,
17 // preprint https://arxiv.org/pdf/1705.03123.p 17 // preprint https://arxiv.org/pdf/1705.03123.pdf
18 // 18 //
19 // The code is loosely based on the Assembly i 19 // The code is loosely based on the Assembly implementation by A. Sibidanov
20 // available at https://github.com/sibidanov/r 20 // available at https://github.com/sibidanov/ranluxpp/.
21 // 21 //
22 // Compared to the original generator, this im 22 // Compared to the original generator, this implementation contains a fix to ensure
23 // that the modulo operation of the LCG always 23 // that the modulo operation of the LCG always returns the smallest value congruent
24 // to the modulus (based on notes by M. Lüsch 24 // to the modulus (based on notes by M. Lüscher). Also, the generator converts the
25 // LCG state back to RANLUX numbers (implement 25 // LCG state back to RANLUX numbers (implementation based on notes by M. Lüscher).
26 // This avoids a bias in the generated numbers 26 // This avoids a bias in the generated numbers because the upper bits of the LCG
27 // state, that is smaller than the modulus \f$ 27 // state, that is smaller than the modulus \f$ m = 2^{576} - 2^{240} + 1 \f$ (not
28 // a power of 2!), have a higher probability o 28 // a power of 2!), have a higher probability of being 0 than 1. And finally, this
29 // implementation draws 48 random bits for eac 29 // implementation draws 48 random bits for each generated floating point number
30 // (instead of 52 bits as in the original gene 30 // (instead of 52 bits as in the original generator) to maintain the theoretical
31 // properties from understanding the original 31 // properties from understanding the original transition function of RANLUX as a
32 // chaotic dynamical system. 32 // chaotic dynamical system.
33 // 33 //
34 // These modifications and the portable implem 34 // These modifications and the portable implementation in general are described in
35 // J. Hahnfeld, L. Moneta, *A Portable Impleme 35 // J. Hahnfeld, L. Moneta, *A Portable Implementation of RANLUX++*, vCHEP2021
36 // preprint https://arxiv.org/pdf/2106.02504.p 36 // preprint https://arxiv.org/pdf/2106.02504.pdf
37 37
38 #include "CLHEP/Random/RanluxppEngine.h" 38 #include "CLHEP/Random/RanluxppEngine.h"
39 39
40 #include "CLHEP/Random/engineIDulong.h" 40 #include "CLHEP/Random/engineIDulong.h"
41 #include "CLHEP/Utility/atomic_int.h" 41 #include "CLHEP/Utility/atomic_int.h"
42 42
43 #include "CLHEP/Random/ranluxpp/mulmod.h" 43 #include "CLHEP/Random/ranluxpp/mulmod.h"
44 #include "CLHEP/Random/ranluxpp/ranlux_lcg.h" 44 #include "CLHEP/Random/ranluxpp/ranlux_lcg.h"
45 45
46 #include <cassert> 46 #include <cassert>
47 #include <fstream> 47 #include <fstream>
48 #include <ios> 48 #include <ios>
49 #include <cstdint> <<
50 49
51 namespace CLHEP { 50 namespace CLHEP {
52 51
53 namespace { 52 namespace {
54 // Number of instances with automatic seed sel 53 // Number of instances with automatic seed selection.
55 CLHEP_ATOMIC_INT_TYPE numberOfEngines(0); 54 CLHEP_ATOMIC_INT_TYPE numberOfEngines(0);
56 55
57 const uint64_t kA_2048[] = { 56 const uint64_t kA_2048[] = {
58 0xed7faa90747aaad9, 0x4cec2c78af55c101, 0x 57 0xed7faa90747aaad9, 0x4cec2c78af55c101, 0xe64dcb31c48228ec,
59 0x6d8a15a13bee7cb0, 0x20b2ca60cb78c509, 0x 58 0x6d8a15a13bee7cb0, 0x20b2ca60cb78c509, 0x256c3d3c662ea36c,
60 0xff74e54107684ed2, 0x492edfcc0cc8e753, 0x 59 0xff74e54107684ed2, 0x492edfcc0cc8e753, 0xb48c187cf5b22097,
61 }; 60 };
62 } // namespace 61 } // namespace
63 62
64 RanluxppEngine::RanluxppEngine() : HepRandomEn 63 RanluxppEngine::RanluxppEngine() : HepRandomEngine() {
65 int numEngines = ++numberOfEngines; 64 int numEngines = ++numberOfEngines;
66 setSeed(numEngines); 65 setSeed(numEngines);
67 } 66 }
68 67
69 RanluxppEngine::RanluxppEngine(long seed) : He 68 RanluxppEngine::RanluxppEngine(long seed) : HepRandomEngine() {
70 theSeed = seed; 69 theSeed = seed;
71 setSeed(seed); 70 setSeed(seed);
72 } 71 }
73 72
74 RanluxppEngine::RanluxppEngine(std::istream &i 73 RanluxppEngine::RanluxppEngine(std::istream &is) : HepRandomEngine() {
75 get(is); 74 get(is);
76 } 75 }
77 76
78 RanluxppEngine::~RanluxppEngine() = default; 77 RanluxppEngine::~RanluxppEngine() = default;
79 78
80 static constexpr int kMaxPos = 9 * 64; 79 static constexpr int kMaxPos = 9 * 64;
81 static constexpr int kBits = 48; 80 static constexpr int kBits = 48;
82 81
83 void RanluxppEngine::advance() { 82 void RanluxppEngine::advance() {
84 uint64_t lcg[9]; 83 uint64_t lcg[9];
85 to_lcg(fState, fCarry, lcg); 84 to_lcg(fState, fCarry, lcg);
86 mulmod(kA_2048, lcg); 85 mulmod(kA_2048, lcg);
87 to_ranlux(lcg, fState, fCarry); 86 to_ranlux(lcg, fState, fCarry);
88 fPosition = 0; 87 fPosition = 0;
89 } 88 }
90 89
91 uint64_t RanluxppEngine::nextRandomBits() { 90 uint64_t RanluxppEngine::nextRandomBits() {
92 if (fPosition + kBits > kMaxPos) { 91 if (fPosition + kBits > kMaxPos) {
93 advance(); 92 advance();
94 } 93 }
95 94
96 int idx = fPosition / 64; 95 int idx = fPosition / 64;
97 int offset = fPosition % 64; 96 int offset = fPosition % 64;
98 int numBits = 64 - offset; 97 int numBits = 64 - offset;
99 98
100 uint64_t bits = fState[idx] >> offset; 99 uint64_t bits = fState[idx] >> offset;
101 if (numBits < kBits) { 100 if (numBits < kBits) {
102 bits |= fState[idx + 1] << numBits; 101 bits |= fState[idx + 1] << numBits;
103 } 102 }
104 bits &= ((uint64_t(1) << kBits) - 1); 103 bits &= ((uint64_t(1) << kBits) - 1);
105 104
106 fPosition += kBits; 105 fPosition += kBits;
107 assert(fPosition <= kMaxPos && "position out 106 assert(fPosition <= kMaxPos && "position out of range!");
108 107
109 return bits; 108 return bits;
110 } 109 }
111 110
112 double RanluxppEngine::flat() { 111 double RanluxppEngine::flat() {
113 // RandomEngine wants a "double random value 112 // RandomEngine wants a "double random values ranging between ]0,1[", so
114 // exclude all zero bits. 113 // exclude all zero bits.
115 uint64_t random; 114 uint64_t random;
116 do { 115 do {
117 random = nextRandomBits(); 116 random = nextRandomBits();
118 } while (random == 0); 117 } while (random == 0);
119 118
120 static constexpr double div = 1.0 / (uint64_ 119 static constexpr double div = 1.0 / (uint64_t(1) << kBits);
121 return random * div; 120 return random * div;
122 } 121 }
123 122
124 void RanluxppEngine::flatArray(const int size, 123 void RanluxppEngine::flatArray(const int size, double *vect) {
125 for (int i = 0; i < size; i++) { 124 for (int i = 0; i < size; i++) {
126 vect[i] = flat(); 125 vect[i] = flat();
127 } 126 }
128 } 127 }
129 128
130 void RanluxppEngine::setSeed(long seed, int) { 129 void RanluxppEngine::setSeed(long seed, int) {
131 theSeed = seed; 130 theSeed = seed;
132 131
133 uint64_t lcg[9]; 132 uint64_t lcg[9];
134 lcg[0] = 1; 133 lcg[0] = 1;
135 for (int i = 1; i < 9; i++) { 134 for (int i = 1; i < 9; i++) {
136 lcg[i] = 0; 135 lcg[i] = 0;
137 } 136 }
138 137
139 uint64_t a_seed[9]; 138 uint64_t a_seed[9];
140 // Skip 2 ** 96 states. 139 // Skip 2 ** 96 states.
141 powermod(kA_2048, a_seed, uint64_t(1) << 48) 140 powermod(kA_2048, a_seed, uint64_t(1) << 48);
142 powermod(a_seed, a_seed, uint64_t(1) << 48); 141 powermod(a_seed, a_seed, uint64_t(1) << 48);
143 // Skip more states according to seed. 142 // Skip more states according to seed.
144 powermod(a_seed, a_seed, seed); 143 powermod(a_seed, a_seed, seed);
145 mulmod(a_seed, lcg); 144 mulmod(a_seed, lcg);
146 145
147 to_ranlux(lcg, fState, fCarry); 146 to_ranlux(lcg, fState, fCarry);
148 fPosition = 0; 147 fPosition = 0;
149 } 148 }
150 149
151 void RanluxppEngine::setSeeds(const long *seed 150 void RanluxppEngine::setSeeds(const long *seeds, int) {
152 theSeeds = seeds; 151 theSeeds = seeds;
153 setSeed(*seeds, 0); 152 setSeed(*seeds, 0);
154 } 153 }
155 154
156 void RanluxppEngine::skip(uint64_t n) { 155 void RanluxppEngine::skip(uint64_t n) {
157 int left = (kMaxPos - fPosition) / kBits; 156 int left = (kMaxPos - fPosition) / kBits;
158 assert(left >= 0 && "position was out of ran 157 assert(left >= 0 && "position was out of range!");
159 if (n < (uint64_t)left) { 158 if (n < (uint64_t)left) {
160 // Just skip the next few entries in the c 159 // Just skip the next few entries in the currently available bits.
161 fPosition += n * kBits; 160 fPosition += n * kBits;
162 return; 161 return;
163 } 162 }
164 163
165 n -= left; 164 n -= left;
166 // Need to advance and possibly skip over bl 165 // Need to advance and possibly skip over blocks.
167 int nPerState = kMaxPos / kBits; 166 int nPerState = kMaxPos / kBits;
168 int skip = int(n / nPerState); << 167 int skip = (n / nPerState);
169 168
170 uint64_t a_skip[9]; 169 uint64_t a_skip[9];
171 powermod(kA_2048, a_skip, skip + 1); 170 powermod(kA_2048, a_skip, skip + 1);
172 171
173 uint64_t lcg[9]; 172 uint64_t lcg[9];
174 to_lcg(fState, fCarry, lcg); 173 to_lcg(fState, fCarry, lcg);
175 mulmod(a_skip, lcg); 174 mulmod(a_skip, lcg);
176 to_ranlux(lcg, fState, fCarry); 175 to_ranlux(lcg, fState, fCarry);
177 176
178 // Potentially skip numbers in the freshly g 177 // Potentially skip numbers in the freshly generated block.
179 int remaining = int(n - skip * nPerState); << 178 int remaining = n - skip * nPerState;
180 assert(remaining >= 0 && "should not end up 179 assert(remaining >= 0 && "should not end up at a negative position!");
181 fPosition = remaining * kBits; 180 fPosition = remaining * kBits;
182 assert(fPosition <= kMaxPos && "position out 181 assert(fPosition <= kMaxPos && "position out of range!");
183 } 182 }
184 183
185 void RanluxppEngine::saveStatus(const char fil 184 void RanluxppEngine::saveStatus(const char filename[]) const {
186 std::ofstream os(filename); 185 std::ofstream os(filename);
187 put(os); 186 put(os);
188 os.close(); 187 os.close();
189 } 188 }
190 189
191 void RanluxppEngine::restoreStatus(const char 190 void RanluxppEngine::restoreStatus(const char filename[]) {
192 std::ifstream is(filename); 191 std::ifstream is(filename);
193 get(is); 192 get(is);
194 is.close(); 193 is.close();
195 } 194 }
196 195
197 void RanluxppEngine::showStatus() const { 196 void RanluxppEngine::showStatus() const {
198 std::cout 197 std::cout
199 << "--------------------- RanluxppEngine 198 << "--------------------- RanluxppEngine status --------------------"
200 << std::endl; 199 << std::endl;
201 std::cout << " fState[] = {"; 200 std::cout << " fState[] = {";
202 std::cout << std::hex << std::setfill('0'); 201 std::cout << std::hex << std::setfill('0');
203 for (int i = 0; i < 9; i++) { 202 for (int i = 0; i < 9; i++) {
204 if (i % 3 == 0) { 203 if (i % 3 == 0) {
205 std::cout << std::endl << " "; 204 std::cout << std::endl << " ";
206 } else { 205 } else {
207 std::cout << " "; 206 std::cout << " ";
208 } 207 }
209 std::cout << "0x" << std::setw(16) << fSta 208 std::cout << "0x" << std::setw(16) << fState[i] << ",";
210 } 209 }
211 std::cout << std::endl << " }" << std::endl; 210 std::cout << std::endl << " }" << std::endl;
212 std::cout << std::dec; 211 std::cout << std::dec;
213 std::cout << " fCarry = " << fCarry << ", fP 212 std::cout << " fCarry = " << fCarry << ", fPosition = " << fPosition
214 << std::endl; 213 << std::endl;
215 std::cout 214 std::cout
216 << "------------------------------------ 215 << "----------------------------------------------------------------"
217 << std::endl; 216 << std::endl;
218 } 217 }
219 218
220 std::string RanluxppEngine::name() const { ret 219 std::string RanluxppEngine::name() const { return engineName(); }
221 220
222 std::string RanluxppEngine::engineName() { ret 221 std::string RanluxppEngine::engineName() { return "RanluxppEngine"; }
223 222
224 std::string RanluxppEngine::beginTag() { retur 223 std::string RanluxppEngine::beginTag() { return "RanluxppEngine-begin"; }
225 224
226 std::ostream &RanluxppEngine::put(std::ostream 225 std::ostream &RanluxppEngine::put(std::ostream &os) const {
227 os << beginTag() << "\n"; 226 os << beginTag() << "\n";
228 const std::vector<unsigned long> state = put 227 const std::vector<unsigned long> state = put();
229 for (unsigned long v : state) { 228 for (unsigned long v : state) {
230 os << v << "\n"; 229 os << v << "\n";
231 } 230 }
232 return os; 231 return os;
233 } 232 }
234 233
235 std::istream &RanluxppEngine::get(std::istream 234 std::istream &RanluxppEngine::get(std::istream &is) {
236 std::string tag; 235 std::string tag;
237 is >> tag; 236 is >> tag;
238 if (tag != beginTag()) { 237 if (tag != beginTag()) {
239 is.clear(std::ios::badbit | is.rdstate()); 238 is.clear(std::ios::badbit | is.rdstate());
240 std::cerr << "No RanluxppEngine found at c 239 std::cerr << "No RanluxppEngine found at current position\n";
241 return is; 240 return is;
242 } 241 }
243 return getState(is); 242 return getState(is);
244 } 243 }
245 244
246 std::istream &RanluxppEngine::getState(std::is 245 std::istream &RanluxppEngine::getState(std::istream &is) {
247 std::vector<unsigned long> state; 246 std::vector<unsigned long> state;
248 state.reserve(VECTOR_STATE_SIZE); 247 state.reserve(VECTOR_STATE_SIZE);
249 for (unsigned int i = 0; i < VECTOR_STATE_SI 248 for (unsigned int i = 0; i < VECTOR_STATE_SIZE; i++) {
250 unsigned long v; 249 unsigned long v;
251 is >> v; 250 is >> v;
252 state.push_back(v); 251 state.push_back(v);
253 } 252 }
254 253
255 getState(state); 254 getState(state);
256 return is; 255 return is;
257 } 256 }
258 257
259 std::vector<unsigned long> RanluxppEngine::put 258 std::vector<unsigned long> RanluxppEngine::put() const {
260 std::vector<unsigned long> v; 259 std::vector<unsigned long> v;
261 v.reserve(VECTOR_STATE_SIZE); 260 v.reserve(VECTOR_STATE_SIZE);
262 v.push_back(engineIDulong<RanluxppEngine>()) 261 v.push_back(engineIDulong<RanluxppEngine>());
263 262
264 // unsigned long is only guaranteed to be 32 263 // unsigned long is only guaranteed to be 32 bit wide, so chop up the 64 bit
265 // values in fState. 264 // values in fState.
266 for (int i = 0; i < 9; i++) { 265 for (int i = 0; i < 9; i++) {
267 unsigned long lower = static_cast<uint32_t 266 unsigned long lower = static_cast<uint32_t>(fState[i]);
268 v.push_back(lower); 267 v.push_back(lower);
269 unsigned long upper = static_cast<uint32_t 268 unsigned long upper = static_cast<uint32_t>(fState[i] >> 32);
270 v.push_back(upper); 269 v.push_back(upper);
271 } 270 }
272 271
273 v.push_back(fCarry); 272 v.push_back(fCarry);
274 v.push_back(fPosition); 273 v.push_back(fPosition);
275 return v; 274 return v;
276 } 275 }
277 276
278 bool RanluxppEngine::get(const std::vector<uns 277 bool RanluxppEngine::get(const std::vector<unsigned long> &v) {
279 if (v[0] != engineIDulong<RanluxppEngine>()) 278 if (v[0] != engineIDulong<RanluxppEngine>()) {
280 std::cerr << "RanluxppEngine::get(): " 279 std::cerr << "RanluxppEngine::get(): "
281 << "vector has wrong ID word - s 280 << "vector has wrong ID word - state unchanged" << std::endl;
282 return false; 281 return false;
283 } 282 }
284 return getState(v); 283 return getState(v);
285 } 284 }
286 285
287 bool RanluxppEngine::getState(const std::vecto 286 bool RanluxppEngine::getState(const std::vector<unsigned long> &v) {
288 if (v.size() != VECTOR_STATE_SIZE) { 287 if (v.size() != VECTOR_STATE_SIZE) {
289 std::cerr << "RanluxppEngine::getState(): 288 std::cerr << "RanluxppEngine::getState(): "
290 << "vector has wrong length - st 289 << "vector has wrong length - state unchanged" << std::endl;
291 return false; 290 return false;
292 } 291 }
293 292
294 // Assemble the state vector (see RanluxppEn 293 // Assemble the state vector (see RanluxppEngine::put).
295 for (int i = 0; i < 9; i++) { 294 for (int i = 0; i < 9; i++) {
296 uint64_t lower = v[2 * i + 1]; 295 uint64_t lower = v[2 * i + 1];
297 uint64_t upper = v[2 * i + 2]; 296 uint64_t upper = v[2 * i + 2];
298 fState[i] = (upper << 32) + lower; 297 fState[i] = (upper << 32) + lower;
299 } 298 }
300 fCarry = (unsigned int)v[19]; << 299 fCarry = v[19];
301 fPosition = (int)v[20]; << 300 fPosition = v[20];
302 301
303 return true; 302 return true;
304 } 303 }
305 304
306 } // namespace CLHEP 305 } // namespace CLHEP
307 306