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