Age Owner TLA Line data Source code
1 : /*-------------------------------------------------------------------------
2 : *
3 : * Pseudo-Random Number Generator
4 : *
5 : * We use Blackman and Vigna's xoroshiro128** 1.0 algorithm
6 : * to have a small, fast PRNG suitable for generating reasonably
7 : * good-quality 64-bit data. This should not be considered
8 : * cryptographically strong, however.
9 : *
10 : * About these generators: https://prng.di.unimi.it/
11 : * See also https://en.wikipedia.org/wiki/List_of_random_number_generators
12 : *
13 : * Copyright (c) 2021-2023, PostgreSQL Global Development Group
14 : *
15 : * src/common/pg_prng.c
16 : *
17 : *-------------------------------------------------------------------------
18 : */
19 :
20 : #include "c.h"
21 :
22 : #include <math.h>
23 :
24 : #include "common/pg_prng.h"
25 : #include "port/pg_bitutils.h"
26 :
27 : /* X/Open (XSI) requires <math.h> to provide M_PI, but core POSIX does not */
28 : #ifndef M_PI
29 : #define M_PI 3.14159265358979323846
30 : #endif
31 :
32 :
33 : /* process-wide state vector */
34 : pg_prng_state pg_global_prng_state;
35 :
36 :
37 : /*
38 : * 64-bit rotate left
39 : */
40 : static inline uint64
497 tgl 41 GIC 8188140 : rotl(uint64 x, int bits)
42 : {
43 8188140 : return (x << bits) | (x >> (64 - bits));
44 : }
45 :
46 : /*
497 tgl 47 ECB : * The basic xoroshiro128** algorithm.
48 : * Generates and returns a 64-bit uniformly distributed number,
49 : * updating the state vector for next time.
50 : *
51 : * Note: the state vector must not be all-zeroes, as that is a fixed point.
52 : */
53 : static uint64
497 tgl 54 GIC 2729380 : xoroshiro128ss(pg_prng_state *state)
55 : {
56 2729380 : uint64 s0 = state->s0,
57 2729380 : sx = state->s1 ^ s0,
58 2729380 : val = rotl(s0 * 5, 7) * 9;
59 :
497 tgl 60 ECB : /* update state */
497 tgl 61 GIC 2729380 : state->s0 = rotl(s0, 24) ^ sx ^ (sx << 16);
497 tgl 62 CBC 2729380 : state->s1 = rotl(sx, 37);
497 tgl 63 ECB :
497 tgl 64 CBC 2729380 : return val;
65 : }
66 :
497 tgl 67 ECB : /*
68 : * We use this generator just to fill the xoroshiro128** state vector
69 : * from a 64-bit seed.
70 : */
71 : static uint64
497 tgl 72 GIC 147286 : splitmix64(uint64 *state)
73 : {
74 : /* state update */
75 147286 : uint64 val = (*state += UINT64CONST(0x9E3779B97f4A7C15));
76 :
77 : /* value extraction */
497 tgl 78 CBC 147286 : val = (val ^ (val >> 30)) * UINT64CONST(0xBF58476D1CE4E5B9);
497 tgl 79 GIC 147286 : val = (val ^ (val >> 27)) * UINT64CONST(0x94D049BB133111EB);
80 :
497 tgl 81 CBC 147286 : return val ^ (val >> 31);
82 : }
83 :
497 tgl 84 ECB : /*
85 : * Initialize the PRNG state from a 64-bit integer,
86 : * taking care that we don't produce all-zeroes.
87 : */
88 : void
497 tgl 89 GIC 73643 : pg_prng_seed(pg_prng_state *state, uint64 seed)
90 : {
91 73643 : state->s0 = splitmix64(&seed);
92 73643 : state->s1 = splitmix64(&seed);
93 : /* Let's just make sure we didn't get all-zeroes */
94 73643 : (void) pg_prng_seed_check(state);
497 tgl 95 CBC 73643 : }
96 :
497 tgl 97 ECB : /*
98 : * Initialize the PRNG state from a double in the range [-1.0, 1.0],
99 : * taking care that we don't produce all-zeroes.
100 : */
101 : void
497 tgl 102 GIC 7 : pg_prng_fseed(pg_prng_state *state, double fseed)
103 : {
104 : /* Assume there's about 52 mantissa bits; the sign contributes too. */
105 7 : int64 seed = ((double) ((UINT64CONST(1) << 52) - 1)) * fseed;
106 :
107 7 : pg_prng_seed(state, (uint64) seed);
497 tgl 108 CBC 7 : }
109 :
110 : /*
497 tgl 111 ECB : * Validate a PRNG seed value.
112 : */
113 : bool
497 tgl 114 CBC 88495 : pg_prng_seed_check(pg_prng_state *state)
115 : {
116 : /*
117 : * If the seeding mechanism chanced to produce all-zeroes, insert
118 : * something nonzero. Anything would do; use Knuth's LCG parameters.
119 : */
120 88495 : if (unlikely(state->s0 == 0 && state->s1 == 0))
121 : {
497 tgl 122 UIC 0 : state->s0 = UINT64CONST(0x5851F42D4C957F2D);
123 0 : state->s1 = UINT64CONST(0x14057B7EF767814F);
124 : }
125 :
497 tgl 126 ECB : /* As a convenience for the pg_prng_strong_seed macro, return true */
497 tgl 127 GIC 88495 : return true;
497 tgl 128 EUB : }
129 :
130 : /*
131 : * Select a random uint64 uniformly from the range [0, PG_UINT64_MAX].
132 : */
497 tgl 133 ECB : uint64
497 tgl 134 GIC 1488 : pg_prng_uint64(pg_prng_state *state)
135 : {
136 1488 : return xoroshiro128ss(state);
137 : }
138 :
139 : /*
497 tgl 140 ECB : * Select a random uint64 uniformly from the range [rmin, rmax].
141 : * If the range is empty, rmin is always produced.
142 : */
143 : uint64
497 tgl 144 GIC 1046541 : pg_prng_uint64_range(pg_prng_state *state, uint64 rmin, uint64 rmax)
145 : {
146 : uint64 val;
147 :
148 1046541 : if (likely(rmax > rmin))
149 : {
497 tgl 150 ECB : /*
151 : * Use bitmask rejection method to generate an offset in 0..range.
152 : * Each generated val is less than twice "range", so on average we
153 : * should not have to iterate more than twice.
154 : */
497 tgl 155 GIC 1046377 : uint64 range = rmax - rmin;
156 1046377 : uint32 rshift = 63 - pg_leftmost_one_pos64(range);
157 :
158 : do
159 : {
160 1639607 : val = xoroshiro128ss(state) >> rshift;
497 tgl 161 CBC 1639607 : } while (val > range);
497 tgl 162 ECB : }
163 : else
497 tgl 164 GIC 164 : val = 0;
165 :
497 tgl 166 CBC 1046541 : return rmin + val;
497 tgl 167 ECB : }
168 :
169 : /*
170 : * Select a random int64 uniformly from the range [PG_INT64_MIN, PG_INT64_MAX].
171 : */
172 : int64
497 tgl 173 UIC 0 : pg_prng_int64(pg_prng_state *state)
174 : {
175 0 : return (int64) xoroshiro128ss(state);
176 : }
177 :
178 : /*
497 tgl 179 EUB : * Select a random int64 uniformly from the range [0, PG_INT64_MAX].
180 : */
181 : int64
497 tgl 182 UIC 0 : pg_prng_int64p(pg_prng_state *state)
183 : {
184 0 : return (int64) (xoroshiro128ss(state) & UINT64CONST(0x7FFFFFFFFFFFFFFF));
185 : }
186 :
187 : /*
497 tgl 188 EUB : * Select a random uint32 uniformly from the range [0, PG_UINT32_MAX].
189 : */
190 : uint32
497 tgl 191 GIC 68273 : pg_prng_uint32(pg_prng_state *state)
192 : {
193 : /*
194 : * Although xoroshiro128** is not known to have any weaknesses in
195 : * randomness of low-order bits, we prefer to use the upper bits of its
196 : * result here and below.
497 tgl 197 ECB : */
497 tgl 198 GIC 68273 : uint64 v = xoroshiro128ss(state);
199 :
200 68273 : return (uint32) (v >> 32);
201 : }
202 :
203 : /*
497 tgl 204 ECB : * Select a random int32 uniformly from the range [PG_INT32_MIN, PG_INT32_MAX].
205 : */
206 : int32
497 tgl 207 UIC 0 : pg_prng_int32(pg_prng_state *state)
208 : {
209 0 : uint64 v = xoroshiro128ss(state);
210 :
211 0 : return (int32) (v >> 32);
212 : }
497 tgl 213 EUB :
214 : /*
215 : * Select a random int32 uniformly from the range [0, PG_INT32_MAX].
216 : */
217 : int32
497 tgl 218 GIC 1 : pg_prng_int32p(pg_prng_state *state)
219 : {
220 1 : uint64 v = xoroshiro128ss(state);
221 :
222 1 : return (int32) (v >> 33);
223 : }
497 tgl 224 ECB :
225 : /*
226 : * Select a random double uniformly from the range [0.0, 1.0).
227 : *
228 : * Note: if you want a result in the range (0.0, 1.0], the standard way
229 : * to get that is "1.0 - pg_prng_double(state)".
230 : */
231 : double
497 tgl 232 GIC 664894 : pg_prng_double(pg_prng_state *state)
233 : {
234 664894 : uint64 v = xoroshiro128ss(state);
235 :
236 : /*
237 : * As above, assume there's 52 mantissa bits in a double. This result
497 tgl 238 ECB : * could round to 1.0 if double's precision is less than that; but we
239 : * assume IEEE float arithmetic elsewhere in Postgres, so this seems OK.
240 : */
497 tgl 241 GIC 664894 : return ldexp((double) (v >> (64 - 52)), -52);
242 : }
243 :
244 : /*
245 : * Select a random double from the normal distribution with
246 : * mean = 0.0 and stddev = 1.0.
247 : *
248 : * To get a result from a different normal distribution use
249 : * STDDEV * pg_prng_double_normal() + MEAN
250 : *
251 : * Uses https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform
252 : */
253 : double
90 tgl 254 GNC 6663 : pg_prng_double_normal(pg_prng_state *state)
255 : {
256 : double u1,
257 : u2,
258 : z0;
259 :
260 : /*
261 : * pg_prng_double generates [0, 1), but for the basic version of the
262 : * Box-Muller transform the two uniformly distributed random numbers are
263 : * expected to be in (0, 1]; in particular we'd better not compute log(0).
264 : */
265 6663 : u1 = 1.0 - pg_prng_double(state);
266 6663 : u2 = 1.0 - pg_prng_double(state);
267 :
268 : /* Apply Box-Muller transform to get one normal-valued output */
269 6663 : z0 = sqrt(-2.0 * log(u1)) * sin(2.0 * M_PI * u2);
270 6663 : return z0;
271 : }
272 :
273 : /*
274 : * Select a random boolean value.
275 : */
497 tgl 276 ECB : bool
497 tgl 277 GIC 355117 : pg_prng_bool(pg_prng_state *state)
278 : {
279 355117 : uint64 v = xoroshiro128ss(state);
280 :
281 355117 : return (bool) (v >> 63);
282 : }
|
