From 11a54d769b65b31a8a5870157079db193195ea00 Mon Sep 17 00:00:00 2001 From: Tyge Løvset Date: Wed, 5 Aug 2020 22:09:17 +0200 Subject: Some API changes in crandom. Added crandom_normal_f64() - normal distributed distribution. --- examples/list.c | 2 +- examples/priority.c | 2 +- examples/rngbirthday.c | 2 +- examples/rngtest.c | 6 ++-- stc/crandom.h | 80 ++++++++++++++++++++++++++++---------------------- 5 files changed, 51 insertions(+), 41 deletions(-) diff --git a/examples/list.c b/examples/list.c index b981b1ec..3018a1f1 100644 --- a/examples/list.c +++ b/examples/list.c @@ -8,7 +8,7 @@ int main() { int k, n = 100000; clist_fx list = clist_init; crandom_eng64_t eng = crandom_eng64_init(time(NULL)); - crandom_uniform_f64_t dist = crandom_uniform_f64_init(0.0f, n); + crandom_distrib_f64_t dist = crandom_uniform_f64_init(0.0f, n); for (int i = 0; i < 100000; ++i) clist_fx_push_back(&list, crandom_uniform_f64(&eng, dist)); diff --git a/examples/priority.c b/examples/priority.c index 70d89f2a..780e5809 100644 --- a/examples/priority.c +++ b/examples/priority.c @@ -10,7 +10,7 @@ declare_cvec_pqueue(i, >); // min-heap (increasing values) int main() { crandom_eng32_t pcg = crandom_eng32_init(time(NULL)); - crandom_uniform_i32_t dist = crandom_uniform_i32_init(0, 100000000); + crandom_distrib_i32_t dist = crandom_uniform_i32_init(0, 100000000); cvec_i heap = cvec_init; // Push ten million random numbers to priority queue diff --git a/examples/rngbirthday.c b/examples/rngbirthday.c index 595bacdb..22ecdfae 100644 --- a/examples/rngbirthday.c +++ b/examples/rngbirthday.c @@ -38,7 +38,7 @@ void distribution(void) const size_t N = 1ull << 28, M = 1ull << 9; // 1ull << 10; cmap_x map = cmap_x_with_capacity(M); clock_t now = clock(); - crandom_uniform_i32_t dist = crandom_uniform_i32_init(0, M); + crandom_distrib_i32_t dist = crandom_uniform_i32_init(0, M); for (size_t i = 0; i < N; ++i) { ++cmap_x_insert(&map, crandom_uniform_i32(&rng, dist), 0)->value; } diff --git a/examples/rngtest.c b/examples/rngtest.c index e19147c7..ca037292 100644 --- a/examples/rngtest.c +++ b/examples/rngtest.c @@ -32,15 +32,15 @@ int main(void) difference = clock() - before; printf("sfc64: %.02f, %zu\n", (float) difference / CLOCKS_PER_SEC, v); - crandom_uniform_i32_t i32dist = crandom_uniform_i32_init(10, 20); + crandom_distrib_i32_t i32dist = crandom_uniform_i32_init(10, 20); for (int i=0; i<8; ++i) printf("%d ", crandom_uniform_i32(&pcg, i32dist)); puts(""); - crandom_uniform_f32_t f32dist = crandom_uniform_f32_init(10, 20); + crandom_distrib_f32_t f32dist = crandom_uniform_f32_init(10, 20); for (int i=0; i<8; ++i) printf("%f ", crandom_uniform_f32(&pcg, f32dist)); puts(""); - crandom_uniform_f64_t fdist = crandom_uniform_f64_init(10, 20); + crandom_distrib_f64_t fdist = crandom_uniform_f64_init(10, 20); for (int i=0; i<8; ++i) printf("%f ", crandom_uniform_f64(&sfc, fdist)); puts(""); } \ No newline at end of file diff --git a/stc/crandom.h b/stc/crandom.h index 5d547bed..b73dd11b 100644 --- a/stc/crandom.h +++ b/stc/crandom.h @@ -26,10 +26,11 @@ #include "cdefs.h" #include +#include /* crandom_eng32_t eng = crandom_eng32_init(seed); - crandom_uniform_f32_t fdist = crandom_uniform_f32_init(1.0f, 6.0f); - crandom_uniform_i32_t idist = crandom_uniform_i32_init(1, 6); + crandom_distrib_f32_t fdist = crandom_uniform_f32_init(1.0f, 6.0f); + crandom_distrib_i32_t idist = crandom_uniform_i32_init(1, 6); uint32_t i = crandom_i32(&eng); int j = crandom_uniform_i32(&eng, idist); @@ -37,8 +38,8 @@ */ typedef struct {uint64_t state[2];} crandom_eng32_t; -typedef struct {int32_t min, range;} crandom_uniform_i32_t; -typedef struct {float min, range;} crandom_uniform_f32_t; +typedef struct {int32_t offset, range;} crandom_distrib_i32_t; +typedef struct {float offset, range;} crandom_distrib_f32_t; /* 32 bit random number generator engine */ STC_API crandom_eng32_t crandom_eng32_with_seq(uint64_t seed, uint64_t seq); @@ -55,24 +56,24 @@ STC_INLINE float crandom_f32(crandom_eng32_t* rng) { } /* int random number generator in range [low, high] */ -STC_INLINE crandom_uniform_i32_t crandom_uniform_i32_init(int32_t low, int32_t high) { - crandom_uniform_i32_t dist = {low, high - low + 1}; return dist; +STC_INLINE crandom_distrib_i32_t crandom_uniform_i32_init(int32_t low, int32_t high) { + crandom_distrib_i32_t dist = {low, high - low + 1}; return dist; } -STC_INLINE int32_t crandom_uniform_i32(crandom_eng32_t* rng, crandom_uniform_i32_t dist) { - return dist.min + (int32_t) (((uint64_t) crandom_i32(rng) * dist.range) >> 32); +STC_INLINE int32_t crandom_uniform_i32(crandom_eng32_t* rng, crandom_distrib_i32_t dist) { + return dist.offset + (int32_t) (((uint64_t) crandom_i32(rng) * dist.range) >> 32); } /* float random number in range [low, high). Note: 23 bit resolution. */ -STC_INLINE crandom_uniform_f32_t crandom_uniform_f32_init(float low, float high) { - crandom_uniform_f32_t dist = {low, high - low}; return dist; +STC_INLINE crandom_distrib_f32_t crandom_uniform_f32_init(float low, float high) { + crandom_distrib_f32_t dist = {low, high - low}; return dist; } -STC_INLINE float crandom_uniform_f32(crandom_eng32_t* rng, crandom_uniform_f32_t dist) { - return dist.min + crandom_f32(rng) * dist.range; +STC_INLINE float crandom_uniform_f32(crandom_eng32_t* rng, crandom_distrib_f32_t dist) { + return dist.offset + crandom_f32(rng) * dist.range; } typedef struct {uint64_t state[4];} crandom_eng64_t; -typedef struct {double min, range;} crandom_uniform_f64_t; +typedef struct {double offset, range;} crandom_distrib_f64_t; /* 64 bit random number generator engine */ STC_API crandom_eng64_t crandom_eng64_with_seq(uint64_t seed, uint64_t seq); @@ -82,18 +83,23 @@ STC_INLINE crandom_eng64_t crandom_eng64_init(uint64_t seed) { /* int random number generator, range [0, 2^64) */ STC_API uint64_t crandom_i64(crandom_eng64_t* rng); +/* double random number in range [low, high). 52 bit resolution. */ STC_INLINE double crandom_f64(crandom_eng64_t* rng) { union {uint64_t i; double f;} u = {0x3FF0000000000000ull | (crandom_i64(rng) >> 12)}; return u.f - 1.0; } -/* double random number in range [low, high). 52 bit resolution. */ -STC_INLINE crandom_uniform_f64_t crandom_uniform_f64_init(float low, float high) { - crandom_uniform_f64_t dist = {low, high - low}; return dist; +STC_INLINE crandom_distrib_f64_t crandom_uniform_f64_init(double low, double high) { + crandom_distrib_f64_t dist = {low, high - low}; return dist; +} +STC_INLINE double crandom_uniform_f64(crandom_eng64_t* rng, crandom_distrib_f64_t dist) { + return dist.offset + crandom_f64(rng) * dist.range; } -STC_INLINE double crandom_uniform_f64(crandom_eng64_t* rng, crandom_uniform_f64_t dist) { - return dist.min + crandom_f64(rng) * dist.range; + +STC_INLINE crandom_distrib_f64_t crandom_normal_f64_init(double mean, double std_dev) { + crandom_distrib_f64_t dist = {mean, std_dev}; return dist; } +STC_API double crandom_normal_f64(crandom_eng64_t* rng, crandom_distrib_f64_t dist); #if !defined(STC_HEADER) || defined(STC_IMPLEMENTATION) @@ -122,22 +128,10 @@ STC_API crandom_eng64_t crandom_eng64_with_seq(uint64_t seed, uint64_t seq) { return rng; } -#ifdef STC_USE_SFC64 -/* SFC64 random number generator: http://pracrand.sourceforge.net */ -STC_API uint64_t crandom_i64(crandom_eng64_t* rng) { - enum {LROT = 24, RSHIFT = 11, LSHIFT = 3}; - uint64_t *s = rng->state; - const uint64_t result = s[0] + s[1] + s[3]++; - s[0] = s[1] ^ (s[1] >> RSHIFT); - s[1] = s[2] + (s[2] << LSHIFT); - s[2] = ((s[2] << LROT) | (s[2] >> (64 - LROT))) + result; - return result; -} -#else -/* My own PRNG inspired by SFC64: Faster and has Weyl-sequence parameter. */ -/* Copyright Tyge Løvset, 2020 */ -/* Faster: updates only 192bit state. Parallel: Ensures unique sequence per seq (2^63) */ -/* Minimum period is 2^64 per seq, average ~ 2^127 per seq */ +/* PRNG copyright Tyge Løvset, NORCE Research, 2020 */ +/* Extremely fast PRNG suited for parallel usage with Weyl-sequence parameter. */ +/* Updates only 192bit state. Parallel: Ensures unique sequence per seq (2^63) */ +/* Minimum period is 2^64 per seq, but high average per Weyl seq. */ STC_API uint64_t crandom_i64(crandom_eng64_t* rng) { enum {LROT = 24, RSHIFT = 11, LSHIFT = 3}; uint64_t *s = rng->state; @@ -146,7 +140,23 @@ STC_API uint64_t crandom_i64(crandom_eng64_t* rng) { s[1] = ((b << LROT) | (b >> (64 - LROT))) + result; return result; } -#endif + +STC_API double crandom_normal_f64(crandom_eng64_t* rng, crandom_distrib_f64_t dist) { + static bool spare = false; /* Marsaglia polar method: */ + static double u2; double u1, s, m; + if (spare) { + spare = false; + return u2 * dist.range + dist.offset; + } + do { + u1 = 2.0 * crandom_f64(rng) - 1.0; + u2 = 2.0 * crandom_f64(rng) - 1.0; + s = u1*u1 + u2*u2; + } while (s >= 1.0 || s == 0.0); + m = sqrt(-2.0 * log(s) / s); + u2 *= m; spare = true; + return (u1 * m) * dist.range + dist.offset; +} #endif #endif -- cgit v1.2.3