Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
RanluxppEngine.cc
Go to the documentation of this file.
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
39
42
45
46#include <cassert>
47#include <fstream>
48#include <ios>
49#include <cstdint>
50
51namespace CLHEP {
52
53namespace {
54// Number of instances with automatic seed selection.
55CLHEP_ATOMIC_INT_TYPE numberOfEngines(0);
56
57const uint64_t kA_2048[] = {
58 0xed7faa90747aaad9, 0x4cec2c78af55c101, 0xe64dcb31c48228ec,
59 0x6d8a15a13bee7cb0, 0x20b2ca60cb78c509, 0x256c3d3c662ea36c,
60 0xff74e54107684ed2, 0x492edfcc0cc8e753, 0xb48c187cf5b22097,
61};
62} // namespace
63
65 int numEngines = ++numberOfEngines;
66 setSeed(numEngines);
67}
68
70 theSeed = seed;
71 setSeed(seed);
72}
73
75 get(is);
76}
77
79
80static constexpr int kMaxPos = 9 * 64;
81static constexpr int kBits = 48;
82
83void 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
91uint64_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
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
124void RanluxppEngine::flatArray(const int size, double *vect) {
125 for (int i = 0; i < size; i++) {
126 vect[i] = flat();
127 }
128}
129
130void 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
151void RanluxppEngine::setSeeds(const long *seeds, int) {
152 theSeeds = seeds;
153 setSeed(*seeds, 0);
154}
155
156void 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 = (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 = 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
185void RanluxppEngine::saveStatus(const char filename[]) const {
186 std::ofstream os(filename);
187 put(os);
188 os.close();
189}
190
191void RanluxppEngine::restoreStatus(const char filename[]) {
192 std::ifstream is(filename);
193 get(is);
194 is.close();
195}
196
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
220std::string RanluxppEngine::name() const { return engineName(); }
221
222std::string RanluxppEngine::engineName() { return "RanluxppEngine"; }
223
224std::string RanluxppEngine::beginTag() { return "RanluxppEngine-begin"; }
225
226std::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
235std::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
246std::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
259std::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
278bool 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
287bool 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 = v[19];
301 fPosition = v[20];
302
303 return true;
304}
305
306} // namespace CLHEP
#define CLHEP_ATOMIC_INT_TYPE
Definition: atomic_int.h:14
void showStatus() const override
void saveStatus(const char filename[]="Ranluxpp.conf") const override
std::string name() const override
void skip(uint64_t n)
std::istream & get(std::istream &is) override
static std::string engineName()
double flat() override
std::vector< unsigned long > put() const override
void setSeeds(const long *seeds, int dummy=0) override
void restoreStatus(const char filename[]="Ranluxpp.conf") override
static const unsigned int VECTOR_STATE_SIZE
void flatArray(const int size, double *vect) override
void setSeed(long seed, int dummy=0) override
std::istream & getState(std::istream &is) override
static std::string beginTag()
Definition: DoubConv.h:17
Definition: xmlparse.c:284