1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
|
/*
* Tiny self-contained version of the PCG Random Number Generation for C++
* put together from pieces of the much larger C/C++ codebase.
* Wenzel Jakob, February 2015
*
* The PCG random number generator was developed by Melissa O'Neill <oneill@pcg-random.org>
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
* For additional information about the PCG random number generation scheme,
* including its license and other licensing options, visit
*
* http://www.pcg-random.org
*/
#ifndef __PCG32_H
#define __PCG32_H 1
#define PCG32_DEFAULT_STATE 0x853c49e6748fea9bULL
#define PCG32_DEFAULT_STREAM 0xda3e39cb94b95bdbULL
#define PCG32_MULT 0x5851f42d4c957f2dULL
#include <inttypes.h>
#include <cmath>
#include <cassert>
#include <algorithm>
/// PCG32 Pseudorandom number generator
struct pcg32 {
/// Initialize the pseudorandom number generator with default seed
pcg32() : state(PCG32_DEFAULT_STATE), inc(PCG32_DEFAULT_STREAM) {}
/// Initialize the pseudorandom number generator with the \ref seed() function
pcg32(uint64_t initstate, uint64_t initseq = 1u) { seed(initstate, initseq); }
/**
* \brief Seed the pseudorandom number generator
*
* Specified in two parts: a state initializer and a sequence selection
* constant (a.k.a. stream id)
*/
void seed(uint64_t initstate, uint64_t initseq = 1) {
state = 0U;
inc = (initseq << 1u) | 1u;
nextUInt();
state += initstate;
nextUInt();
}
/// Generate a uniformly distributed unsigned 32-bit random number
uint32_t nextUInt() {
uint64_t oldstate = state;
state = oldstate * PCG32_MULT + inc;
uint32_t xorshifted = (uint32_t) (((oldstate >> 18u) ^ oldstate) >> 27u);
uint32_t rot = (uint32_t) (oldstate >> 59u);
return (xorshifted >> rot) | (xorshifted << ((~rot + 1u) & 31));
}
/// Generate a uniformly distributed number, r, where 0 <= r < bound
uint32_t nextUInt(uint32_t bound) {
// To avoid bias, we need to make the range of the RNG a multiple of
// bound, which we do by dropping output less than a threshold.
// A naive scheme to calculate the threshold would be to do
//
// uint32_t threshold = 0x100000000ull % bound;
//
// but 64-bit div/mod is slower than 32-bit div/mod (especially on
// 32-bit platforms). In essence, we do
//
// uint32_t threshold = (0x100000000ull-bound) % bound;
//
// because this version will calculate the same modulus, but the LHS
// value is less than 2^32.
uint32_t threshold = (~bound+1u) % bound;
// Uniformity guarantees that this loop will terminate. In practice, it
// should usually terminate quickly; on average (assuming all bounds are
// equally likely), 82.25% of the time, we can expect it to require just
// one iteration. In the worst case, someone passes a bound of 2^31 + 1
// (i.e., 2147483649), which invalidates almost 50% of the range. In
// practice, bounds are typically small and only a tiny amount of the range
// is eliminated.
for (;;) {
uint32_t r = nextUInt();
if (r >= threshold)
return r % bound;
}
}
/// Generate a single precision floating point value on the interval [0, 1)
float nextFloat() {
/* Trick from MTGP: generate an uniformly distributed
single precision number in [1,2) and subtract 1. */
union {
uint32_t u;
float f;
} x;
x.u = (nextUInt() >> 9) | 0x3f800000UL;
return x.f - 1.0f;
}
/**
* \brief Generate a double precision floating point value on the interval [0, 1)
*
* \remark Since the underlying random number generator produces 32 bit output,
* only the first 32 mantissa bits will be filled (however, the resolution is still
* finer than in \ref nextFloat(), which only uses 23 mantissa bits)
*/
double nextDouble() {
/* Trick from MTGP: generate an uniformly distributed
double precision number in [1,2) and subtract 1. */
union {
uint64_t u;
double d;
} x;
x.u = ((uint64_t) nextUInt() << 20) | 0x3ff0000000000000ULL;
return x.d - 1.0;
}
/**
* \brief Multi-step advance function (jump-ahead, jump-back)
*
* The method used here is based on Brown, "Random Number Generation
* with Arbitrary Stride", Transactions of the American Nuclear
* Society (Nov. 1994). The algorithm is very similar to fast
* exponentiation.
*/
void advance(int64_t delta_) {
uint64_t
cur_mult = PCG32_MULT,
cur_plus = inc,
acc_mult = 1u,
acc_plus = 0u;
/* Even though delta is an unsigned integer, we can pass a signed
integer to go backwards, it just goes "the long way round". */
uint64_t delta = (uint64_t) delta_;
while (delta > 0) {
if (delta & 1) {
acc_mult *= cur_mult;
acc_plus = acc_plus * cur_mult + cur_plus;
}
cur_plus = (cur_mult + 1) * cur_plus;
cur_mult *= cur_mult;
delta /= 2;
}
state = acc_mult * state + acc_plus;
}
/**
* \brief Draw uniformly distributed permutation and permute the
* given STL container
*
* From: Knuth, TAoCP Vol. 2 (3rd 3d), Section 3.4.2
*/
template <typename Iterator> void shuffle(Iterator begin, Iterator end) {
if (begin <= end) return;
for (Iterator it = end - 1; it > begin; --it)
std::iter_swap(it, begin + nextUInt((uint32_t) (it - begin + 1)));
}
/// Compute the distance between two PCG32 pseudorandom number generators
int64_t operator-(const pcg32 &other) const {
assert(inc == other.inc);
uint64_t
cur_mult = PCG32_MULT,
cur_plus = inc,
cur_state = other.state,
the_bit = 1u,
distance = 0u;
while (state != cur_state) {
if ((state & the_bit) != (cur_state & the_bit)) {
cur_state = cur_state * cur_mult + cur_plus;
distance |= the_bit;
}
assert((state & the_bit) == (cur_state & the_bit));
the_bit <<= 1;
cur_plus = (cur_mult + 1ULL) * cur_plus;
cur_mult *= cur_mult;
}
return (int64_t) distance;
}
/// Equality operator
bool operator==(const pcg32 &other) const { return state == other.state && inc == other.inc; }
/// Inequality operator
bool operator!=(const pcg32 &other) const { return state != other.state || inc != other.inc; }
uint64_t state; // RNG state. All values are possible.
uint64_t inc; // Controls which RNG sequence (stream) is selected. Must *always* be odd.
};
#endif // __PCG32_H
|
