From 724a087be71b03f7dd4047b3efd3ed091227f928 Mon Sep 17 00:00:00 2001 From: Tyge Løvset Date: Wed, 6 Jan 2021 11:57:16 +0100 Subject: Some rand additions. May swap names of stc64_rand(rng) and stc64_random(void) in future. --- benchmarks/crand_benchmark2.cpp | 25 ++++++++++++++++++++++++- docs/clist_api.md | 1 + docs/crand_api.md | 6 ++++++ stc/crand.h | 17 ++++++++++------- 4 files changed, 41 insertions(+), 8 deletions(-) diff --git a/benchmarks/crand_benchmark2.cpp b/benchmarks/crand_benchmark2.cpp index 966a3675..2692e760 100644 --- a/benchmarks/crand_benchmark2.cpp +++ b/benchmarks/crand_benchmark2.cpp @@ -4,6 +4,28 @@ #include "stc/crand.h" #include "others/pcg_random.hpp" +static struct stc32_state { stc64_t rng; uint64_t spare; unsigned n; } stc32_global = + {{0x7a5fed, 0x8e3f52, 0x9bc713, 0x6a09e667a7541669}, 0, 0}; + +STC_INLINE void stc32_srandom(uint64_t seed) { stc32_global.rng = stc64_init(seed); } +STC_INLINE uint32_t stc32_random(void) { + return (uint32_t) (++stc32_global.n & 1 ? (stc32_global.spare = stc64_rand(&stc32_global.rng)) + : (stc32_global.spare >> 32)); +} + +static unsigned long myrand_next = 1; + +/* RAND_MAX assumed to be 32767 */ +int myrand(void) { + myrand_next = myrand_next * 214013 + 2531011; + return (myrand_next >> 16) & 0x7fff; +} + +void mysrand(unsigned seed) { + myrand_next = seed; +} + + enum {N = 1000000000}; void test1(void) @@ -88,7 +110,8 @@ void test3(void) before = clock(); sum = 0; c_forrange (N) { - sum += stc64_rand(&rng); + //sum += stc64_rand(&rng); + sum += rand(); } diff = clock() - before; printf("stc64_random:\t\t%.02f, %zu sz:%zu\n", (float) diff / CLOCKS_PER_SEC, sum, sizeof rng); diff --git a/docs/clist_api.md b/docs/clist_api.md index deaec104..b056629d 100644 --- a/docs/clist_api.md +++ b/docs/clist_api.md @@ -95,6 +95,7 @@ clist_X_iter_t clist_X_begin(const clist_X* self); clist_X_iter_t clist_X_end(const clist_X* self); void clist_X_next(clist_X_iter_t* it); clist_X_value_t* clist_X_itval(clist_X_iter_t it); +clist_X_iter_t clist_X_fwd(clist_X_iter it, size_t n); clist_X_value_t clist_X_value_clone(clist_X_value_t val); ``` diff --git a/docs/crand_api.md b/docs/crand_api.md index 46d6ef1c..995c5fa0 100644 --- a/docs/crand_api.md +++ b/docs/crand_api.md @@ -40,14 +40,20 @@ All cstr definitions and prototypes may be included in your C source file by inc ## Methods ```c + void stc64_srandom(uint64_t seed); + uint64_t stc64_random(void); + 1) stc64_t stc64_init(uint64_t seed); 2) stc64_t stc64_with_seq(uint64_t seed, uint64_t seq); + 3) uint64_t stc64_rand(stc64_t* rng); 4) double stc64_randf(stc64_t* rng); + 5) stc64_uniform_t stc64_uniform_init(int64_t low, int64_t high); 6) int64_t stc64_uniform(stc64_t* rng, stc64_uniform_t* dist); 7) stc64_uniformf_t stc64_uniformf_init(double low, double high); 8) double stc64_uniformf(stc64_t* rng, stc64_uniformf_t* dist); + 9) stc64_normalf_t stc64_normalf_init(double mean, double stddev); 10) double stc64_normalf(stc64_t* rng, stc64_normalf_t* dist); ``` diff --git a/stc/crand.h b/stc/crand.h index 3d010eb2..929634ee 100644 --- a/stc/crand.h +++ b/stc/crand.h @@ -46,7 +46,7 @@ int main() { typedef struct {uint64_t state[4];} stc64_t; typedef struct {int64_t lower; uint64_t range, threshold;} stc64_uniform_t; typedef struct {double lower, range;} stc64_uniformf_t; -typedef struct {double mean, stddev, next; bool has_next;} stc64_normalf_t; +typedef struct {double mean, stddev, next; unsigned has_next;} stc64_normalf_t; /* Stc64: random number generator, range [0, 2^64). PRNG copyright Tyge Løvset, NORCE Research, 2020 */ @@ -61,6 +61,11 @@ STC_INLINE uint64_t stc64_rand(stc64_t* rng) { return result; } +/* Global random() */ +static stc64_t stc64_global = {{0x26aa069ea2fb1a4d, 0x70c72c95cd592d04, 0x504f333d3aa0b359, 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 [low, high). */ STC_INLINE double stc64_randf(stc64_t* rng) { union {uint64_t i; double f;} u = {0x3FF0000000000000ull | (stc64_rand(rng) >> 12)}; @@ -100,7 +105,7 @@ STC_INLINE int64_t stc64_uniform(stc64_t* rng, stc64_uniform_t* d) { /* Normal distributed RNG, Float64. */ STC_INLINE stc64_normalf_t stc64_normalf_init(double mean, double stddev) { - stc64_normalf_t dist = {mean, stddev, 0.0, false}; return dist; + stc64_normalf_t dist = {mean, stddev, 0.0, 0}; return dist; } STC_API double stc64_normalf(stc64_t* rng, stc64_normalf_t* dist); @@ -118,7 +123,7 @@ STC_API double stc64_normalf(stc64_t* rng, stc64_normalf_t* dist); */ STC_DEF stc64_t stc64_init(uint64_t seed) { - return stc64_with_seq(seed, 0x3504f333d3aa0b34); + return stc64_with_seq(seed, seed + 0x3504f333d3aa0b34); } STC_DEF stc64_t stc64_with_seq(uint64_t seed, uint64_t seq) { stc64_t rng = {{seed, seed, seed, (seq << 1u) | 1u}}; @@ -136,17 +141,15 @@ STC_DEF stc64_uniform_t stc64_uniform_init(int64_t low, int64_t high) { /* Marsaglia polar method for gaussian/normal distribution. */ STC_DEF double stc64_normalf(stc64_t* rng, stc64_normalf_t* dist) { double u1, u2, s, m; - if (dist->has_next) { - dist->has_next = false; + 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, dist->has_next = true; + dist->next = u2 * m; return (u1 * m) * dist->stddev + dist->mean; } -- cgit v1.2.3