ROOT   Reference Guide
Searching...
No Matches
RanluxppEngineImpl.cxx
Go to the documentation of this file.
1// @(#)root/mathcore:$Id$
2// Author: Jonas Hahnfeld 11/2020
3
4/*************************************************************************
7 * *
8 * For the licensing terms see $ROOTSYS/LICENSE. * 9 * For the list of contributors see$ROOTSYS/README/CREDITS. *
10 *************************************************************************/
11
12/** \class ROOT::Math::RanluxppEngine
13Implementation of the RANLUX++ generator
14
15RANLUX++ is an LCG equivalent of RANLUX using 576 bit numbers.
16
17The idea of the generator (such as the initialization method) and the algorithm
18for the modulo operation are described in
19A. Sibidanov, *A revision of the subtract-with-borrow random numbergenerators*,
20*Computer Physics Communications*, 221(2017), 299-303,
21preprint https://arxiv.org/pdf/1705.03123.pdf
22
23The code is loosely based on the Assembly implementation by A. Sibidanov
24available at https://github.com/sibidanov/ranluxpp/.
25
26Compared to the original generator, this implementation contains a fix to ensure
27that the modulo operation of the LCG always returns the smallest value congruent
28to the modulus (based on notes by M. Lüscher). Also, the generator converts the
29LCG state back to RANLUX numbers (implementation based on notes by M. Lüscher).
30This avoids a bias in the generated numbers because the upper bits of the LCG
31state, that is smaller than the modulus \f$m = 2^{576} - 2^{240} + 1 \f$ (not
32a power of 2!), have a higher probability of being 0 than 1. And finally, this
33implementation draws 48 random bits for each generated floating point number
34(instead of 52 bits as in the original generator) to maintain the theoretical
35properties from understanding the original transition function of RANLUX as a
36chaotic dynamical system.
37*/
38
39#include "Math/RanluxppEngine.h"
40
41#include "ranluxpp/mulmod.h"
42#include "ranluxpp/ranlux_lcg.h"
43
44#include <cassert>
45#include <cstdint>
46
47namespace {
48
49// Variable templates are a feature of C++14, use the older technique of having
50// a static member in a template class.
51
52template <int p>
53struct RanluxppData;
54
55template <>
56struct RanluxppData<24> {
57 static const uint64_t kA[9];
58};
59const uint64_t RanluxppData<24>::kA[] = {
60 0x0000000000000000, 0x0000000000000000, 0x0000000000010000, 0xfffe000000000000, 0xffffffffffffffff,
61 0xffffffffffffffff, 0xffffffffffffffff, 0xfffffffeffffffff, 0xffffffffffffffff,
62};
63
64template <>
65struct RanluxppData<2048> {
66 static const uint64_t kA[9];
67};
68const uint64_t RanluxppData<2048>::kA[] = {
69 0xed7faa90747aaad9, 0x4cec2c78af55c101, 0xe64dcb31c48228ec, 0x6d8a15a13bee7cb0, 0x20b2ca60cb78c509,
70 0x256c3d3c662ea36c, 0xff74e54107684ed2, 0x492edfcc0cc8e753, 0xb48c187cf5b22097,
71};
72
73} // end anonymous namespace
74
75namespace ROOT {
76namespace Math {
77
78template <int w, int p>
80
81private:
82 uint64_t fState[9]; ///< RANLUX state of the generator
83 unsigned fCarry; ///< Carry bit of the RANLUX state
84 int fPosition = 0; ///< Current position in bits
85
86 static constexpr const uint64_t *kA = RanluxppData<p>::kA;
87 static constexpr int kMaxPos = 9 * 64;
88
89 /// Produce next block of random bits
91 {
92 uint64_t lcg[9];
93 to_lcg(fState, fCarry, lcg);
94 mulmod(kA, lcg);
95 to_ranlux(lcg, fState, fCarry);
96 fPosition = 0;
97 }
98
99public:
100 /// Return the next random bits, generate a new block if necessary
101 uint64_t NextRandomBits()
102 {
103 if (fPosition + w > kMaxPos) {
105 }
106
107 int idx = fPosition / 64;
108 int offset = fPosition % 64;
109 int numBits = 64 - offset;
110
111 uint64_t bits = fState[idx] >> offset;
112 if (numBits < w) {
113 bits |= fState[idx + 1] << numBits;
114 }
115 bits &= ((uint64_t(1) << w) - 1);
116
117 fPosition += w;
118 assert(fPosition <= kMaxPos && "position out of range!");
119
120 return bits;
121 }
122
123 /// Return a floating point number, converted from the next random bits.
125 {
126 static constexpr double div = 1.0 / (uint64_t(1) << w);
127 uint64_t bits = NextRandomBits();
128 return bits * div;
129 }
130
131 /// Initialize and seed the state of the generator
132 void SetSeed(uint64_t s)
133 {
134 uint64_t lcg[9];
135 lcg[0] = 1;
136 for (int i = 1; i < 9; i++) {
137 lcg[i] = 0;
138 }
139
140 uint64_t a_seed[9];
141 // Skip 2 ** 96 states.
142 powermod(kA, a_seed, uint64_t(1) << 48);
143 powermod(a_seed, a_seed, uint64_t(1) << 48);
144 // Skip another s states.
145 powermod(a_seed, a_seed, s);
146 mulmod(a_seed, lcg);
147
148 to_ranlux(lcg, fState, fCarry);
149 fPosition = 0;
150 }
151
152 /// Skip n random numbers without generating them
153 void Skip(uint64_t n)
154 {
155 int left = (kMaxPos - fPosition) / w;
156 assert(left >= 0 && "position was out of range!");
157 if (n < (uint64_t)left) {
158 // Just skip the next few entries in the currently available bits.
159 fPosition += n * w;
160 assert(fPosition <= kMaxPos && "position out of range!");
161 return;
162 }
163
164 n -= left;
165 // Need to advance and possibly skip over blocks.
166 int nPerState = kMaxPos / w;
167 int skip = (n / nPerState);
168
169 uint64_t a_skip[9];
170 powermod(kA, a_skip, skip + 1);
171
172 uint64_t lcg[9];
173 to_lcg(fState, fCarry, lcg);
174 mulmod(a_skip, lcg);
175 to_ranlux(lcg, fState, fCarry);
176
177 // Potentially skip numbers in the freshly generated block.
178 int remaining = n - skip * nPerState;
179 assert(remaining >= 0 && "should not end up at a negative position!");
180 fPosition = remaining * w;
181 assert(fPosition <= kMaxPos && "position out of range!");
182 }
183};
184
185template <int p>
187{
188 fImpl->SetSeed(seed);
189}
190
191template <int p>
193
194template <int p>
196{
197 return (*this)();
198}
199
200template <int p>
202{
203 return fImpl->NextRandomFloat();
204}
205
206template <int p>
208{
209 return fImpl->NextRandomBits();
210}
211
212template <int p>
213void RanluxppEngine<p>::SetSeed(uint64_t seed)
214{
215 fImpl->SetSeed(seed);
216}
217
218template <int p>
220{
221 fImpl->Skip(n);
222}
223
224template class RanluxppEngine<24>;
225template class RanluxppEngine<2048>;
226
227} // end namespace Math
228} // end namespace ROOT
unsigned fCarry
Carry bit of the RANLUX state.
uint64_t fState[9]
RANLUX state of the generator.
void Skip(uint64_t n)
Skip n random numbers without generating them.
static constexpr const uint64_t * kA
uint64_t NextRandomBits()
Return the next random bits, generate a new block if necessary.
void SetSeed(uint64_t s)
Initialize and seed the state of the generator.
int fPosition
Current position in bits.
double NextRandomFloat()
Return a floating point number, converted from the next random bits.
Produce next block of random bits.
Implementation of the RANLUX++ generator.
uint64_t IntRndm()
Generate a random integer value with 48 bits.
std::unique_ptr< RanluxppEngineImpl< 48, p > > fImpl
double Rndm() override
Generate a double-precision random number with 48 bits of randomness.
void Skip(uint64_t n)
Skip n random numbers without generating them.
void SetSeed(uint64_t seed)
Initialize and seed the state of the generator.
double operator()()
Generate a double-precision random number (non-virtual method)
RanluxppEngine(uint64_t seed=314159265)
const Int_t n
Definition legend1.C:16
static void mulmod(const uint64_t *in1, uint64_t *inout)
Combine multiply9x9 and mod_m with internal temporary storage.
Definition mulmod.h:215
static void powermod(const uint64_t *base, uint64_t *res, uint64_t n)
Compute base to the n modulo m.
Definition mulmod.h:229
Namespace for new Math classes and functions.
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...
static void to_lcg(const uint64_t *ranlux, unsigned c, uint64_t *lcg)
Convert RANLUX numbers to an LCG state.
Definition ranlux_lcg.h:26
static void to_ranlux(const uint64_t *lcg, uint64_t *ranlux, unsigned &c_out)
Convert an LCG state to RANLUX numbers.
Definition ranlux_lcg.h:58