PostgreSQL Source Code  git master
pg_prng.c
Go to the documentation of this file.
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-2024, 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 */
35 
36 
37 /*
38  * 64-bit rotate left
39  */
40 static inline uint64
41 rotl(uint64 x, int bits)
42 {
43  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
55 {
56  uint64 s0 = state->s0,
57  sx = state->s1 ^ s0,
58  val = rotl(s0 * 5, 7) * 9;
59 
60  /* update state */
61  state->s0 = rotl(s0, 24) ^ sx ^ (sx << 16);
62  state->s1 = rotl(sx, 37);
63 
64  return val;
65 }
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 splitmix64(uint64 *state)
73 {
74  /* state update */
75  uint64 val = (*state += UINT64CONST(0x9E3779B97f4A7C15));
76 
77  /* value extraction */
78  val = (val ^ (val >> 30)) * UINT64CONST(0xBF58476D1CE4E5B9);
79  val = (val ^ (val >> 27)) * UINT64CONST(0x94D049BB133111EB);
80 
81  return val ^ (val >> 31);
82 }
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
90 {
91  state->s0 = splitmix64(&seed);
92  state->s1 = splitmix64(&seed);
93  /* Let's just make sure we didn't get all-zeroes */
94  (void) pg_prng_seed_check(state);
95 }
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
103 {
104  /* Assume there's about 52 mantissa bits; the sign contributes too. */
105  int64 seed = ((double) ((UINT64CONST(1) << 52) - 1)) * fseed;
106 
107  pg_prng_seed(state, (uint64) seed);
108 }
109 
110 /*
111  * Validate a PRNG seed value.
112  */
113 bool
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  if (unlikely(state->s0 == 0 && state->s1 == 0))
121  {
122  state->s0 = UINT64CONST(0x5851F42D4C957F2D);
123  state->s1 = UINT64CONST(0x14057B7EF767814F);
124  }
125 
126  /* As a convenience for the pg_prng_strong_seed macro, return true */
127  return true;
128 }
129 
130 /*
131  * Select a random uint64 uniformly from the range [0, PG_UINT64_MAX].
132  */
133 uint64
135 {
136  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 pg_prng_uint64_range(pg_prng_state *state, uint64 rmin, uint64 rmax)
145 {
146  uint64 val;
147 
148  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  uint64 range = rmax - rmin;
156  uint32 rshift = 63 - pg_leftmost_one_pos64(range);
157 
158  do
159  {
160  val = xoroshiro128ss(state) >> rshift;
161  } while (val > range);
162  }
163  else
164  val = 0;
165 
166  return rmin + val;
167 }
168 
169 /*
170  * Select a random int64 uniformly from the range [PG_INT64_MIN, PG_INT64_MAX].
171  */
172 int64
174 {
175  return (int64) xoroshiro128ss(state);
176 }
177 
178 /*
179  * Select a random int64 uniformly from the range [0, PG_INT64_MAX].
180  */
181 int64
183 {
184  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 pg_prng_int64_range(pg_prng_state *state, int64 rmin, int64 rmax)
193 {
194  int64 val;
195 
196  if (likely(rmax > rmin))
197  {
198  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  uval = (uint64) rmin +
205  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  if (uval > PG_INT64_MAX)
213  val = (int64) (uval - PG_INT64_MIN) + PG_INT64_MIN;
214  else
215  val = (int64) uval;
216  }
217  else
218  val = rmin;
219 
220  return val;
221 }
222 
223 /*
224  * Select a random uint32 uniformly from the range [0, PG_UINT32_MAX].
225  */
226 uint32
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  uint64 v = xoroshiro128ss(state);
235 
236  return (uint32) (v >> 32);
237 }
238 
239 /*
240  * Select a random int32 uniformly from the range [PG_INT32_MIN, PG_INT32_MAX].
241  */
242 int32
244 {
245  uint64 v = xoroshiro128ss(state);
246 
247  return (int32) (v >> 32);
248 }
249 
250 /*
251  * Select a random int32 uniformly from the range [0, PG_INT32_MAX].
252  */
253 int32
255 {
256  uint64 v = xoroshiro128ss(state);
257 
258  return (int32) (v >> 33);
259 }
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
269 {
270  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  return ldexp((double) (v >> (64 - 52)), -52);
278 }
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
291 {
292  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  u1 = 1.0 - pg_prng_double(state);
302  u2 = 1.0 - pg_prng_double(state);
303 
304  /* Apply Box-Muller transform to get one normal-valued output */
305  z0 = sqrt(-2.0 * log(u1)) * sin(2.0 * M_PI * u2);
306  return z0;
307 }
308 
309 /*
310  * Select a random boolean value.
311  */
312 bool
314 {
315  uint64 v = xoroshiro128ss(state);
316 
317  return (bool) (v >> 63);
318 }
unsigned int uint32
Definition: c.h:509
#define likely(x)
Definition: c.h:313
signed int int32
Definition: c.h:497
#define PG_INT64_MAX
Definition: c.h:595
#define PG_INT64_MIN
Definition: c.h:594
#define unlikely(x)
Definition: c.h:314
long val
Definition: informix.c:689
int x
Definition: isn.c:71
static int pg_leftmost_one_pos64(uint64 word)
Definition: pg_bitutils.h:72
int64 pg_prng_int64_range(pg_prng_state *state, int64 rmin, int64 rmax)
Definition: pg_prng.c:192
double pg_prng_double(pg_prng_state *state)
Definition: pg_prng.c:268
uint64 pg_prng_uint64_range(pg_prng_state *state, uint64 rmin, uint64 rmax)
Definition: pg_prng.c:144
int32 pg_prng_int32(pg_prng_state *state)
Definition: pg_prng.c:243
static uint64 splitmix64(uint64 *state)
Definition: pg_prng.c:72
uint32 pg_prng_uint32(pg_prng_state *state)
Definition: pg_prng.c:227
uint64 pg_prng_uint64(pg_prng_state *state)
Definition: pg_prng.c:134
int32 pg_prng_int32p(pg_prng_state *state)
Definition: pg_prng.c:254
int64 pg_prng_int64(pg_prng_state *state)
Definition: pg_prng.c:173
void pg_prng_seed(pg_prng_state *state, uint64 seed)
Definition: pg_prng.c:89
pg_prng_state pg_global_prng_state
Definition: pg_prng.c:34
int64 pg_prng_int64p(pg_prng_state *state)
Definition: pg_prng.c:182
static uint64 rotl(uint64 x, int bits)
Definition: pg_prng.c:41
bool pg_prng_seed_check(pg_prng_state *state)
Definition: pg_prng.c:114
double pg_prng_double_normal(pg_prng_state *state)
Definition: pg_prng.c:290
#define M_PI
Definition: pg_prng.c:29
bool pg_prng_bool(pg_prng_state *state)
Definition: pg_prng.c:313
void pg_prng_fseed(pg_prng_state *state, double fseed)
Definition: pg_prng.c:102
static uint64 xoroshiro128ss(pg_prng_state *state)
Definition: pg_prng.c:54
static struct cvec * range(struct vars *v, chr a, chr b, int cases)
Definition: regc_locale.c:412
Definition: regguts.h:323