LCOV - code coverage report
Current view: top level - src/common - pg_prng.c (source / functions) Hit Total Coverage
Test: PostgreSQL 18devel Lines: 62 71 87.3 %
Date: 2024-07-19 16:11:31 Functions: 14 17 82.4 %
Legend: Lines: hit not hit

          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-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    18600486 : rotl(uint64 x, int bits)
      42             : {
      43    18600486 :     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     6200162 : xoroshiro128ss(pg_prng_state *state)
      55             : {
      56     6200162 :     uint64      s0 = state->s0,
      57     6200162 :                 sx = state->s1 ^ s0,
      58     6200162 :                 val = rotl(s0 * 5, 7) * 9;
      59             : 
      60             :     /* update state */
      61     6200162 :     state->s0 = rotl(s0, 24) ^ sx ^ (sx << 16);
      62     6200162 :     state->s1 = rotl(sx, 37);
      63             : 
      64     6200162 :     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       56176 : splitmix64(uint64 *state)
      73             : {
      74             :     /* state update */
      75       56176 :     uint64      val = (*state += UINT64CONST(0x9E3779B97f4A7C15));
      76             : 
      77             :     /* value extraction */
      78       56176 :     val = (val ^ (val >> 30)) * UINT64CONST(0xBF58476D1CE4E5B9);
      79       56176 :     val = (val ^ (val >> 27)) * UINT64CONST(0x94D049BB133111EB);
      80             : 
      81       56176 :     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       28088 : pg_prng_seed(pg_prng_state *state, uint64 seed)
      90             : {
      91       28088 :     state->s0 = splitmix64(&seed);
      92       28088 :     state->s1 = splitmix64(&seed);
      93             :     /* Let's just make sure we didn't get all-zeroes */
      94       28088 :     (void) pg_prng_seed_check(state);
      95       28088 : }
      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          14 : pg_prng_fseed(pg_prng_state *state, double fseed)
     103             : {
     104             :     /* Assume there's about 52 mantissa bits; the sign contributes too. */
     105          14 :     int64       seed = ((double) ((UINT64CONST(1) << 52) - 1)) * fseed;
     106             : 
     107          14 :     pg_prng_seed(state, (uint64) seed);
     108          14 : }
     109             : 
     110             : /*
     111             :  * Validate a PRNG seed value.
     112             :  */
     113             : bool
     114       60558 : 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       60558 :     if (unlikely(state->s0 == 0 && state->s1 == 0))
     121             :     {
     122           0 :         state->s0 = UINT64CONST(0x5851F42D4C957F2D);
     123           0 :         state->s1 = UINT64CONST(0x14057B7EF767814F);
     124             :     }
     125             : 
     126             :     /* As a convenience for the pg_prng_strong_seed macro, return true */
     127       60558 :     return true;
     128             : }
     129             : 
     130             : /*
     131             :  * Select a random uint64 uniformly from the range [0, PG_UINT64_MAX].
     132             :  */
     133             : uint64
     134      403116 : pg_prng_uint64(pg_prng_state *state)
     135             : {
     136      403116 :     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     2242380 : pg_prng_uint64_range(pg_prng_state *state, uint64 rmin, uint64 rmax)
     145             : {
     146             :     uint64      val;
     147             : 
     148     2242380 :     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     2242052 :         uint64      range = rmax - rmin;
     156     2242052 :         uint32      rshift = 63 - pg_leftmost_one_pos64(range);
     157             : 
     158             :         do
     159             :         {
     160     3463226 :             val = xoroshiro128ss(state) >> rshift;
     161     3463226 :         } while (val > range);
     162             :     }
     163             :     else
     164         328 :         val = 0;
     165             : 
     166     2242380 :     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           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       96192 : pg_prng_int64_range(pg_prng_state *state, int64 rmin, int64 rmax)
     193             : {
     194             :     int64       val;
     195             : 
     196       96192 :     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       96180 :         uval = (uint64) rmin +
     205       96180 :             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       96180 :         if (uval > PG_INT64_MAX)
     213       34574 :             val = (int64) (uval - PG_INT64_MIN) + PG_INT64_MIN;
     214             :         else
     215       61606 :             val = (int64) uval;
     216             :     }
     217             :     else
     218          12 :         val = rmin;
     219             : 
     220       96192 :     return val;
     221             : }
     222             : 
     223             : /*
     224             :  * Select a random uint32 uniformly from the range [0, PG_UINT32_MAX].
     225             :  */
     226             : uint32
     227       65876 : 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       65876 :     uint64      v = xoroshiro128ss(state);
     235             : 
     236       65876 :     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           0 : pg_prng_int32(pg_prng_state *state)
     244             : {
     245           0 :     uint64      v = xoroshiro128ss(state);
     246             : 
     247           0 :     return (int32) (v >> 32);
     248             : }
     249             : 
     250             : /*
     251             :  * Select a random int32 uniformly from the range [0, PG_INT32_MAX].
     252             :  */
     253             : int32
     254           2 : pg_prng_int32p(pg_prng_state *state)
     255             : {
     256           2 :     uint64      v = xoroshiro128ss(state);
     257             : 
     258           2 :     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     1439846 : pg_prng_double(pg_prng_state *state)
     269             : {
     270     1439846 :     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     1439846 :     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       25326 : 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       25326 :     u1 = 1.0 - pg_prng_double(state);
     302       25326 :     u2 = 1.0 - pg_prng_double(state);
     303             : 
     304             :     /* Apply Box-Muller transform to get one normal-valued output */
     305       25326 :     z0 = sqrt(-2.0 * log(u1)) * sin(2.0 * M_PI * u2);
     306       25326 :     return z0;
     307             : }
     308             : 
     309             : /*
     310             :  * Select a random boolean value.
     311             :  */
     312             : bool
     313      828096 : pg_prng_bool(pg_prng_state *state)
     314             : {
     315      828096 :     uint64      v = xoroshiro128ss(state);
     316             : 
     317      828096 :     return (bool) (v >> 63);
     318             : }

Generated by: LCOV version 1.14