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
|
/*
* PCG Random Number Generation for C.
*
* Copyright 2014 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
*/
/*
* This code is derived from the full C implementation, which is in turn
* derived from the canonical C++ PCG implementation. The C++ version
* has many additional features and is preferable if you can use C++ in
* your project.
*/
#include "pcg_basic.h"
// state for global RNGs
static pcg32_random_t pcg32_global = PCG32_INITIALIZER;
// pcg32_srandom(initstate, initseq)
// pcg32_srandom_r(rng, initstate, initseq):
// Seed the rng. Specified in two parts, state initializer and a
// sequence selection constant (a.k.a. stream id)
void pcg32_srandom_r(pcg32_random_t* rng, uint64_t initstate, uint64_t initseq)
{
rng->state = 0U;
rng->inc = (initseq << 1u) | 1u;
pcg32_random_r(rng);
rng->state += initstate;
pcg32_random_r(rng);
}
void pcg32_srandom(uint64_t seed, uint64_t seq)
{
pcg32_srandom_r(&pcg32_global, seed, seq);
}
// pcg32_random()
// pcg32_random_r(rng)
// Generate a uniformly distributed 32-bit random number
uint32_t pcg32_random_r(pcg32_random_t* rng)
{
uint64_t oldstate = rng->state;
rng->state = oldstate * 6364136223846793005ULL + rng->inc;
uint32_t xorshifted = (uint32_t)(((oldstate >> 18u) ^ oldstate) >> 27u);
uint32_t rot = oldstate >> 59u;
return (xorshifted >> rot) | (xorshifted << ((-rot) & 31));
}
uint32_t pcg32_random()
{
return pcg32_random_r(&pcg32_global);
}
// pcg32_boundedrand(bound):
// pcg32_boundedrand_r(rng, bound):
// Generate a uniformly distributed number, r, where 0 <= r < bound
uint32_t pcg32_boundedrand_r(pcg32_random_t* rng, 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 % 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 = pcg32_random_r(rng);
if (r >= threshold)
return r % bound;
}
}
uint32_t pcg32_boundedrand(uint32_t bound)
{
return pcg32_boundedrand_r(&pcg32_global, bound);
}
|
