Branch data 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-2026, 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
41 : 2136078 : rotl(uint64 x, int bits)
42 : : {
43 : 2136078 : return (x << bits) | (x >> (64 - bits));
44 : : }
45 : :
46 : : /*
47 : : * 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
54 : 712026 : xoroshiro128ss(pg_prng_state *state)
55 : : {
56 : 1424052 : uint64 s0 = state->s0,
57 : 712026 : sx = state->s1 ^ s0,
58 : 712026 : val = rotl(s0 * 5, 7) * 9;
59 : :
60 : : /* update state */
61 : 712026 : state->s0 = rotl(s0, 24) ^ sx ^ (sx << 16);
62 : 712026 : state->s1 = rotl(sx, 37);
63 : :
64 : 1424052 : return val;
65 : 712026 : }
66 : :
67 : : /*
68 : : * We use this generator just to fill the xoroshiro128** state vector
69 : : * from a 64-bit seed.
70 : : */
71 : : static uint64
72 : 4018 : splitmix64(uint64 *state)
73 : : {
74 : : /* state update */
75 : 4018 : uint64 val = (*state += UINT64CONST(0x9E3779B97f4A7C15));
76 : :
77 : : /* value extraction */
78 : 4018 : val = (val ^ (val >> 30)) * UINT64CONST(0xBF58476D1CE4E5B9);
79 : 4018 : val = (val ^ (val >> 27)) * UINT64CONST(0x94D049BB133111EB);
80 : :
81 : 8036 : return val ^ (val >> 31);
82 : 4018 : }
83 : :
84 : : /*
85 : : * Initialize the PRNG state from a 64-bit integer,
86 : : * taking care that we don't produce all-zeroes.
87 : : */
88 : : void
89 : 2009 : pg_prng_seed(pg_prng_state *state, uint64 seed)
90 : : {
91 : 2009 : state->s0 = splitmix64(&seed);
92 : 2009 : state->s1 = splitmix64(&seed);
93 : : /* Let's just make sure we didn't get all-zeroes */
94 : 2009 : (void) pg_prng_seed_check(state);
95 : 2009 : }
96 : :
97 : : /*
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
102 : 8 : pg_prng_fseed(pg_prng_state *state, double fseed)
103 : : {
104 : : /* Assume there's about 52 mantissa bits; the sign contributes too. */
105 : 8 : int64 seed = ((double) ((UINT64CONST(1) << 52) - 1)) * fseed;
106 : :
107 : 8 : pg_prng_seed(state, (uint64) seed);
108 : 8 : }
109 : :
110 : : /*
111 : : * Validate a PRNG seed value.
112 : : */
113 : : bool
114 : 2843 : 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 [ + - + - ]: 2843 : if (unlikely(state->s0 == 0 && state->s1 == 0))
121 : : {
122 : 0 : state->s0 = UINT64CONST(0x5851F42D4C957F2D);
123 : 0 : state->s1 = UINT64CONST(0x14057B7EF767814F);
124 : 0 : }
125 : :
126 : : /* As a convenience for the pg_prng_strong_seed macro, return true */
127 : 2843 : return true;
128 : : }
129 : :
130 : : /*
131 : : * Select a random uint64 uniformly from the range [0, PG_UINT64_MAX].
132 : : */
133 : : uint64
134 : 0 : pg_prng_uint64(pg_prng_state *state)
135 : : {
136 : 0 : return xoroshiro128ss(state);
137 : : }
138 : :
139 : : /*
140 : : * Select a random uint64 uniformly from the range [rmin, rmax].
141 : : * If the range is empty, rmin is always produced.
142 : : */
143 : : uint64
144 : 48480 : pg_prng_uint64_range(pg_prng_state *state, uint64 rmin, uint64 rmax)
145 : : {
146 : 48480 : uint64 val;
147 : :
148 [ + + ]: 48480 : if (likely(rmax > rmin))
149 : : {
150 : : /*
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 : : */
155 : 48462 : uint64 range = rmax - rmin;
156 : 48462 : uint32 rshift = 63 - pg_leftmost_one_pos64(range);
157 : :
158 : 48462 : do
159 : : {
160 : 53764 : val = xoroshiro128ss(state) >> rshift;
161 [ + + ]: 53764 : } while (val > range);
162 : 48462 : }
163 : : else
164 : 18 : val = 0;
165 : :
166 : 96960 : return rmin + val;
167 : 48480 : }
168 : :
169 : : /*
170 : : * Select a random int64 uniformly from the range [PG_INT64_MIN, PG_INT64_MAX].
171 : : */
172 : : int64
173 : 0 : pg_prng_int64(pg_prng_state *state)
174 : : {
175 : 0 : return (int64) xoroshiro128ss(state);
176 : : }
177 : :
178 : : /*
179 : : * Select a random int64 uniformly from the range [0, PG_INT64_MAX].
180 : : */
181 : : int64
182 : 0 : pg_prng_int64p(pg_prng_state *state)
183 : : {
184 : 0 : return (int64) (xoroshiro128ss(state) & UINT64CONST(0x7FFFFFFFFFFFFFFF));
185 : : }
186 : :
187 : : /*
188 : : * Select a random int64 uniformly from the range [rmin, rmax].
189 : : * If the range is empty, rmin is always produced.
190 : : */
191 : : int64
192 : 16045 : pg_prng_int64_range(pg_prng_state *state, int64 rmin, int64 rmax)
193 : : {
194 : 16045 : int64 val;
195 : :
196 [ + + ]: 16045 : if (likely(rmax > rmin))
197 : : {
198 : 16040 : uint64 uval;
199 : :
200 : : /*
201 : : * Use pg_prng_uint64_range(). Can't simply pass it rmin and rmax,
202 : : * since (uint64) rmin will be larger than (uint64) rmax if rmin < 0.
203 : : */
204 : 32080 : uval = (uint64) rmin +
205 : 16040 : pg_prng_uint64_range(state, 0, (uint64) rmax - (uint64) rmin);
206 : :
207 : : /*
208 : : * Safely convert back to int64, avoiding implementation-defined
209 : : * behavior for values larger than PG_INT64_MAX. Modern compilers
210 : : * will reduce this to a simple assignment.
211 : : */
212 [ + + ]: 16040 : if (uval > PG_INT64_MAX)
213 : 5763 : val = (int64) (uval - PG_INT64_MIN) + PG_INT64_MIN;
214 : : else
215 : 10277 : val = (int64) uval;
216 : 16040 : }
217 : : else
218 : 5 : val = rmin;
219 : :
220 : 32090 : return val;
221 : 16045 : }
222 : :
223 : : /*
224 : : * Select a random uint32 uniformly from the range [0, PG_UINT32_MAX].
225 : : */
226 : : uint32
227 : 259645 : pg_prng_uint32(pg_prng_state *state)
228 : : {
229 : : /*
230 : : * Although xoroshiro128** is not known to have any weaknesses in
231 : : * randomness of low-order bits, we prefer to use the upper bits of its
232 : : * result here and below.
233 : : */
234 : 259645 : uint64 v = xoroshiro128ss(state);
235 : :
236 : 519290 : return (uint32) (v >> 32);
237 : 259645 : }
238 : :
239 : : /*
240 : : * Select a random int32 uniformly from the range [PG_INT32_MIN, PG_INT32_MAX].
241 : : */
242 : : int32
243 : 0 : pg_prng_int32(pg_prng_state *state)
244 : : {
245 : 0 : uint64 v = xoroshiro128ss(state);
246 : :
247 : 0 : return (int32) (v >> 32);
248 : 0 : }
249 : :
250 : : /*
251 : : * Select a random int32 uniformly from the range [0, PG_INT32_MAX].
252 : : */
253 : : int32
254 : 0 : pg_prng_int32p(pg_prng_state *state)
255 : : {
256 : 0 : uint64 v = xoroshiro128ss(state);
257 : :
258 : 0 : return (int32) (v >> 33);
259 : 0 : }
260 : :
261 : : /*
262 : : * Select a random double uniformly from the range [0.0, 1.0).
263 : : *
264 : : * Note: if you want a result in the range (0.0, 1.0], the standard way
265 : : * to get that is "1.0 - pg_prng_double(state)".
266 : : */
267 : : double
268 : 351824 : pg_prng_double(pg_prng_state *state)
269 : : {
270 : 351824 : uint64 v = xoroshiro128ss(state);
271 : :
272 : : /*
273 : : * As above, assume there's 52 mantissa bits in a double. This result
274 : : * could round to 1.0 if double's precision is less than that; but we
275 : : * assume IEEE float arithmetic elsewhere in Postgres, so this seems OK.
276 : : */
277 : 703648 : return ldexp((double) (v >> (64 - 52)), -52);
278 : 351824 : }
279 : :
280 : : /*
281 : : * Select a random double from the normal distribution with
282 : : * mean = 0.0 and stddev = 1.0.
283 : : *
284 : : * To get a result from a different normal distribution use
285 : : * STDDEV * pg_prng_double_normal() + MEAN
286 : : *
287 : : * Uses https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform
288 : : */
289 : : double
290 : 4220 : pg_prng_double_normal(pg_prng_state *state)
291 : : {
292 : 4220 : double u1,
293 : : u2,
294 : : z0;
295 : :
296 : : /*
297 : : * pg_prng_double generates [0, 1), but for the basic version of the
298 : : * Box-Muller transform the two uniformly distributed random numbers are
299 : : * expected to be in (0, 1]; in particular we'd better not compute log(0).
300 : : */
301 : 4220 : u1 = 1.0 - pg_prng_double(state);
302 : 4220 : u2 = 1.0 - pg_prng_double(state);
303 : :
304 : : /* Apply Box-Muller transform to get one normal-valued output */
305 : 4220 : z0 = sqrt(-2.0 * log(u1)) * sin(2.0 * M_PI * u2);
306 : 8440 : return z0;
307 : 4220 : }
308 : :
309 : : /*
310 : : * Select a random boolean value.
311 : : */
312 : : bool
313 : 46793 : pg_prng_bool(pg_prng_state *state)
314 : : {
315 : 46793 : uint64 v = xoroshiro128ss(state);
316 : :
317 : 93586 : return (bool) (v >> 63);
318 : 46793 : }
|
