Logo ROOT  
Reference Guide
RanluxppEngineImpl.cxx
Go to the documentation of this file.
1 // @(#)root/mathcore:$Id$
2 // Author: Jonas Hahnfeld 11/2020
3 
4 /*************************************************************************
5  * Copyright (C) 1995-2020, Rene Brun and Fons Rademakers. *
6  * All rights reserved. *
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
13 Implementation of the RANLUX++ generator
14 
15 RANLUX++ is an LCG equivalent of RANLUX using 576 bit numbers.
16 
17 Described in
18 A. Sibidanov, *A revision of the subtract-with-borrow random numbergenerators*,
19 *Computer Physics Communications*, 221(2017), 299-303,
20 preprint https://arxiv.org/pdf/1705.03123.pdf
21 
22 The code is loosely based on the Assembly implementation by A. Sibidanov
23 available at https://github.com/sibidanov/ranluxpp/.
24 */
25 
26 #include "Math/RanluxppEngine.h"
27 
28 #include "mulmod.h"
29 
30 #include <cassert>
31 #include <cstdint>
32 
33 namespace {
34 
35 // Variable templates are a feature of C++14, use the older technique of having
36 // a static member in a template class.
37 
38 template <int p>
39 struct RanluxppData;
40 
41 template <>
42 struct RanluxppData<24> {
43  static const uint64_t kA[9];
44 };
45 const uint64_t RanluxppData<24>::kA[] = {
46  0x0000000000000000, 0x0000000000000000, 0x0000000000010000, 0xfffe000000000000, 0xffffffffffffffff,
47  0xffffffffffffffff, 0xffffffffffffffff, 0xfffffffeffffffff, 0xffffffffffffffff,
48 };
49 
50 template <>
51 struct RanluxppData<2048> {
52  static const uint64_t kA[9];
53 };
54 const uint64_t RanluxppData<2048>::kA[] = {
55  0xed7faa90747aaad9, 0x4cec2c78af55c101, 0xe64dcb31c48228ec, 0x6d8a15a13bee7cb0, 0x20b2ca60cb78c509,
56  0x256c3d3c662ea36c, 0xff74e54107684ed2, 0x492edfcc0cc8e753, 0xb48c187cf5b22097,
57 };
58 
59 } // end anonymous namespace
60 
61 namespace ROOT {
62 namespace Math {
63 
64 template <int w, int p>
66 
67 private:
68  uint64_t fState[9]; ///< State of the generator
69  int fPosition = 0; ///< Current position in bits
70 
71  static constexpr const uint64_t *kA = RanluxppData<p>::kA;
72  static constexpr int kMaxPos = 9 * 64;
73 
74  /// Produce next block of random bits
75  void Advance()
76  {
77  mulmod(kA, fState);
78  fPosition = 0;
79  }
80 
81 public:
82  /// Return the next random bits, generate a new block if necessary
83  uint64_t NextRandomBits()
84  {
85  if (fPosition + w > kMaxPos) {
86  Advance();
87  }
88 
89  int idx = fPosition / 64;
90  int offset = fPosition % 64;
91  int numBits = 64 - offset;
92 
93  uint64_t bits = fState[idx] >> offset;
94  if (numBits < w) {
95  bits |= fState[idx + 1] << numBits;
96  }
97  bits &= ((uint64_t(1) << w) - 1);
98 
99  fPosition += w;
100  assert(fPosition <= kMaxPos && "position out of range!");
101 
102  return bits;
103  }
104 
105  /// Initialize and seed the state of the generator
106  void SetSeed(uint64_t s)
107  {
108  fState[0] = 1;
109  for (int i = 1; i < 9; i++) {
110  fState[i] = 0;
111  }
112 
113  uint64_t a_seed[9];
114  // Skip 2 ** 96 states.
115  powermod(kA, a_seed, uint64_t(1) << 48);
116  powermod(a_seed, a_seed, uint64_t(1) << 48);
117  // Skip another s states.
118  powermod(a_seed, a_seed, s);
119  mulmod(a_seed, fState);
120 
121  fPosition = 0;
122  }
123 
124  /// Skip `n` random numbers without generating them
125  void Skip(uint64_t n)
126  {
127  int left = (kMaxPos - fPosition) / w;
128  assert(left >= 0 && "position was out of range!");
129  if (n < (uint64_t)left) {
130  // Just skip the next few entries in the currently available bits.
131  fPosition += n * w;
132  assert(fPosition <= kMaxPos && "position out of range!");
133  return;
134  }
135 
136  n -= left;
137  // Need to advance and possibly skip over blocks.
138  int nPerState = kMaxPos / w;
139  int skip = (n / nPerState);
140 
141  uint64_t a_skip[9];
142  powermod(kA, a_skip, skip + 1);
143  mulmod(a_skip, fState);
144 
145  // Potentially skip numbers in the freshly generated block.
146  int remaining = n - skip * nPerState;
147  assert(remaining >= 0 && "should not end up at a negative position!");
148  fPosition = remaining * w;
149  assert(fPosition <= kMaxPos && "position out of range!");
150  }
151 };
152 
153 template <int p>
154 RanluxppEngine<p>::RanluxppEngine(uint64_t seed) : fImpl(new RanluxppEngineImpl<52, p>)
155 {
156  fImpl->SetSeed(seed);
157 }
158 
159 template <int p>
161 
162 template <int p>
164 {
165  return (*this)();
166 }
167 
168 template <int p>
170 {
171  // Get 52 bits of randomness.
172  uint64_t bits = fImpl->NextRandomBits();
173 
174  // Construct the double in [1, 2), using the random bits as mantissa.
175  static constexpr uint64_t exp = 0x3ff0000000000000;
176  union {
177  double dRandom;
178  uint64_t iRandom;
179  };
180  iRandom = exp | bits;
181 
182  // Shift to the right interval of [0, 1).
183  return dRandom - 1;
184 }
185 
186 template <int p>
188 {
189  return fImpl->NextRandomBits();
190 }
191 
192 template <int p>
193 void RanluxppEngine<p>::SetSeed(uint64_t seed)
194 {
195  fImpl->SetSeed(seed);
196 }
197 
198 template <int p>
200 {
201  fImpl->Skip(n);
202 }
203 
204 template class RanluxppEngine<24>;
205 template class RanluxppEngine<2048>;
206 
207 } // end namespace Math
208 } // end namespace ROOT
n
const Int_t n
Definition: legend1.C:16
ROOT::Math::RanluxppEngineImpl::kA
static constexpr const uint64_t * kA
Definition: RanluxppEngineImpl.cxx:71
ROOT::Math::RanluxppEngineImpl::fState
uint64_t fState[9]
State of the generator.
Definition: RanluxppEngineImpl.cxx:68
ROOT::Math::RanluxppEngineImpl::NextRandomBits
uint64_t NextRandomBits()
Return the next random bits, generate a new block if necessary.
Definition: RanluxppEngineImpl.cxx:83
exp
double exp(double)
TGeant4Unit::s
static constexpr double s
Definition: TGeant4SystemOfUnits.h:162
ROOT::Math::RanluxppEngineImpl::Advance
void Advance()
Produce next block of random bits.
Definition: RanluxppEngineImpl.cxx:75
ROOT::Math::RanluxppEngine::~RanluxppEngine
virtual ~RanluxppEngine()
ROOT::Math::RanluxppEngine::Skip
void Skip(uint64_t n)
Skip n random numbers without generating them.
Definition: RanluxppEngineImpl.cxx:199
ROOT::Math::RanluxppEngine::SetSeed
void SetSeed(uint64_t seed)
Initialize and seed the state of the generator.
Definition: RanluxppEngineImpl.cxx:193
ROOT::Math::RanluxppEngine::Rndm
double Rndm() override
Generate a double-precision random number with 52 bits of randomness.
Definition: RanluxppEngineImpl.cxx:163
mulmod.h
ROOT::Math::RanluxppEngineImpl::SetSeed
void SetSeed(uint64_t s)
Initialize and seed the state of the generator.
Definition: RanluxppEngineImpl.cxx:106
RanluxppEngine.h
ROOT::Math::RanluxppEngineImpl::kMaxPos
static constexpr int kMaxPos
Definition: RanluxppEngineImpl.cxx:72
ROOT::Math::RanluxppEngineImpl::fPosition
int fPosition
Current position in bits.
Definition: RanluxppEngineImpl.cxx:69
ROOT::Math::RanluxppEngine::IntRndm
uint64_t IntRndm()
Generate a random integer value with 52 bits.
Definition: RanluxppEngineImpl.cxx:187
ROOT::Math::RanluxppEngineImpl::Skip
void Skip(uint64_t n)
Skip n random numbers without generating them.
Definition: RanluxppEngineImpl.cxx:125
ROOT::Math::RanluxppEngine::fImpl
std::unique_ptr< RanluxppEngineImpl< 52, p > > fImpl
Definition: RanluxppEngine.h:30
ROOT::Math::RanluxppEngine::RanluxppEngine
RanluxppEngine(uint64_t seed=314159265)
Definition: RanluxppEngineImpl.cxx:154
ROOT::Math::RanluxppEngine::operator()
double operator()()
Generate a double-precision random number (non-virtual method)
Definition: RanluxppEngineImpl.cxx:169
ROOT::Math::RanluxppEngineImpl
Definition: RanluxppEngineImpl.cxx:65
ROOT
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...
Definition: EExecutionPolicy.hxx:4
Math
Namespace for new Math classes and functions.
ROOT::Math::RanluxppEngine
Implementation of the RANLUX++ generator.
Definition: RanluxppEngine.h:27