From b3a1581b7dabd0fe8989f605bdcf1fee93f657d9 Mon Sep 17 00:00:00 2001 From: Tyge Løvset Date: Sun, 13 Sep 2020 11:43:18 +0200 Subject: reverted back to separate random engine and distribution parameters. --- examples/birthday.c | 4 +-- examples/ex_gaussian.c | 4 +-- examples/list.c | 4 +-- examples/priority.c | 6 ++--- examples/queue.c | 7 ++--- examples/random.c | 16 ++++++------ examples/rngtest.c | 20 +++++++------- stc/crandom.h | 71 +++++++++++++++++++++++++------------------------- 8 files changed, 66 insertions(+), 66 deletions(-) diff --git a/examples/birthday.c b/examples/birthday.c index bf807e67..8b6abd77 100644 --- a/examples/birthday.c +++ b/examples/birthday.c @@ -38,9 +38,9 @@ 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(); - crand_uniform_i32_t dist = crand_uniform_i32_init(rng, 0, M); + crand_uniform_i32_t dist = crand_uniform_i32_init(0, M); for (size_t i = 0; i < N; ++i) { - ++cmap_x_emplace(&map, crand_uniform_i32(&dist), 0).item->value; + ++cmap_x_emplace(&map, crand_uniform_i32(&rng, &dist), 0).item->value; } float diff = (float) (clock() - now) / CLOCKS_PER_SEC; diff --git a/examples/ex_gaussian.c b/examples/ex_gaussian.c index d2da77d9..cd44e78a 100644 --- a/examples/ex_gaussian.c +++ b/examples/ex_gaussian.c @@ -26,12 +26,12 @@ int main() // Setup random engine with normal distribution. uint64_t seed = time(NULL); crand_rng64_t rng = crand_rng64_init(seed); - crand_normal_f64_t dist = crand_normal_f64_init(rng, Mean, StdDev); + crand_normal_f64_t dist = crand_normal_f64_init(Mean, StdDev); // Create histogram map cmap_i mhist = cmap_ini; for (size_t i = 0; i < N; ++i) { - int index = round( crand_normal_f64(&dist) ); + int index = round( crand_normal_f64(&rng, &dist) ); cmap_i_emplace(&mhist, index, 0).item->value += 1; } diff --git a/examples/list.c b/examples/list.c index 24dfc82f..dd8bba2f 100644 --- a/examples/list.c +++ b/examples/list.c @@ -8,10 +8,10 @@ int main() { int k, n = 100000; clist_fx list = clist_ini; crand_rng64_t eng = crand_rng64_init(time(NULL)); - crand_uniform_f64_t dist = crand_uniform_f64_init(eng, 0.0f, n); + crand_uniform_f64_t dist = crand_uniform_f64_init(0.0f, n); for (int i = 0; i < 100000; ++i) - clist_fx_push_back(&list, crand_uniform_f64(&dist)); + clist_fx_push_back(&list, crand_uniform_f64(&eng, &dist)); k = 0; c_foreach (i, clist_fx, list) if (++k <= 10) printf("%8d: %10f\n", k, i.item->value); else break; diff --git a/examples/priority.c b/examples/priority.c index 82e78621..e857761c 100644 --- a/examples/priority.c +++ b/examples/priority.c @@ -12,18 +12,18 @@ declare_cpqueue(i, cvec_i, >); // min-heap (increasing values) int main() { size_t N = 10000000; crand_rng64_t pcg = crand_rng64_init(time(NULL)); - crand_uniform_i64_t dist = crand_uniform_i64_init(pcg, 0, N * 10); + crand_uniform_i64_t dist = crand_uniform_i64_init(0, N * 10); cpqueue_i heap = cpqueue_i_init(); // Push ten million random numbers to priority queue for (int i=0; i0; --i) { - int r = crand_uniform_i32(&dist); + int r = crand_uniform_i32(&rng, &dist); if (r & 1) ++n, cqueue_i_push(&queue, r); else diff --git a/examples/random.c b/examples/random.c index 16632373..b81f1826 100644 --- a/examples/random.c +++ b/examples/random.c @@ -13,47 +13,47 @@ int main() uint64_t seed = time(NULL); crand_rng32_t pcg = crand_rng32_init(seed); uint32_t range = crand_i32(&pcg) & ((1u << 28) - 1); - crand_uniform_i32_t dist0 = crand_uniform_i32_init(pcg, 0, range); + crand_uniform_i32_t dist0 = crand_uniform_i32_init(0, range); printf("32 uniform: %u\n", dist0.range); double fsum = 0; before = clock(); for (size_t i = 0; i < N; ++i) { - fsum += (double) crand_uniform_i32(&dist0) / dist0.range; + fsum += (double) crand_uniform_i32(&pcg, &dist0) / dist0.range; } difference = clock() - before; printf("%zu %f: %f secs\n", N, fsum / N, (float) difference / CLOCKS_PER_SEC); pcg = crand_rng32_init(seed); - dist0 = crand_uniform_i32_init(pcg, 0, range); + dist0 = crand_uniform_i32_init(0, range); puts("32 unbiased"); fsum = 0; before = clock(); for (size_t i = 0; i < N; ++i) { - fsum += (double) crand_unbiased_i32(&dist0) / dist0.range; + fsum += (double) crand_unbiased_i32(&pcg, &dist0) / dist0.range; } difference = clock() - before; printf("%zu %f: %f secs\n", N, fsum / N, (float) difference / CLOCKS_PER_SEC); puts("64 uniform"); crand_rng64_t stc = crand_rng64_init(seed); - crand_uniform_i64_t dist1 = crand_uniform_i64_init(stc, 0, N); + crand_uniform_i64_t dist1 = crand_uniform_i64_init(0, N); sum = 0; before = clock(); for (size_t i = 0; i < N; ++i) { - sum += crand_uniform_i64(&dist1); + sum += crand_uniform_i64(&stc, &dist1); } difference = clock() - before; printf("%zu %f: %f secs\n", N, (double) sum / N, (float) difference / CLOCKS_PER_SEC); puts("normal distribution"); - crand_normal_f64_t dist2 = crand_normal_f64_init(stc, R / 2.0, R / 6.0); + crand_normal_f64_t dist2 = crand_normal_f64_init(R / 2.0, R / 6.0); size_t N2 = 10000000; int hist[R] = {0}; sum = 0; for (size_t i = 0; i < N2; ++i) { - int n = (int) (crand_normal_f64(&dist2) + 0.5); + int n = (int) (crand_normal_f64(&stc, &dist2) + 0.5); sum += n; if (n >= 0 && n < R) ++hist[n]; } diff --git a/examples/rngtest.c b/examples/rngtest.c index 6aa783e2..dbe272c0 100644 --- a/examples/rngtest.c +++ b/examples/rngtest.c @@ -14,21 +14,21 @@ int main(void) uint64_t v; crand_rng64_t stc = crand_rng64_init(time(NULL)); - crand_uniform_i64_t idist = crand_uniform_i64_init(stc, 10, 20); - crand_uniform_f64_t fdist = crand_uniform_f64_init(stc, 10, 20); + crand_uniform_i64_t idist = crand_uniform_i64_init(10, 20); + crand_uniform_f64_t fdist = crand_uniform_f64_init(10, 20); - for (int i=0; i<30; ++i) printf("%02zd ", crand_uniform_i64(&idist)); + for (int i=0; i<30; ++i) printf("%02zd ", crand_uniform_i64(&stc, &idist)); puts(""); crand_rng32_t pcg = crand_rng32_init(time(NULL)); - crand_uniform_i32_t i32dist = crand_uniform_i32_init(pcg, 10, 20); - crand_uniform_f32_t f32dist = crand_uniform_f32_init(pcg, 10, 20); + crand_uniform_i32_t i32dist = crand_uniform_i32_init(10, 20); + crand_uniform_f32_t f32dist = crand_uniform_f32_init(10, 20); before = clock(); \ v = 0; for (size_t i=0; ioffset + (int32_t) (((uint64_t) crand_i32(&dist->rng) * dist->range) >> 32); +STC_INLINE int32_t crand_uniform_i32(crand_rng32_t* rng, crand_uniform_i32_t* dist) { + return dist->offset + (int32_t) (((uint64_t) crand_i32(rng) * dist->range) >> 32); } /* https://github.com/lemire/fastrange */ -STC_API uint32_t crand_unbiased_i32(crand_uniform_i32_t* dist); +STC_API uint32_t crand_unbiased_i32(crand_rng32_t* rng, crand_uniform_i32_t* dist); /* float random number in range [low, high). Note: 23 bit resolution. */ -STC_INLINE crand_uniform_f32_t crand_uniform_f32_init(crand_rng32_t rng, float low, float high) { - crand_uniform_f32_t dist = {rng, low, high - low}; return dist; +STC_INLINE crand_uniform_f32_t crand_uniform_f32_init(float low, float high) { + crand_uniform_f32_t dist = {low, high - low}; return dist; } -STC_INLINE float crand_uniform_f32(crand_uniform_f32_t* dist) { - return dist->offset + crand_f32(&dist->rng) * dist->range; +STC_INLINE float crand_uniform_f32(crand_rng32_t* rng, crand_uniform_f32_t* dist) { + return dist->offset + crand_f32(rng) * dist->range; } /* 64 BIT RANDOM NUMBER GENERATOR */ typedef struct {uint64_t state[4];} crand_rng64_t; -typedef struct {crand_rng64_t rng; int64_t offset; uint64_t range;} crand_uniform_i64_t; -typedef struct {crand_rng64_t rng; double offset, range;} crand_uniform_f64_t; -typedef struct {crand_rng64_t rng; double mean, stddev, next; bool has_next;} crand_normal_f64_t; +typedef struct {int64_t offset; uint64_t range;} crand_uniform_i64_t; +typedef struct {double offset, range;} crand_uniform_f64_t; +typedef struct {double mean, stddev, next; bool has_next;} crand_normal_f64_t; /* engine initializers */ @@ -99,39 +99,40 @@ STC_INLINE double crand_f64(crand_rng64_t* rng) { } /* int random number generator in range [low, high] */ -STC_INLINE crand_uniform_i64_t crand_uniform_i64_init(crand_rng64_t rng, int64_t low, int64_t high) { - crand_uniform_i64_t dist = {rng, low, (uint64_t) (high - low + 1)}; return dist; +STC_INLINE crand_uniform_i64_t crand_uniform_i64_init(int64_t low, int64_t high) { + crand_uniform_i64_t dist = {low, (uint64_t) (high - low + 1)}; return dist; } #if defined(_MSC_VER) && defined(_WIN64) #include #endif -STC_INLINE int64_t crand_uniform_i64(crand_uniform_i64_t* dist) { +STC_INLINE int64_t crand_uniform_i64(crand_rng64_t* rng, crand_uniform_i64_t* dist) { #if defined(__SIZEOF_INT128__) - return dist->offset + (int64_t) (((__uint128_t) crand_i64(&dist->rng) * dist->range) >> 64); + return dist->offset + (int64_t) (((__uint128_t) crand_i64(rng) * dist->range) >> 64); #elif defined(_MSC_VER) && defined(_WIN64) - uint64_t hi; _umul128(crand_i64(&dist->rng), dist->range, &hi); return dist->offset + hi; + uint64_t hi; _umul128(crand_i64(rng), dist->range, &hi); return dist->offset + hi; #else - return dist->offset + crand_i64(&dist->rng) % dist->range; // slower + return dist->offset + crand_i64(rng) % dist->range; // slower #endif } -STC_INLINE crand_uniform_f64_t crand_uniform_f64_init(crand_rng64_t rng, double low, double high) { - crand_uniform_f64_t dist = {rng, low, high - low}; return dist; +STC_INLINE crand_uniform_f64_t crand_uniform_f64_init(double low, double high) { + crand_uniform_f64_t dist = {low, high - low}; return dist; } -STC_INLINE double crand_uniform_f64(crand_uniform_f64_t* dist) { - return dist->offset + crand_f64(&dist->rng) * dist->range; +STC_INLINE double crand_uniform_f64(crand_rng64_t* rng, crand_uniform_f64_t* dist) { + return dist->offset + crand_f64(rng) * dist->range; } -STC_INLINE crand_normal_f64_t crand_normal_f64_init(crand_rng64_t rng, double mean, double stddev) { - crand_normal_f64_t dist = {rng, mean, stddev, 0.0, false}; return dist; +STC_INLINE crand_normal_f64_t crand_normal_f64_init(double mean, double stddev) { + crand_normal_f64_t dist = {mean, stddev, 0.0, false}; return dist; } -STC_API double crand_normal_f64(crand_normal_f64_t* dist); +STC_API double crand_normal_f64(crand_rng64_t* rng, crand_normal_f64_t* dist); #if !defined(STC_HEADER) || defined(STC_IMPLEMENTATION) +/* PRNG PCG32 https://www.pcg-random.org/download.html */ STC_API crand_rng32_t crand_rng32_with_seq(uint64_t seed, uint64_t seq) { crand_rng32_t rng = {{0u, (seq << 1u) | 1u}}; /* inc must be odd */ crand_i32(&rng); @@ -139,8 +140,6 @@ STC_API crand_rng32_t crand_rng32_with_seq(uint64_t seed, uint64_t seq) { crand_i32(&rng); return rng; } - -/* PCG32 https://www.pcg-random.org/download.html */ STC_API uint32_t crand_i32(crand_rng32_t* rng) { uint64_t old = rng->state[0]; rng->state[0] = old * 6364136223846793005ull + rng->state[1]; @@ -149,7 +148,7 @@ STC_API uint32_t crand_i32(crand_rng32_t* rng) { return (xors >> rot) | (xors << ((-rot) & 31)); } -/* PRNG copyright Tyge Løvset, NORCE Research, 2020 */ +/* PRNG STC64: copyright Tyge Løvset, NORCE Research, 2020 */ /* Extremely fast PRNG suited for parallel usage with Weyl-sequence parameter. */ /* Even faster than sfc64: updates only 192bit state. Better for parallel processing: */ /* Guarantees 2^63 unique threads with minimum 2^64 period length ~ 2^160 average period. */ @@ -169,28 +168,28 @@ STC_API uint64_t crand_i64(crand_rng64_t* rng) { } /* Unbiased uniform https://github.com/lemire/fastrange */ -STC_API uint32_t crand_unbiased_i32(crand_uniform_i32_t* dist) { +STC_API uint32_t crand_unbiased_i32(crand_rng32_t* rng, crand_uniform_i32_t* dist) { uint32_t r = dist->range; - uint64_t m = (uint64_t) crand_i32(&dist->rng) * r; + uint64_t m = (uint64_t) crand_i32(rng) * r; uint32_t l = (uint32_t) m; if (l < r) { uint32_t t = -r; if (t >= r) if ((t -= r) >= r) t %= r; - while (l < t) l = (uint32_t) (m = (uint64_t) crand_i32(&dist->rng) * r); + while (l < t) l = (uint32_t) (m = (uint64_t) crand_i32(rng) * r); } return dist->offset + (m >> 32); } /* Marsaglia polar method for gaussian distribution. */ -STC_API double crand_normal_f64(crand_normal_f64_t* dist) { +STC_API double crand_normal_f64(crand_rng64_t* rng, crand_normal_f64_t* dist) { double u1, u2, s, m; if (dist->has_next) { dist->has_next = false; return dist->next * dist->stddev + dist->mean; } do { - u1 = 2.0 * crand_f64(&dist->rng) - 1.0; - u2 = 2.0 * crand_f64(&dist->rng) - 1.0; + u1 = 2.0 * crand_f64(rng) - 1.0; + u2 = 2.0 * crand_f64(rng) - 1.0; s = u1*u1 + u2*u2; } while (s >= 1.0 || s == 0.0); m = sqrt(-2.0 * log(s) / s); -- cgit v1.2.3