1// *Really* minimal PCG32 code / (c) 2014 M.E. O'Neill / pcg-random.org
2// Licensed under Apache License 2.0 (NO WARRANTY, etc. see website)
3
4#include "pcg.h"
5
6uint32_t pcg32_random_r(pcg32_random_t* rng)
7{
8 uint64_t oldstate = rng->state;
9 // Advance internal state
10 rng->state = oldstate * 6364136223846793005ULL + (rng->inc|1);
11 // Calculate output function (XSH RR), uses old state for max ILP
12 uint32_t xorshifted = ((oldstate >> 18u) ^ oldstate) >> 27u;
13 uint32_t rot = oldstate >> 59u;
14 return (xorshifted >> rot) | (xorshifted << ((-rot) & 31));
15}
16
17// Source from http://www.pcg-random.org/downloads/pcg-c-basic-0.9.zip
18void pcg32_srandom_r(pcg32_random_t* rng, uint64_t initstate, uint64_t initseq)
19{
20 rng->state = 0U;
21 rng->inc = (initseq << 1u) | 1u;
22 pcg32_random_r(rng);
23 rng->state += initstate;
24 pcg32_random_r(rng);
25}
26
27// Source from https://github.com/imneme/pcg-c-basic/blob/master/pcg_basic.c
28// pcg32_boundedrand_r(rng, bound):
29// Generate a uniformly distributed number, r, where 0 <= r < bound
30uint32_t pcg32_boundedrand_r(pcg32_random_t *rng, uint32_t bound) {
31 // To avoid bias, we need to make the range of the RNG a multiple of
32 // bound, which we do by dropping output less than a threshold.
33 // A naive scheme to calculate the threshold would be to do
34 //
35 // uint32_t threshold = 0x100000000ull % bound;
36 //
37 // but 64-bit div/mod is slower than 32-bit div/mod (especially on
38 // 32-bit platforms). In essence, we do
39 //
40 // uint32_t threshold = (0x100000000ull-bound) % bound;
41 //
42 // because this version will calculate the same modulus, but the LHS
43 // value is less than 2^32.
44 uint32_t threshold = -bound % bound;
45
46 // Uniformity guarantees that this loop will terminate. In practice, it
47 // should usually terminate quickly; on average (assuming all bounds are
48 // equally likely), 82.25% of the time, we can expect it to require just
49 // one iteration. In the worst case, someone passes a bound of 2^31 + 1
50 // (i.e., 2147483649), which invalidates almost 50% of the range. In
51 // practice, bounds are typically small and only a tiny amount of the range
52 // is eliminated.
53 for (;;) {
54 uint32_t r = pcg32_random_r(rng);
55 if (r >= threshold)
56 return r % bound;
57 }
58}
59