From d3ee4ed5d2a76d12e7c905b7aa7e1e7c87a49119 Mon Sep 17 00:00:00 2001 From: Tyge Lovset Date: Sat, 17 Jul 2021 12:19:04 +0200 Subject: Some refactoring. --- include/stc/crandom.h | 72 +++++++++++++++++++++++++++++++++------------------ 1 file changed, 47 insertions(+), 25 deletions(-) (limited to 'include/stc/crandom.h') diff --git a/include/stc/crandom.h b/include/stc/crandom.h index b6bd0c93..3e911d93 100644 --- a/include/stc/crandom.h +++ b/include/stc/crandom.h @@ -49,10 +49,10 @@ typedef struct {double lower, range;} stc64_uniformf_t; typedef struct {double mean, stddev, next; unsigned has_next;} stc64_normalf_t; /* PRNG stc64. - * Extremely fast PRNG suited for parallel usage with Weyl-sequence parameter. - * 320bit state, 256bit mutable. + * Very fast PRNG suited for parallel usage with Weyl-sequence parameter. + * 320-bit state, 256 bit is mutable. * Noticable faster than xoshiro and pcg, but slighly slower than wyrand64, - * which only has a single 64-bit state and is unfit when longer periods or + * which only has a single 64-bit state and therefore unfit when * multiple threads are required. * stc64 supports 2^63 unique threads with a minimum 2^64 period lengths each. * Passes PractRand tested up to 8TB output, Vigna's Hamming weight test, @@ -60,36 +60,37 @@ typedef struct {double mean, stddev, next; unsigned has_next;} stc64_normalf_t; * The 16-bit version with LR=6, RS=5, LS=3 passes PractRand to multiple TB input. */ -STC_API stc64_t stc64_init(uint64_t seed); -STC_API stc64_t stc64_with_seq(uint64_t seed, uint64_t seq); +/* Global STC64 PRNG */ +STC_API void stc64_srandom(uint64_t seed); +STC_API uint64_t stc64_random(void); +/* STC64 PRNG state */ +STC_API stc64_t stc64_init(uint64_t seed); +STC_API stc64_t stc64_with_seq(uint64_t seed, uint64_t seq); +/* Int64 uniform distributed RNG, range [low, high]. */ +STC_API stc64_uniform_t stc64_uniform_init(int64_t low, int64_t high); +/* Float64 uniform distributed RNG, range [low, high). */ +STC_API stc64_uniformf_t stc64_uniformf_init(double low, double high); +/* Normal distribution (double) */ +STC_API stc64_normalf_t stc64_normalf_init(double mean, double stddev); +STC_API double stc64_normalf(stc64_t* rng, stc64_normalf_t* dist); + STC_INLINE uint64_t stc64_rand(stc64_t* rng) { uint64_t *s = rng->state; enum {LR=24, RS=11, LS=3}; - const uint64_t result = (s[0] ^ (s[3] += s[4]|1)) + s[1]; + const uint64_t result = (s[0] ^ (s[3] += s[4])) + s[1]; s[0] = s[1] ^ (s[1] >> RS); s[1] = s[2] + (s[2] << LS); s[2] = ((s[2] << LR) | (s[2] >> (64 - LR))) + result; return result; } -/* Global random() */ -static stc64_t stc64_global = {{0x26aa069ea2fb1a4d, 0x70c72c95cd592d04, 0x504f333d3aa0b359, 0x26aa069ea2fb1a4d, 0x6a09e667a754166b}}; -STC_INLINE void stc64_srandom(uint64_t seed) { stc64_global = stc64_init(seed); } -STC_INLINE uint64_t stc64_random(void) { return stc64_rand(&stc64_global); } - /* Float64 random number in range [0.0, 1.0). */ STC_INLINE double stc64_randf(stc64_t* rng) { union {uint64_t i; double f;} u = {0x3FF0000000000000ull | (stc64_rand(rng) >> 12)}; return u.f - 1.0; } -/* Int64 uniform distributed RNG, range [low, high]. */ -STC_API stc64_uniform_t stc64_uniform_init(int64_t low, int64_t high); - /* Float64 uniform distributed RNG, range [low, high). */ -STC_INLINE stc64_uniformf_t stc64_uniformf_init(double low, double high) { - return c_make(stc64_uniformf_t){low, high - low}; -} STC_INLINE double stc64_uniformf(stc64_t* rng, stc64_uniformf_t* dist) { return stc64_randf(rng)*dist->range + dist->lower; } @@ -105,32 +106,53 @@ STC_INLINE int64_t stc64_uniform(stc64_t* rng, stc64_uniform_t* d) { return d->lower + hi; } -/* Normal distributed RNG, Float64. */ -STC_INLINE stc64_normalf_t stc64_normalf_init(double mean, double stddev) { - return c_make(stc64_normalf_t){mean, stddev, 0.0, 0}; -} -STC_API double stc64_normalf(stc64_t* rng, stc64_normalf_t* dist); - +/* -------------------------- IMPLEMENTATION ------------------------- */ #if !defined(STC_HEADER) || defined(STC_IMPLEMENTATION) +/* Global random() */ +static stc64_t stc64_global = {{ + 0x26aa069ea2fb1a4d, 0x70c72c95cd592d04, + 0x504f333d3aa0b359, 0x9e3779b97f4a7c15, + 0x6a09e667a754166b +}}; + +STC_DEF void stc64_srandom(uint64_t seed) { + stc64_global = stc64_init(seed); +} + +STC_DEF uint64_t stc64_random(void) { + return stc64_rand(&stc64_global); +} + STC_DEF stc64_t stc64_init(uint64_t seed) { return stc64_with_seq(seed, seed + 0x3504f333d3aa0b34); } + +/* rng.state[4] must be odd */ STC_DEF stc64_t stc64_with_seq(uint64_t seed, uint64_t seq) { stc64_t rng = {{seed, seed, seed, seed, (seq << 1) | 1}}; for (int i = 0; i < 12; ++i) stc64_rand(&rng); return rng; } -/* Very fast unbiased uniform int RNG with bounds [low, high] */ +/* Init unbiased uniform uint64 RNG with bounds [low, high] */ STC_DEF stc64_uniform_t stc64_uniform_init(int64_t low, int64_t high) { stc64_uniform_t dist = {low, (uint64_t) (high - low + 1)}; dist.threshold = (uint64_t)(-dist.range) % dist.range; return dist; } -/* Marsaglia polar method for gaussian/normal distribution. */ +/* Init uniform distributed float64 RNG, range [low, high). */ +STC_DEF stc64_uniformf_t stc64_uniformf_init(double low, double high) { + return c_make(stc64_uniformf_t){low, high - low}; +} + +/* Marsaglia polar method for gaussian/normal distribution, float64. */ +STC_DEF stc64_normalf_t stc64_normalf_init(double mean, double stddev) { + return c_make(stc64_normalf_t){mean, stddev, 0.0, 0}; +} + STC_DEF double stc64_normalf(stc64_t* rng, stc64_normalf_t* dist) { double u1, u2, s, m; if (dist->has_next++ & 1) -- cgit v1.2.3