From b39c5fef6ebbe1b2db8ffc324679278382789c99 Mon Sep 17 00:00:00 2001 From: Tyge Løvset Date: Sun, 13 Feb 2022 22:14:49 +0100 Subject: Refactored crandom. Removed 32-bit version. Made normal dist inline - does not require math lib if not used. Added stc64_randomf() --- include/stc/crandom.h | 124 +++++++++++++++++--------------------------------- 1 file changed, 42 insertions(+), 82 deletions(-) (limited to 'include') diff --git a/include/stc/crandom.h b/include/stc/crandom.h index 7543f7b9..97e0275d 100644 --- a/include/stc/crandom.h +++ b/include/stc/crandom.h @@ -44,9 +44,7 @@ int main() { #include typedef struct stc64 { uint64_t state[5]; } stc64_t; -typedef struct stc32 { uint32_t state[5]; } stc32_t; typedef struct stc64_uniform { int64_t lower; uint64_t range, threshold; } stc64_uniform_t; -typedef struct stc32_uniform { int32_t lower; uint32_t range, threshold; } stc32_uniform_t; typedef struct stc64_uniformf { double lower, range; } stc64_uniformf_t; typedef struct stc64_normalf { double mean, stddev, next; unsigned has_next; } stc64_normalf_t; @@ -61,37 +59,28 @@ typedef struct stc64_normalf { double mean, stddev, next; unsigned has_next; } s * PractRand to multiple TB input. */ -/* Global STC64 PRNG */ +/* Global STC64 PRNGs */ STC_API void stc64_srandom(uint64_t seed); STC_API uint64_t stc64_random(void); - -STC_API uint64_t stc64_rand(stc64_t* rng); -STC_API uint32_t stc32_rand(stc32_t* rng); +STC_API double stc64_randomf(void); /* Init with sequence number */ -STC_API stc64_t stc64_with_seq(uint64_t seed, uint64_t seq); -STC_API stc32_t stc32_with_seq(uint32_t seed, uint32_t seq); - -/* Int uniform distributed RNG, range [low, high]. */ -STC_API stc64_uniform_t stc64_uniform_init(int64_t low, int64_t high); -STC_API stc32_uniform_t stc32_uniform_init(int32_t low, int32_t high); - -/* Float64 uniform distributed RNG, range [low, high). */ -STC_API stc64_uniformf_t stc64_uniformf_init(double low, double high); +STC_API stc64_t stc64_with_seq(uint64_t seed, uint64_t seq); -/* 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); - -/* Unbiased bounded uniform distribution. */ -STC_API int64_t stc64_uniform(stc64_t* rng, stc64_uniform_t* d); +/* Unbiased bounded uniform distribution. range [low, high] */ +STC_API int64_t stc64_uniform(stc64_t* rng, stc64_uniform_t* dist); STC_INLINE stc64_t stc64_init(uint64_t seed) { return stc64_with_seq(seed, seed + 0x3504f333d3aa0b37); } -STC_INLINE stc32_t stc32_init(uint32_t seed) { - return stc32_with_seq(seed, seed + 0xd3aa0b37); +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])) + 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; } /* Float64 random number in range [0.0, 1.0). */ @@ -105,10 +94,29 @@ STC_INLINE double stc64_uniformf(stc64_t* rng, stc64_uniformf_t* dist) { return stc64_randf(rng)*dist->range + dist->lower; } -STC_INLINE int32_t stc32_uniform(stc32_t* rng, stc32_uniform_t* d) { - uint64_t val; - do { val = stc32_rand(rng) * (uint64_t)d->range; } while ((uint32_t)val < d->threshold); - return d->lower + (val >> 32); +/* Init uniform distributed float64 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}; +} + +/* Marsaglia polar method for gaussian/normal distribution, float64. */ +STC_INLINE stc64_normalf_t stc64_normalf_init(double mean, double stddev) { + return c_make(stc64_normalf_t){mean, stddev, 0.0, 0}; +} + +/* Normal distribution PRNG */ +STC_INLINE double stc64_normalf(stc64_t* rng, stc64_normalf_t* dist) { + double u1, u2, s, m; + if (dist->has_next++ & 1) + return dist->next * dist->stddev + dist->mean; + do { + u1 = 2.0 * stc64_randf(rng) - 1.0; + u2 = 2.0 * stc64_randf(rng) - 1.0; + s = u1*u1 + u2*u2; + } while (s >= 1.0 || s == 0.0); + m = sqrt(-2.0 * log(s) / s); + dist->next = u2 * m; + return (u1 * m) * dist->stddev + dist->mean; } /* -------------------------- IMPLEMENTATION ------------------------- */ @@ -129,20 +137,18 @@ STC_DEF uint64_t stc64_random(void) { return stc64_rand(&stc64_global); } +STC_DEF double stc64_randomf(void) { + return stc64_randf(&stc64_global); +} + /* rng.state[4] must be odd */ STC_DEF stc64_t stc64_with_seq(uint64_t seed, uint64_t seq) { - stc64_t rng = {{seed, seed+0x26aa069ea2fb1a4d, seed+0x70c72c95cd592d04, - seed+0x504f333d3aa0b359, (seq << 1) | 1}}; + stc64_t rng = {{seed+0x26aa069ea2fb1a4d, seed+0x70c72c95cd592d04, + seed+0x504f333d3aa0b359, seed, seed<<1 | 1}}; for (int i = 0; i < 6; ++i) stc64_rand(&rng); return rng; } -STC_DEF stc32_t stc32_with_seq(uint32_t seed, uint32_t seq) { - stc32_t rng = {{seed, seed+0x26aa069e, seed+0xa2fb1a4d, - seed+0x70c72c95, (seq << 1) | 1}}; - for (int i = 0; i < 6; ++i) stc32_rand(&rng); - return rng; -} #ifdef _MSC_VER #pragma warning(disable: 4146) // unary minus operator applied to unsigned type #endif @@ -152,34 +158,11 @@ STC_DEF stc64_uniform_t stc64_uniform_init(int64_t low, int64_t high) { dist.threshold = (uint64_t)(-dist.range) % dist.range; return dist; } - -STC_DEF stc32_uniform_t stc32_uniform_init(int32_t low, int32_t high) { - stc32_uniform_t dist = {low, (uint32_t) (high - low + 1)}; - dist.threshold = (uint32_t)(-dist.range) % dist.range; - return dist; -} #ifdef _MSC_VER #pragma warning(default: 4146) #endif -STC_DEF 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])) + 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; -} - -STC_DEF uint32_t stc32_rand(stc32_t* rng) { - uint32_t *s = rng->state; enum {LR=21, RS=9, LS=3}; - const uint32_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] >> (32 - LR))) + result; - return result; -} - +/* Int uniform distributed RNG, range [low, high]. */ STC_DEF int64_t stc64_uniform(stc64_t* rng, stc64_uniform_t* d) { #ifdef c_umul128 uint64_t lo, hi; @@ -195,29 +178,6 @@ STC_DEF int64_t stc64_uniform(stc64_t* rng, stc64_uniform_t* d) { #endif } -/* 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) - return dist->next * dist->stddev + dist->mean; - do { - u1 = 2.0 * stc64_randf(rng) - 1.0; - u2 = 2.0 * stc64_randf(rng) - 1.0; - s = u1*u1 + u2*u2; - } while (s >= 1.0 || s == 0.0); - m = sqrt(-2.0 * log(s) / s); - dist->next = u2 * m; - return (u1 * m) * dist->stddev + dist->mean; -} #endif #endif #undef i_opt -- cgit v1.2.3