src/common/pg_prng.c

   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 */
  34 pg_prng_state pg_global_prng_state;
  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
  54 xoroshiro128ss(pg_prng_state *state)
  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
  89 pg_prng_seed(pg_prng_state *state, uint64 seed)
  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
 102 pg_prng_fseed(pg_prng_state *state, double fseed)
 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
 114 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         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
 134 pg_prng_uint64(pg_prng_state *state)
 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
 173 pg_prng_int64(pg_prng_state *state)
 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
 182 pg_prng_int64p(pg_prng_state *state)
 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
 227 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         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
 243 pg_prng_int32(pg_prng_state *state)
 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
 254 pg_prng_int32p(pg_prng_state *state)
 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
 268 pg_prng_double(pg_prng_state *state)
 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
 290 pg_prng_double_normal(pg_prng_state *state)
 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
 313 pg_prng_bool(pg_prng_state *state)
 314 {
 315         uint64          v = xoroshiro128ss(state);
 316
 317         return (bool) (v >> 63);
 318 }