Hello Tom,
Maybe on OpenBSD pg should switch srandom to srandom_deterministic?
Dunno. I'm fairly annoyed by their idea that they're smarter than POSIX.
Hmmm. I'm afraid that is not that hard.
However, for most of our uses of srandom, this behavior isn't awful;
it's only pgbench that has an expectation that the platform random()
can be made to behave deterministically. And TBH I think that's just
an expectation that's going to bite us.
I'd suggest that maybe we should get rid of the use of both random()
and srandom() in pgbench, and go over to letting set_random_seed()
fill the pg_erand48 state directly. In the integer-seed case you
could use something equivalent to pg_srand48. (In the other cases
probably you could do better, certainly the strong-random case could
just fill all 6 bytes directly.) That would get us to a place where
the behavior of --random-seed=N is not only deterministic but
platform-independent, which seems like an improvement.
That's a point. Althought I'm not found of round48, indeed having something
platform independent for testing makes definite sense.
I'll look into it.
Here is a POC which defines an internal interface for a PRNG, and use it
within pgbench, with several possible implementations which default to
rand48.
I must admit that I have a grudge against standard rand48:
- it is a known poor PRNG which was designed at a time when LCG where
basically the only low cost PRNG available. Newer designs were very
recent when the standard was set.
- it is a LCG, i.e. its low bits cycle quickly, so should not be used.
- so the 48 bit state size is relevant for generating 32 bits ints
and floats.
- however it eis used to generate more bits...
- the double function uses all 48 bits, whereas 52 need to be filled...
- and it is used to generate integers, which means that for large range
some values are inaccessible.
- 3 * 16 bits integers state looks silly on 32/64 bit architectures.
- ...
Given that postgres needs doubles (52 bits mantissa) and possibly 64 bits
integers, IMO the internal state should be 64 bits as a bare minimum,
which anyway is also the minimal bite on 64 bit architectures, which is
what is encoutered in practice.
--
Fabien.
diff --git a/src/bin/pgbench/pgbench.c b/src/bin/pgbench/pgbench.c
index b5c75ce1c6..06f1083d5a 100644
--- a/src/bin/pgbench/pgbench.c
+++ b/src/bin/pgbench/pgbench.c
@@ -185,7 +185,7 @@ int64 latency_limit = 0;
char *tablespace = NULL;
char *index_tablespace = NULL;
-/* random seed used when calling srandom() */
+/* random seed used to initialize the PRNG */
int64 random_seed = -1;
/*
@@ -279,14 +279,6 @@ typedef struct StatsData
SimpleStats lag;
} StatsData;
-/*
- * Struct to keep random state.
- */
-typedef struct RandomState
-{
- unsigned short xseed[3];
-} RandomState;
-
/*
* Connection state machine states.
*/
@@ -383,7 +375,7 @@ typedef struct
* Separate randomness for each client. This is used for random functions
* PGBENCH_RANDOM_* during the execution of the script.
*/
- RandomState cs_func_rs;
+ pg_prng_state cs_func_rs;
int use_file; /* index in sql_script for this client */
int command; /* command number in script */
@@ -450,9 +442,9 @@ typedef struct
* random state to make all of them independent of each other and therefore
* deterministic at the thread level.
*/
- RandomState ts_choose_rs; /* random state for selecting a script */
- RandomState ts_throttle_rs; /* random state for transaction throttling */
- RandomState ts_sample_rs; /* random state for log sampling */
+ pg_prng_state ts_choose_rs; /* random state for selecting a script */
+ pg_prng_state ts_throttle_rs; /* random state for transaction throttling */
+ pg_prng_state ts_sample_rs; /* random state for log sampling */
int64 throttle_trigger; /* previous/next throttling (us) */
FILE *logfile; /* where to log, or NULL */
@@ -832,30 +824,34 @@ strtodouble(const char *str, bool errorOK, double *dv)
}
/*
- * Initialize a random state struct.
+ * Initialize a pseudo-random state struct from main PRNG.
+ * This done at the beginning of process allows to have deterministic
+ * pgbench runs.
*/
static void
-initRandomState(RandomState *random_state)
+pg_prng_state_init(pg_prng_state *rs)
{
- random_state->xseed[0] = random();
- random_state->xseed[1] = random();
- random_state->xseed[2] = random();
+ pg_prng_setstate_r(rs, pg_prng_u64());
}
/* random number generator: uniform distribution from min to max inclusive */
static int64
-getrand(RandomState *random_state, int64 min, int64 max)
+getrand(pg_prng_state *random_state, int64 min, int64 max)
{
/*
* Odd coding is so that min and max have approximately the same chance of
* being selected as do numbers between them.
*
- * pg_erand48() is thread-safe and concurrent, which is why we use it
+ * pg_prng_double_r() is thread-safe and concurrent, which is why we use it
* rather than random(), which in glibc is non-reentrant, and therefore
* protected by a mutex, and therefore a bottleneck on machines with many
* CPUs.
+ *
+ * Note: pg_prng_double returns 52 bits (and pg_erand48 only 48 bits),
+ * so this does not really work if more is needed between min & max.
+ * Should rather use a modulo? Would be biased a little close to 64 bits.
*/
- return min + (int64) ((max - min + 1) * pg_erand48(random_state->xseed));
+ return min + (int64) ((max - min + 1) * pg_prng_double_r(random_state));
}
/*
@@ -864,7 +860,7 @@ getrand(RandomState *random_state, int64 min, int64 max)
* value is exp(-parameter).
*/
static int64
-getExponentialRand(RandomState *random_state, int64 min, int64 max,
+getExponentialRand(pg_prng_state *random_state, int64 min, int64 max,
double parameter)
{
double cut,
@@ -875,7 +871,7 @@ getExponentialRand(RandomState *random_state, int64 min, int64 max,
Assert(parameter > 0.0);
cut = exp(-parameter);
/* erand in [0, 1), uniform in (0, 1] */
- uniform = 1.0 - pg_erand48(random_state->xseed);
+ uniform = 1.0 - pg_prng_double_r(random_state);
/*
* inner expression in (cut, 1] (if parameter > 0), rand in [0, 1)
@@ -888,7 +884,7 @@ getExponentialRand(RandomState *random_state, int64 min, int64 max,
/* random number generator: gaussian distribution from min to max inclusive */
static int64
-getGaussianRand(RandomState *random_state, int64 min, int64 max,
+getGaussianRand(pg_prng_state *random_state, int64 min, int64 max,
double parameter)
{
double stdev;
@@ -912,13 +908,13 @@ getGaussianRand(RandomState *random_state, int64 min, int64 max,
do
{
/*
- * pg_erand48 generates [0,1), but for the basic version of the
+ * pseudo-random output in [0,1), but for the basic version of the
* Box-Muller transform the two uniformly distributed random numbers
* are expected in (0, 1] (see
* https://en.wikipedia.org/wiki/Box-Muller_transform)
*/
- double rand1 = 1.0 - pg_erand48(random_state->xseed);
- double rand2 = 1.0 - pg_erand48(random_state->xseed);
+ double rand1 = 1.0 - pg_prng_double_r(random_state);
+ double rand2 = 1.0 - pg_prng_double_r(random_state);
/* Box-Muller basic form transform */
double var_sqrt = sqrt(-2.0 * log(rand1));
@@ -948,7 +944,7 @@ getGaussianRand(RandomState *random_state, int64 min, int64 max,
* not be one.
*/
static int64
-getPoissonRand(RandomState *random_state, double center)
+getPoissonRand(pg_prng_state *random_state, double center)
{
/*
* Use inverse transform sampling to generate a value > 0, such that the
@@ -957,7 +953,7 @@ getPoissonRand(RandomState *random_state, double center)
double uniform;
/* erand in [0, 1), uniform in (0, 1] */
- uniform = 1.0 - pg_erand48(random_state->xseed);
+ uniform = 1.0 - pg_prng_double_r(random_state);
return (int64) (-log(uniform) * center + 0.5);
}
@@ -1035,7 +1031,7 @@ zipfFindOrCreateCacheCell(ZipfCache *cache, int64 n, double s)
* Luc Devroye, p. 550-551, Springer 1986.
*/
static int64
-computeIterativeZipfian(RandomState *random_state, int64 n, double s)
+computeIterativeZipfian(pg_prng_state *random_state, int64 n, double s)
{
double b = pow(2.0, s - 1.0);
double x,
@@ -1046,8 +1042,8 @@ computeIterativeZipfian(RandomState *random_state, int64 n, double s)
while (true)
{
/* random variates */
- u = pg_erand48(random_state->xseed);
- v = pg_erand48(random_state->xseed);
+ u = pg_prng_double_r(random_state);
+ v = pg_prng_double_r(random_state);
x = floor(pow(u, -1.0 / (s - 1.0)));
@@ -1065,11 +1061,11 @@ computeIterativeZipfian(RandomState *random_state, int64 n, double s)
* Jim Gray et al, SIGMOD 1994
*/
static int64
-computeHarmonicZipfian(ZipfCache *zipf_cache, RandomState *random_state,
+computeHarmonicZipfian(ZipfCache *zipf_cache, pg_prng_state *random_state,
int64 n, double s)
{
ZipfCell *cell = zipfFindOrCreateCacheCell(zipf_cache, n, s);
- double uniform = pg_erand48(random_state->xseed);
+ double uniform = pg_prng_double_r(random_state);
double uz = uniform * cell->harmonicn;
if (uz < 1.0)
@@ -1081,7 +1077,7 @@ computeHarmonicZipfian(ZipfCache *zipf_cache, RandomState *random_state,
/* random number generator: zipfian distribution from min to max inclusive */
static int64
-getZipfianRand(ZipfCache *zipf_cache, RandomState *random_state, int64 min,
+getZipfianRand(ZipfCache *zipf_cache, pg_prng_state *random_state, int64 min,
int64 max, double s)
{
int64 n = max - min + 1;
@@ -3565,7 +3561,7 @@ doLog(TState *thread, CState *st,
* to the random sample.
*/
if (sample_rate != 0.0 &&
- pg_erand48(thread->ts_sample_rs.xseed) > sample_rate)
+ pg_prng_double_r(&thread->ts_sample_rs) > sample_rate)
return;
/* should we aggregate the results or not? */
@@ -5126,12 +5122,11 @@ printResults(TState *threads, StatsData *total, instr_time total_time,
}
}
-/* call srandom based on some seed. NULL triggers the default behavior. */
+/* init prng based on some seed. NULL triggers the default behavior. */
static bool
set_random_seed(const char *seed)
{
- /* srandom expects an unsigned int */
- unsigned int iseed;
+ uint64 iseed;
if (seed == NULL || strcmp(seed, "time") == 0)
{
@@ -5139,7 +5134,7 @@ set_random_seed(const char *seed)
instr_time now;
INSTR_TIME_SET_CURRENT(now);
- iseed = (unsigned int) INSTR_TIME_GET_MICROSEC(now);
+ iseed = INSTR_TIME_GET_MICROSEC(now);
}
else if (strcmp(seed, "rand") == 0)
{
@@ -5155,7 +5150,7 @@ set_random_seed(const char *seed)
/* parse seed unsigned int value */
char garbage;
- if (sscanf(seed, "%u%c", &iseed, &garbage) != 1)
+ if (sscanf(seed, UINT64_FORMAT "%c", &iseed, &garbage) != 1)
{
fprintf(stderr,
"unrecognized random seed option \"%s\": expecting an unsigned integer, \"time\" or \"rand\"\n",
@@ -5165,9 +5160,12 @@ set_random_seed(const char *seed)
}
if (seed != NULL)
- fprintf(stderr, "setting random seed to %u\n", iseed);
- srandom(iseed);
- /* no precision loss: 32 bit unsigned int cast to 64 bit int */
+ fprintf(stderr, "setting random seed to " UINT64_FORMAT "\n", iseed);
+
+ /* initialize main PRNG */
+ pg_prng_setstate(iseed);
+
+ /* keep value, casted as signed */
random_seed = iseed;
return true;
}
@@ -5775,7 +5773,7 @@ main(int argc, char **argv)
for (i = 0; i < nclients; i++)
{
state[i].cstack = conditional_stack_create();
- initRandomState(&state[i].cs_func_rs);
+ pg_prng_state_init(&state[i].cs_func_rs);
}
if (debug)
@@ -5862,10 +5860,7 @@ main(int argc, char **argv)
/* set default seed for hash functions */
if (lookupVariable(&state[0], "default_seed") == NULL)
{
- uint64 seed = ((uint64) (random() & 0xFFFF) << 48) |
- ((uint64) (random() & 0xFFFF) << 32) |
- ((uint64) (random() & 0xFFFF) << 16) |
- (uint64) (random() & 0xFFFF);
+ uint64 seed = pg_prng_u64();
for (i = 0; i < nclients; i++)
if (!putVariableInt(&state[i], "startup", "default_seed", (int64) seed))
@@ -5909,9 +5904,9 @@ main(int argc, char **argv)
thread->state = &state[nclients_dealt];
thread->nstate =
(nclients - nclients_dealt + nthreads - i - 1) / (nthreads - i);
- initRandomState(&thread->ts_choose_rs);
- initRandomState(&thread->ts_throttle_rs);
- initRandomState(&thread->ts_sample_rs);
+ pg_prng_state_init(&thread->ts_choose_rs);
+ pg_prng_state_init(&thread->ts_throttle_rs);
+ pg_prng_state_init(&thread->ts_sample_rs);
thread->logfile = NULL; /* filled in later */
thread->latency_late = 0;
thread->zipf_cache.nb_cells = 0;
diff --git a/src/bin/pgbench/t/001_pgbench_with_server.pl b/src/bin/pgbench/t/001_pgbench_with_server.pl
index 17323e6ca5..6e5e2cf997 100644
--- a/src/bin/pgbench/t/001_pgbench_with_server.pl
+++ b/src/bin/pgbench/t/001_pgbench_with_server.pl
@@ -260,10 +260,10 @@ pgbench(
qr{setting random seed to 5432\b},
# After explicit seeding, the four * random checks (1-3,20) should be
- # deterministic, but not necessarily portable.
- qr{command=1.: int 1\d\b}, # uniform random: 12 on linux
- qr{command=2.: int 1\d\d\b}, # exponential random: 106 on linux
- qr{command=3.: int 1\d\d\d\b}, # gaussian random: 1462 on linux
+ # deterministic on all platforms
+ qr{command=1.: int 18\b}, # uniform random
+ qr{command=2.: int 101\b}, # exponential random
+ qr{command=3.: int 1509\b}, # gaussian random
qr{command=4.: int 4\b},
qr{command=5.: int 5\b},
qr{command=6.: int 6\b},
@@ -276,7 +276,7 @@ pgbench(
qr{command=15.: double 15\b},
qr{command=16.: double 16\b},
qr{command=17.: double 17\b},
- qr{command=20.: int \d\b}, # zipfian random: 1 on linux
+ qr{command=20.: int 3\b}, # zipfian random
qr{command=21.: double -27\b},
qr{command=22.: double 1024\b},
qr{command=23.: double 1\b},
diff --git a/src/include/port.h b/src/include/port.h
index a55c473262..7597f1a78b 100644
--- a/src/include/port.h
+++ b/src/include/port.h
@@ -348,6 +348,57 @@ extern long pg_lrand48(void);
extern long pg_jrand48(unsigned short xseed[3]);
extern void pg_srand48(long seed);
+/*----
+ * Postgres internal pseudo-random number generator (PRNG)
+ *
+ * 0: rand48 (pg historical linear congruential generator)
+ * 1: Knuth's MMIX LCG
+ * 2: splitmix64
+ * 3: xoshiro128+
+ * 4: xoshiro256** (default)
+ */
+#ifndef PG_PRNG_ALGORITHM
+#define PG_PRNG_ALGORITHM 0
+#endif
+
+/* set PRNG state byte size depending on algorithm */
+#if PG_PRNG_ALGORITHM == 0
+#define PG_PRNG_STATE_SIZE 6
+#elif PG_PRNG_ALGORITHM == 1 || PG_PRNG_ALGORITHM == 2
+#define PG_PRNG_STATE_SIZE 8
+#elif PG_PRNG_ALGORITHM == 3
+#define PG_PRNG_STATE_SIZE 16
+#elif PG_PRNG_ALGORITHM == 4
+#define PG_PRNG_STATE_SIZE 32
+#else
+#error "unexpected PG_PRNG_ALGORITHM value, must be in 0..4"
+#endif
+
+/* closed type seen as a byte array */
+typedef struct {
+ unsigned char data[PG_PRNG_STATE_SIZE];
+} pg_prng_state;
+
+/* interface with a reference to an external state */
+extern double pg_prng_double_r(pg_prng_state *state);
+extern int32 pg_prng_u31_r(pg_prng_state *state);
+extern int32 pg_prng_s32_r(pg_prng_state *state);
+extern uint32 pg_prng_u32_r(pg_prng_state *state);
+extern int64 pg_prng_s64_r(pg_prng_state *state);
+extern uint64 pg_prng_u64_r(pg_prng_state *state);
+extern void pg_prng_setstate_r(pg_prng_state *state, uint64 u);
+extern bool pg_prng_strong_setstate_r(pg_prng_state *state);
+
+/* interface with a shared internal state */
+extern double pg_prng_double(void);
+extern int32 pg_prng_u31(void);
+extern int32 pg_prng_s32(void);
+extern uint32 pg_prng_u32(void);
+extern int64 pg_prng_s64(void);
+extern uint64 pg_prng_u64(void);
+extern void pg_prng_setstate(uint64 u);
+extern bool pg_prng_strong_setstate(void);
+
#ifndef HAVE_FLS
extern int fls(int mask);
#endif
diff --git a/src/port/Makefile b/src/port/Makefile
index 9cfc0f9279..afa4161dde 100644
--- a/src/port/Makefile
+++ b/src/port/Makefile
@@ -37,7 +37,7 @@ LIBS += $(PTHREAD_LIBS)
OBJS = $(LIBOBJS) $(PG_CRC32C_OBJS) chklocale.o erand48.o inet_net_ntop.o \
noblock.o path.o pgcheckdir.o pgmkdirp.o pgsleep.o \
- pg_strong_random.o pgstrcasecmp.o pgstrsignal.o pqsignal.o \
+ pg_strong_random.o pgstrcasecmp.o pgstrsignal.o pqsignal.o prng.o \
qsort.o qsort_arg.o quotes.o snprintf.o sprompt.o strerror.o \
tar.o thread.o
diff --git a/src/port/prng.c b/src/port/prng.c
new file mode 100644
index 0000000000..4b8002270e
--- /dev/null
+++ b/src/port/prng.c
@@ -0,0 +1,438 @@
+/*-------------------------------------------------------------------------
+ *
+ * prng.c
+ *
+ * Postgres internal PRNG (pseudo-random number generator).
+ *
+ * This file supplies several implementations of a PRNG: POSIX-standard
+ * rand48, splitmix64, xoroshiro128+ and xoshiro256**, all wrapped into
+ * the same abstract interface so that changing the function only requires
+ * recompiling, without further changes. The variants with a reference to
+ * the state (*_r) are thread-safe if used with a thread-specific state.
+ *
+ * None of these functions should apply to anything which needs cryptographic
+ * security, even if seeded from strong random data: beware that the internal
+ * state can probably be recovered from a few outputs. If the state is too
+ * small or initialized from too small (random) data, they may be subject to a
+ * state enumeration attack. Most of these functions pass statistical tests,
+ * although they fail under some chosen transformation.
+ *
+ * Portions Copyright (c) 1996-2018, PostgreSQL Global Development Group
+ *
+ * Portions Copyright (c) 1993 Martin Birgmeier
+ * All rights reserved.
+ *
+ * You may redistribute unmodified or modified versions of this source
+ * code provided that the above copyright notice and this and the
+ * following conditions are retained.
+ *
+ * This software is provided ``as is'', and comes with no warranties
+ * of any kind. I shall in no event be liable for anything that happens
+ * to anyone/anything when using this software.
+ *
+ * IDENTIFICATION
+ * src/port/prng.c
+ *
+ *-------------------------------------------------------------------------
+ */
+
+#include "c.h"
+#include "port.h"
+
+#include <math.h>
+
+/*
+ * Left rotation, decent compilers should generate hardware instructions
+ */
+#define ROTL(x, n) ((x << n) | (x >> (64 - n)))
+
+/*
+ * PG PRNG ALGORITHM
+ */
+const char * pg_prng_algorithm =
+#if PG_PRNG_ALGORITHM == 0
+ /* POSIX standard */
+ "rand48"
+#elif PG_PRNG_ALGORITHM == 1
+ /* Donald Knuth's MMIX LCG */
+ "knuth-mmix-lcg"
+#elif PG_PRNG_ALGORITHM == 2
+ /*----
+ * Basic 8 bytes state PRNG from Java 8,
+ * see http://xoshiro.di.unimi.it/splitmix64.c
+ */
+ "splitmix64"
+#elif PG_PRNG_ALGORITHM == 3
+ /*----
+ * Algorithm by David Blackman and Sebastiano Vigna,
+ * see http://xoshiro.di.unimi.it/xoroshiro128plus.c
+ */
+ "xoroshiro128+"
+#elif PG_PRNG_ALGORITHM == 4
+ /*----
+ * Algorithm by David Blackman and Sebastiano Vigna,
+ * see http://xoshiro.di.unimi.it/xoshiro256starstar.c
+ */
+ "xoshiro256**"
+#endif
+ ;
+
+#if PG_PRNG_ALGORITHM >= 2
+/* also used to initialize large states */
+static uint64 _splitmix64(uint64 *state)
+{
+ uint64 v;
+ *state += UINT64CONST(0x9e3779b97f4a7c15);
+ v = *state;
+ v = (v ^ (v >> 30)) * UINT64CONST(0xbf58476d1ce4e5b9);
+ v = (v ^ (v >> 27)) * UINT64CONST(0x94d049bb133111eb);
+ return v ^ (v >> 31);
+}
+#endif
+
+/*
+ * INTERNAL STATE TYPE
+ */
+typedef struct {
+#if PG_PRNG_STATE_SIZE == 6
+ uint16 data[3];
+#else
+ uint64 data[PG_PRNG_STATE_SIZE / 8];
+#endif
+} _pg_prng_state;
+
+/*
+ * INTERNAL STATE DEFAULT
+ *
+ * pi in hexa is 3.243F6A8885 A308D31319 8A2E037073 44A4093822
+ * 299F31D008 2EFA98EC4E 6C89452821 E638D01377...
+ */
+static _pg_prng_state _prng_state = {
+ {
+
+#if PG_PRNG_STATE_SIZE == 6
+
+/* POSIX specifies 0x330e's use in srand48, but the other bits are arbitrary */
+#define RAND48_SEED_0 (0x330E)
+#define RAND48_SEED_1 (0x243F)
+#define RAND48_SEED_2 (0x6A88)
+
+ RAND48_SEED_0,
+ RAND48_SEED_1,
+ RAND48_SEED_2
+
+#elif PG_PRNG_STATE_SIZE == 8
+
+ UINT64CONST(0x243F6A8885A308D3)
+
+#elif PG_PRNG_STATE_SIZE == 16
+
+ UINT64CONST(0x243F6A8885A308D3),
+ UINT64CONST(0x13198A2E03707344)
+
+#else /* PG_PRNG_STATE_SIZE == 32 */
+
+ UINT64CONST(0x243F6A8885A308D3),
+ UINT64CONST(0x13198A2E03707344),
+ UINT64CONST(0xA4093822299F31D0),
+ UINT64CONST(0x082EFA98EC4E6C89)
+
+#endif
+
+ }
+};
+
+/*
+ * NEXT: advance internal state and return 64 bit value
+ */
+static uint64
+_next(_pg_prng_state *state)
+{
+#if PG_PRNG_ALGORITHM == 0
+
+/* These values are specified by POSIX */
+#define RAND48_MUL UINT64CONST(0x0005deece66d)
+#define RAND48_ADD UINT64CONST(0x000b)
+
+ /*
+ * We do the arithmetic in uint64; any type wider than 48 bits would work.
+ */
+ uint64 in;
+ uint64 out;
+
+ in = (uint64) state->data[2] << 32 |
+ (uint64) state->data[1] << 16 |
+ (uint64) state->data[0];
+
+ /* this is really 64 bits */
+ out = in * RAND48_MUL + RAND48_ADD;
+
+ /* truncated back to 48 bits */
+ state->data[0] = out & 0xFFFF;
+ state->data[1] = (out >> 16) & 0xFFFF;
+ state->data[2] = (out >> 32) & 0xFFFF;
+
+ return out;
+
+#elif PG_PRNG_ALGORITHM == 1
+
+ state->data[0] = UINT64CONST(6364136223846793005) * state->data[0] +
+ UINT64CONST(1442695040888963407);
+
+ return state->data[0];
+
+#elif PG_PRNG_ALGORITHM == 2
+
+ return _splitmix64(state->data);
+
+#elif PG_PRNG_ALGORITHM == 3
+
+ const uint64 s0 = state->data[0];
+ uint64 s1 = state->data[1];
+ const uint64 next = s0 + s1;
+
+ s1 ^= s0;
+
+ state->data[0] = ROTL(s0, 24) ^ s1 ^ (s1 << 16);
+ state->data[1] = ROTL(s1, 37);
+
+ return next;
+
+#else /* PG_PRNG_ALGORITHM == 4 */
+
+ const uint64 v = state->data[1] * 5;
+ const uint64 next = ROTL(v, 7) * 9;
+ const uint64 t = state->data[1] << 17;
+
+ state->data[2] ^= state->data[0];
+ state->data[3] ^= state->data[1];
+ state->data[1] ^= state->data[2];
+ state->data[0] ^= state->data[3];
+
+ state->data[2] ^= t;
+ state->data[3] = ROTL(state->data[3], 45);
+
+ return next;
+
+#endif
+}
+
+/*----
+ * Initialize the internal state from a 64-bit seed.
+ *
+ * This makes sense for determinism, eg needed under debug, where the user
+ * can supply an integer for this purpose. If this is not needed, consider
+ * using the strong variant.
+ */
+void
+pg_prng_setstate_r(pg_prng_state *state, uint64 u)
+{
+ _pg_prng_state * s = (_pg_prng_state *) state;
+
+ Assert(sizeof(pg_prng_state) == sizeof(_pg_prng_state));
+
+#if PG_PRNG_STATE_SIZE == 6
+
+ /*
+ * Per POSIX, this uses only 32 bits from "seed" even if "long" is wider.
+ * Hence, the set of possible seed values is smaller than it could be.
+ * Better practice is to use caller-supplied state and initialize it with
+ * random bits obtained from a high-quality source of random bits.
+ *
+ * Note: POSIX specifies a function seed48() that allows all 48 bits
+ * of the internal state to be set, but we don't currently support that.
+ */
+ s->data[0] = RAND48_SEED_0;
+ s->data[1] = (uint16) u;
+ s->data[2] = (uint16) (u >> 16);
+
+#elif PG_PRNG_STATE_SIZE == 8
+
+ s->data[0] = u;
+
+#elif PG_PRNG_STATE_SIZE == 16
+
+ s->data[0] = u;
+ s->data[1] = _splitmix64(&u);
+
+#else /* PG_PRNG_STATE_SIZE == 32 */
+
+ s->data[0] = u;
+ s->data[1] = _splitmix64(&u);
+ s->data[2] = _splitmix64(&u);
+ s->data[3] = _splitmix64(&u);
+
+#endif
+}
+
+/*------------------------------------------------------------------------------
+ * Common pseudo-random data extraction functions
+ */
+
+/*----
+ * Generate a pseudo-random floating-point value using caller-supplied state.
+ * Values are uniformly distributed over the interval [0.0, 1.0).
+ */
+double
+pg_prng_double_r(pg_prng_state *state)
+{
+ uint64 x = _next((_pg_prng_state *) state);
+
+#if PG_PRNG_STATE_SIZE == 6
+ /* not enough bits to fill the mantissa?
+ * in fact we have but that would depart from erand48
+ */
+ return ldexp((double) (x & UINT64CONST(0x0000FFFFFFFFFFFF)), -48);
+#else /* sizes 8, 16 and 32 are ok, prefer higher order bits though */
+ return ldexp((double) (x >> 12), -52);
+#endif
+}
+
+/* Generate a pseudo-random double in [0.0, 1.0) from internal state. */
+double
+pg_prng_double(void)
+{
+ return pg_prng_double_r((pg_prng_state *) &_prng_state);
+}
+
+/*----
+ * Generate a pseudo-random non-negative integral value using caller-supplied
+ * state. Values are uniformly distributed over the interval [0, 2^31).
+ */
+int32
+pg_prng_u31_r(pg_prng_state *state)
+{
+ uint64 x = _next((_pg_prng_state *) state);
+
+#if PG_PRNG_STATE_SIZE == 6
+ /* LCG must avoid low-weight bits */
+ return (x >> 17) & UINT64CONST(0x7FFFFFFF);
+#else /* combine available 64 bits */
+ return (x + (x >> 32)) & UINT64CONST(0x7FFFFFFF);
+#endif
+}
+
+/* Generate pseudo-random non-negative integral from internal state. */
+int32
+pg_prng_u31(void)
+{
+ return pg_prng_u31_r((pg_prng_state *) &_prng_state);
+}
+
+/*---
+ * Generate a pseudo-random signed integral value using caller-supplied state.
+ * Values are uniformly distributed over the interval [-2^31, 2^31).
+ */
+int32
+pg_prng_s32_r(pg_prng_state *state)
+{
+ uint64 x = _next((_pg_prng_state *) state);
+
+#if PG_PRNG_STATE_SIZE == 6
+ /* LCG must avoid low-weight bits */
+ return (int32) ((x >> 16) & UINT64CONST(0xFFFFFFFF));
+#else
+ /* combine available 64 bits */
+ return (int32) ((x + (x >> 32)) & UINT64CONST(0xFFFFFFFF));
+#endif
+}
+
+/* Generate pseudo-random signed integral from internal state. */
+int32
+pg_prng_s32(void)
+{
+ return pg_prng_s32_r((pg_prng_state *) &_prng_state);
+}
+
+/*---
+ * Generate a pseudo-random unsigned integral value using caller-supplied state.
+ * Values are uniformly distributed over the interval [0, 2^32).
+ */
+uint32
+pg_prng_u32_r(pg_prng_state *state)
+{
+ uint64 x = _next((_pg_prng_state *) state);
+
+#if PG_PRNG_STATE_SIZE == 6
+ /* LCG must avoid low-weight bits */
+ return (x >> 16) & UINT64CONST(0xFFFFFFFF);
+#else
+ /* combine available 64 bits */
+ return ((x + (x >> 32)) & UINT64CONST(0xFFFFFFFF));
+#endif
+}
+
+/* Generate pseudo-random signed integral from internal state. */
+uint32
+pg_prng_u32(void)
+{
+ return pg_prng_u32_r((pg_prng_state *) &_prng_state);
+}
+
+/*
+ * Generate a pseudo-random non-negative integral value using caller-supplied
+ * state. Values are uniformly distributed over the interval [-2^63, 2^63).
+ * (note: not really for low-end LCG)
+ */
+int64
+pg_prng_s64_r(pg_prng_state *state)
+{
+#if PG_PRNG_STATE_SIZE == 6 /* rand48 */
+ /* not enough bits, call _next twice and drop some low bits */
+ return (int64) (_next((_pg_prng_state *) state) << 16 ^
+ _next((_pg_prng_state *) state) >> 20);
+#else
+ return (int64) _next((_pg_prng_state *) state);
+#endif
+}
+
+/* Generate a pseudo-random int64 from internal state */
+int64
+pg_prng_s64(void)
+{
+ return pg_prng_s64_r((pg_prng_state *) &_prng_state);
+}
+
+/*
+ * Generate a pseudo-random non-negative integral value using caller-supplied
+ * state. Values are uniformly distributed over the interval [0, 2^64).
+ */
+uint64
+pg_prng_u64_r(pg_prng_state *state)
+{
+#if PG_PRNG_STATE_SIZE == 6 /* rand48 */
+ /* not enough bits, call _next twice and drop some low bits */
+ return _next((_pg_prng_state *) state) << 16 ^
+ _next((_pg_prng_state *) state) >> 20;
+#else
+ return _next((_pg_prng_state *) state);
+#endif
+}
+
+/* Generate a pseudo-random uint64 from internal state */
+uint64
+pg_prng_u64(void)
+{
+ return pg_prng_u64_r((pg_prng_state *) &_prng_state);
+}
+
+/* seed internal state from 64 bit integer */
+void
+pg_prng_setstate(uint64 u)
+{
+ pg_prng_setstate_r((pg_prng_state *) &_prng_state, u);
+}
+
+/* try to initialize prng state from a strong random source */
+bool
+pg_prng_strong_setstate_r(pg_prng_state *state)
+{
+ return pg_strong_random(&state->data, PG_PRNG_STATE_SIZE);
+}
+
+/* try to initialize the internal prng state from a strong random source */
+bool
+pg_prng_strong_setstate(void)
+{
+ return pg_prng_strong_setstate_r((pg_prng_state *) &_prng_state);
+}
