From a51861d02f9d8ed8c69fd243de18960307b0fd63 Mon Sep 17 00:00:00 2001 From: Tyge Løvset Date: Mon, 14 Dec 2020 23:00:15 +0100 Subject: Reworked the crandom.h module. Changed API (sorry!), made uniform distr. unbiased. Removed 32-bit pcg32 PRNG. --- docs/cpqueue_api.md | 6 +- docs/crandom_api.md | 78 ++++++++------------ examples/benchmark.cpp | 6 +- examples/birthday.c | 10 +-- examples/ex_gaussian.c | 6 +- examples/heap.c | 9 +-- examples/list.c | 6 +- examples/priority.c | 8 +- examples/queue.c | 10 +-- examples/random.c | 56 ++++---------- examples/rngtest.c | 44 ++++------- stc/cbitset.h | 25 ++++--- stc/clist.h | 4 +- stc/cpqueue.h | 6 +- stc/cqueue.h | 8 +- stc/crandom.h | 194 +++++++++++++++++-------------------------------- 16 files changed, 182 insertions(+), 294 deletions(-) diff --git a/docs/cpqueue_api.md b/docs/cpqueue_api.md index 39569992..1e92277a 100644 --- a/docs/cpqueue_api.md +++ b/docs/cpqueue_api.md @@ -63,13 +63,13 @@ using_cpqueue(i, cvec_i, >); // adaptor type, '>' = min-heap int main() { size_t N = 10000000; - crand_rng64_t rng = crand_rng64_init(1234); - crand_uniform_i64_t dist = crand_uniform_i64_init(0, N * 10); + cstc64_t rng = cstc64_init(1234); + cstc64_uniform_t dist = cstc64_uniform_init(0, N * 10); cpqueue_i heap = cpqueue_i_init(); // Push ten million random numbers to priority queue, plus some negative ones. c_forrange (N) - cpqueue_i_push(&heap, crand_uniform_i64(&rng, &dist)); + cpqueue_i_push(&heap, cstc64_uniform(&rng, &dist)); c_push_items(&heap, cpqueue_i, {-231, -32, -873, -4, -343}); // Extract and display the fifty smallest. diff --git a/docs/crandom_api.md b/docs/crandom_api.md index 9af69b9e..57697067 100644 --- a/docs/crandom_api.md +++ b/docs/crandom_api.md @@ -1,19 +1,16 @@ -# Module crand: Pseudo Random Number Generators +# Module crandom: Pseudo Random Number Generators This describes the API of module **crand**. It contains *pcg32* created by Melissa O'Neill, and an *64-bit PRNG* created by Tyge Løvset. The PRNG can generate uniform and normal distributions. ## Types -| Name | Type definition | Used to represent... | -|:----------------------|:--------------------------------------------|:-------------------------------------| -| `crand_rng32_t` | `struct {uint64_t state[2];}` | The crandom type | -| `crand_uniform_i32_t` | `struct {int32_t offset; uint32_t range;}` | The crandom element type | -| `crand_uniform_f32_t` | `struct {float offset, range;}` | crandom iterator | -| `crand_rng64_t` | `struct {uint64_t state[4];}` | | -| `crand_uniform_i64_t` | `struct {int64_t offset; uint64_t range;}` | | -| `crand_uniform_f64_t` | `struct {double offset, range;}` | | -| `crand_normal_f64_t` | `struct {double mean, stddev, ...;}` | | +| Name | Type definition | Used to represent... | +|:--------------------|:--------------------------------------------|:-----------------------------| +| `cstc64_t` | `struct {uint64_t state[4];}` | The PRNG engine type | +| `cstc64_uniform_t` | `struct {int64_t offset; uint64_t range;}` | Integer uniform distribution | +| `cstc64_uniformf_t` | `struct {double offset, range;}` | Real number uniform distr. | +| `cstc64_normalf_t` | `struct {double mean, stddev;}` | Normal distribution type | ## Header file @@ -25,41 +22,28 @@ All cstr definitions and prototypes may be included in your C source file by inc ## Methods ```c - 1) crand_rng32_t crand_rng32_init(uint64_t seed); - 2) crand_rng32_t crand_rng32_with_seq(uint64_t seed, uint64_t seq); - 3) uint32_t crand_i32(crand_rng32_t* rng); - 4) float crand_f32(crand_rng32_t* rng); - 5) crand_uniform_i32_t crand_uniform_i32_init(int32_t low, int32_t high); - 6) int32_t crand_uniform_i32(crand_rng32_t* rng, crand_uniform_i32_t* dist); - 7) uint32_t crand_unbiased_i32(crand_rng32_t* rng, crand_uniform_i32_t* dist); - 8) crand_uniform_f32_t crand_uniform_f32_init(float low, float high); /* */ - 9) float crand_uniform_f32(crand_rng32_t* rng, crand_uniform_f32_t* dist); + 1) cstc64_t cstc64_with_seq(uint64_t seed, uint64_t seq); + 2) cstc64_t cstc64_init(uint64_t seed); + 3) uint64_t cstc64_rand(cstc64_t* rng); + 4) double cstc64_randf(cstc64_t* rng); + 5) cstc64_uniform_t cstc64_uniform_init(int64_t low, int64_t high); + 6) int64_t cstc64_uniform(cstc64_t* rng, cstc64_uniform_t* dist); + 7) cstc64_uniformf_t cstc64_uniformf_init(double low, double high); + 8) double cstc64_uniformf(cstc64_t* rng, cstc64_uniformf_t* dist); + 9) cstc64_normalf_t cstc64_normalf_init(double mean, double stddev); +10) double cstc64_normalf(cstc64_t* rng, cstc64_normalf_t* dist); ``` -`1-2)` PRNG 32-bit engine initializers. `3)` Integer RNG with range \[0, 2^32). `4)` Float RNG with range \[0, 1). -`5-6)` Uniform integer RNG with range \[*low*, *high*]. `7)` Unbiased version, see https://github.com/lemire/fastrange. -`8-9)` Uniform float RNG with range \[*low*, *high*). -```c - 1) crand_rng64_t crand_rng64_with_seq(uint64_t seed, uint64_t seq); - 2) crand_rng64_t crand_rng64_init(uint64_t seed); - 3) uint64_t crand_i64(crand_rng64_t* rng); - 4) double crand_f64(crand_rng64_t* rng); - 5) crand_uniform_i64_t crand_uniform_i64_init(int64_t low, int64_t high); - 6) int64_t crand_uniform_i64(crand_rng64_t* rng, crand_uniform_i64_t* dist); - 7) crand_uniform_f64_t crand_uniform_f64_init(double low, double high); - 8) double crand_uniform_f64(crand_rng64_t* rng, crand_uniform_f64_t* dist); - 9) crand_normal_f64_t crand_normal_f64_init(double mean, double stddev); -10) double crand_normal_f64(crand_rng64_t* rng, crand_normal_f64_t* dist); -``` -`1-2)` PRNG 64-bit engine initializers. `3)` Integer generator, range \[0, 2^64). `4)` Double RNG with range \[0, 1). -`5-6)` Uniform integer RNG with range \[*low*, *high*]. `7-8)` Uniform double RNG with range \[*low*, *high*). -`9-10)` Normal-distributed double RNG were around 68% of the values are within the range [*mean* - *stddev, *mean* + *stddev*]. - -The method `crand_i64(crand_rng64_t* rng)` is an extremely fast PRNG suited for parallel usage, featuring -a Weyl-sequence as part of the state. It is faster than *sfc64*, *wyhash64*, *pcg*, and the *xoroshiro* -families of RPNGs. It does not require fast multiplication or 128-bit integer operations. The state is -256-bits, but updates only 192 bit per generated number. It can create 2^63 unique threads with minimum period -length of 2^64 per thread (expected period lengths 2^127). There is no *jump function*, but incrementing -the Weyl-increment by 2, achieves the same goal. Passes *PractRand*, tested up to 8TB output. +`1-2)` PRNG 64-bit engine initializers. `3)` Integer generator, range \[0, 2^64). +`4)` Double RNG with range \[0, 1). `5-6)` Uniform integer RNG with range \[*low*, *high*]. +`7-8)` Uniform double RNG with range \[*low*, *high*). `9-10)` Normal-distributed double +RNG, around 68% of the values are within the range [*mean* - *stddev, *mean* + *stddev*]. + +The method `cstc64_rand(cstc64_t* rng)` is an extremely fast PRNG by Tyge Løvset, suited for parallel usage. +It features a Weyl-sequence as part of the state. It is faster than *sfc64*, *wyhash64*, *pcg64*, and almost +50% faster than *xoshiro256\*\** on GCC. It does not require fast multiplication or 128-bit integer operations. +It has a 256 bit state, but updates only 192 bit per generated number. There is no *jump function*, but incrementing +the Weyl-increment by 2, initializes a new unique 2^64 *minimum* length period. Passes *PractRand*, tested up to +8TB output, Vigna's Hamming weight test, and simple correlation tests, i.e. interleaved streams with one-bit diff state. ## Example ```c @@ -91,13 +75,13 @@ 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(Mean, StdDev); + cstc64_t rng = cstc64_init(seed); + cstc64_normalf_t dist = cstc64_normalf_init(Mean, StdDev); // Create histogram map cmap_i mhist = cmap_i_init(); for (size_t i = 0; i < N; ++i) { - int index = (int) round( crand_normal_f64(&rng, &dist) ); + int index = (int) round( cstc64_normalf(&rng, &dist) ); cmap_i_emplace(&mhist, index, 0).first->second += 1; } diff --git a/examples/benchmark.cpp b/examples/benchmark.cpp index ca744e4d..4830c9ba 100644 --- a/examples/benchmark.cpp +++ b/examples/benchmark.cpp @@ -28,9 +28,9 @@ KHASH_MAP_INIT_INT64(ii, int64_t) size_t seed; static const float max_load_factor = 0.77f; -crand_rng64_t rng; -#define SEED(s) rng = crand_rng64_init(seed) -#define RAND(N) (crand_i64(&rng) & ((1 << N) - 1)) +cstc64_t rng; +#define SEED(s) rng = cstc64_init(seed) +#define RAND(N) (cstc64_rand(&rng) & ((1 << N) - 1)) #define CMAP_SETUP(X, Key, Value) cmap_##X map = cmap_inits \ diff --git a/examples/birthday.c b/examples/birthday.c index 54ffa30e..fedb2c5c 100644 --- a/examples/birthday.c +++ b/examples/birthday.c @@ -14,12 +14,12 @@ const static uint64_t mask = (1ull << 52) - 1; void repeats(void) { - crand_rng64_t rng = crand_rng64_init(seed); + cstc64_t rng = cstc64_init(seed); cmap_ic m = cmap_ic_init(); cmap_ic_reserve(&m, N); clock_t now = clock(); c_forrange (i, N) { - uint64_t k = crand_i64(&rng) & mask; + uint64_t k = cstc64_rand(&rng) & mask; int v = ++cmap_ic_emplace(&m, k, 0).first->second; if (v > 1) printf("%zu: %llx - %d\n", i, k, v); } @@ -33,14 +33,14 @@ using_cvec(x, uint64_t); void distribution(void) { - crand_rng32_t rng = crand_rng32_init(seed); // time(NULL), time(NULL)); + cstc64_t rng = cstc64_init(seed); // time(NULL), time(NULL)); 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(0, M); + cstc64_uniform_t dist = cstc64_uniform_init(0, M); c_forrange (N) { - ++cmap_x_emplace(&map, crand_uniform_i32(&rng, &dist), 0).first->second; + ++cmap_x_emplace(&map, cstc64_uniform(&rng, &dist), 0).first->second; } float diff = (float) (clock() - now) / CLOCKS_PER_SEC; diff --git a/examples/ex_gaussian.c b/examples/ex_gaussian.c index 3ed30cf2..5303dd0c 100644 --- a/examples/ex_gaussian.c +++ b/examples/ex_gaussian.c @@ -25,13 +25,13 @@ 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(Mean, StdDev); + cstc64_t rng = cstc64_init(seed); + cstc64_normalf_t dist = cstc64_normalf_init(Mean, StdDev); // Create histogram map cmap_i mhist = cmap_i_init(); for (size_t i = 0; i < N; ++i) { - int index = (int) round( crand_normal_f64(&rng, &dist) ); + int index = (int) round( cstc64_normalf(&rng, &dist) ); cmap_i_emplace(&mhist, index, 0).first->second += 1; } diff --git a/examples/heap.c b/examples/heap.c index 873077fe..166dd6ba 100644 --- a/examples/heap.c +++ b/examples/heap.c @@ -10,15 +10,15 @@ using_cpqueue(f, cvec_f, >); int main() { uint32_t seed = time(NULL); - crand_rng32_t pcg; + cstc64_t rng; int N = 3000000, M = 100; cpqueue_f pq = cpqueue_f_init(); - pcg = crand_rng32_init(seed); + rng = cstc64_init(seed); clock_t start = clock(); c_forrange (i, int, N) - cvec_f_push_back(&pq, (float) crand_f32(&pcg)*100000); + cvec_f_push_back(&pq, (float) cstc64_randf(&rng)*100000); cpqueue_f_make_heap(&pq); printf("Built priority queue: %f secs\n", (clock() - start) / (float) CLOCKS_PER_SEC); @@ -33,10 +33,9 @@ int main() cpqueue_f_pop(&pq); printf("\n\npopped PQ: %f secs\n", (clock() - start) / (float) CLOCKS_PER_SEC); - pcg = crand_rng32_init(seed); start = clock(); c_forrange (i, int, N) - cpqueue_f_push(&pq, (float) crand_f32(&pcg)*100000); + cpqueue_f_push(&pq, (float) cstc64_randf(&rng)*100000); printf("pushed PQ: %f secs\n", (clock() - start) / (float) CLOCKS_PER_SEC); c_forrange (i, int, M) { diff --git a/examples/list.c b/examples/list.c index 4758b072..fbfa8830 100644 --- a/examples/list.c +++ b/examples/list.c @@ -9,11 +9,11 @@ int main() { const int n = 2000000; clist_fx list = clist_inits; - crand_rng64_t eng = crand_rng64_init(1234); - crand_uniform_f64_t dist = crand_uniform_f64_init(100.0f, n); + cstc64_t eng = cstc64_init(1234); + cstc64_uniformf_t dist = cstc64_uniformf_init(100.0f, n); int m = 0; c_forrange (i, int, n) - clist_fx_push_back(&list, crand_uniform_f64(&eng, &dist)), ++m; + clist_fx_push_back(&list, cstc64_uniformf(&eng, &dist)), ++m; double sum = 0.0; printf("sumarize %d:\n", m); c_foreach (i, clist_fx, list) diff --git a/examples/priority.c b/examples/priority.c index e4a20da9..e81cd8a3 100644 --- a/examples/priority.c +++ b/examples/priority.c @@ -11,19 +11,19 @@ using_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(0, N * 10); + cstc64_t pcg = cstc64_init(time(NULL)); + cstc64_uniform_t dist = cstc64_uniform_init(0, N * 10); cpqueue_i heap = cpqueue_i_init(); // Push ten million random numbers to priority queue c_forrange (N) - cpqueue_i_push(&heap, crand_uniform_i64(&pcg, &dist)); + cpqueue_i_push(&heap, cstc64_uniform(&pcg, &dist)); // push some negative numbers too. c_push_items(&heap, cpqueue_i, {-231, -32, -873, -4, -343}); c_forrange (N) - cpqueue_i_push(&heap, crand_uniform_i64(&pcg, &dist)); + cpqueue_i_push(&heap, cstc64_uniform(&pcg, &dist)); // Extract the hundred smallest. diff --git a/examples/queue.c b/examples/queue.c index c56214a7..0ff29736 100644 --- a/examples/queue.c +++ b/examples/queue.c @@ -7,20 +7,20 @@ using_cqueue(i, clist_i); // min-heap (increasing values) int main() { int n = 10000000; - crand_uniform_i32_t dist; - crand_rng32_t rng = crand_rng32_init(1234); - dist = crand_uniform_i32_init(0, n); + cstc64_uniform_t dist; + cstc64_t rng = cstc64_init(1234); + dist = cstc64_uniform_init(0, n); cqueue_i queue = cqueue_i_init(); // Push ten million random numbers onto the queue. c_forrange (n) - cqueue_i_push(&queue, crand_uniform_i32(&rng, &dist)); + cqueue_i_push(&queue, cstc64_uniform(&rng, &dist)); // Push or pop on the queue ten million times printf("%d\n", n); c_forrange (n) { // range uses initial n only. - int r = crand_uniform_i32(&rng, &dist); + int r = cstc64_uniform(&rng, &dist); if (r & 1) ++n, cqueue_i_push(&queue, r); else diff --git a/examples/random.c b/examples/random.c index 7f1be3f5..4834a974 100644 --- a/examples/random.c +++ b/examples/random.c @@ -6,62 +6,38 @@ int main() { enum {R = 30}; - const size_t N = 1000000000; - clock_t difference, before; + const size_t N = 100000000; + clock_t diff, before; uint64_t sum = 0; - 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(0, range); - - printf("32 uniform: %u\n", dist0.range); - double fsum = 0; - before = clock(); - c_forrange (N) { - 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(0, range); - puts("32 unbiased"); - fsum = 0; - before = clock(); - c_forrange (N) { - 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(0, N); + cstc64_t stc = cstc64_init(seed); + cstc64_uniform_t dist1 = cstc64_uniform_init(0, N); sum = 0; before = clock(); c_forrange (N) { - sum += crand_uniform_i64(&stc, &dist1); + sum += cstc64_uniform(&stc, &dist1); } - difference = clock() - before; - printf("%zu %f: %f secs\n", N, (double) sum / N, (float) difference / CLOCKS_PER_SEC); - + diff = clock() - before; + printf("uniform: %zu %f: %f secs\n", N, (double) sum / N, (float) diff / CLOCKS_PER_SEC); - puts("normal distribution"); - crand_normal_f64_t dist2 = crand_normal_f64_init(R / 2.0, R / 6.0); + cstc64_normalf_t dist2 = cstc64_normalf_init(R / 2.0, R / 6.0); size_t N2 = 10000000; int hist[R] = {0}; sum = 0; + before = clock(); c_forrange (N2) { - int n = (int) (crand_normal_f64(&stc, &dist2) + 0.5); + int n = (int) (cstc64_normalf(&stc, &dist2) + 0.5); sum += n; if (n >= 0 && n < R) ++hist[n]; } - - cstr_t bar = cstr_init(); + diff = clock() - before; + printf("normal : %zu %f: %f secs\n", N, (double) sum / N2, (float) diff / CLOCKS_PER_SEC); + + cstr_t bar = cstr_inits; c_forrange (i, int, R) { - cstr_take(&bar, cstr_with_size(hist[i] * 25ull * R / N2, '*')); - printf("%2d %s\n", i, bar.str); + cstr_resize(&bar, hist[i] * 25ull * R / N2, '*'); + printf("%3d %s\n", i, bar.str); } cstr_del(&bar); } \ No newline at end of file diff --git a/examples/rngtest.c b/examples/rngtest.c index df3aa41c..76666484 100644 --- a/examples/rngtest.c +++ b/examples/rngtest.c @@ -6,49 +6,35 @@ #endif -#define NN 3000000000 +#define NN 1000000000 int main(void) { - clock_t difference, before; + clock_t diff, before; uint64_t v; - crand_rng64_t stc = crand_rng64_init(time(NULL)); - crand_uniform_i64_t idist = crand_uniform_i64_init(10, 20); - crand_uniform_f64_t fdist = crand_uniform_f64_init(10, 20); + cstc64_t stc = cstc64_init(time(NULL)); + cstc64_uniform_t idist = cstc64_uniform_init(10, 20); + cstc64_uniformf_t fdist = cstc64_uniformf_init(10, 20); - c_forrange (30) 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(10, 20); - crand_uniform_f32_t f32dist = crand_uniform_f32_init(10, 20); - - before = clock(); \ + before = clock(); v = 0; c_forrange (NN) { - //v += crand_i32(&pcg); - v += crand_uniform_i32(&pcg, &i32dist); + v += cstc64_rand(&stc); } - difference = clock() - before; - printf("pcg32: %.02f, %zu\n", (float) difference / CLOCKS_PER_SEC, v); + diff = clock() - before; + printf("stc64_rand: %.02f, %zu\n", (float) diff / CLOCKS_PER_SEC, v); - before = clock(); \ + before = clock(); v = 0; c_forrange (NN) { - //v += crand_i64(&stc) & 0xffffffff; - v += crand_uniform_i64(&stc, &idist); + v += cstc64_uniform(&stc, &idist); } - difference = clock() - before; - printf("stc64: %.02f, %zu\n", (float) difference / CLOCKS_PER_SEC, v); + diff = clock() - before; + printf("stc64_uniform: %.02f, %zu\n\n", (float) diff / CLOCKS_PER_SEC, v); - c_forrange (8) printf("%d ", crand_uniform_i32(&pcg, &i32dist)); + c_forrange (30) printf("%02zd ", cstc64_uniform(&stc, &idist)); puts(""); - - - c_forrange (8) printf("%f ", crand_uniform_f32(&pcg, &f32dist)); - puts(""); - - c_forrange (8) printf("%f ", crand_uniform_f64(&stc, &fdist)); + c_forrange (8) printf("%f ", cstc64_uniformf(&stc, &fdist)); puts(""); } \ No newline at end of file diff --git a/stc/cbitset.h b/stc/cbitset.h index 715eddb4..2e1887a2 100644 --- a/stc/cbitset.h +++ b/stc/cbitset.h @@ -51,20 +51,21 @@ int main() { typedef struct cbitset { uint64_t* _arr; size_t size; } cbitset_t; -STC_API cbitset_t cbitset_with_size(size_t size, bool value); -STC_API cbitset_t cbitset_from_str(const char* str); -STC_API char* cbitset_to_str(cbitset_t set, char* str, size_t start, intptr_t stop); -STC_API cbitset_t cbitset_clone(cbitset_t other); -STC_API void cbitset_resize(cbitset_t* self, size_t size, bool value); -STC_API size_t cbitset_count(cbitset_t set); -STC_API bool cbitset_is_disjoint(cbitset_t set, cbitset_t other); -STC_API bool cbitset_is_subset(cbitset_t set, cbitset_t other); -STC_API bool cbitset_is_superset(cbitset_t set, cbitset_t other); - -STC_INLINE size_t cbitset_size(cbitset_t set) {return set.size;} +STC_API cbitset_t cbitset_with_size(size_t size, bool value); +STC_API cbitset_t cbitset_from_str(const char* str); +STC_API char* cbitset_to_str(cbitset_t set, char* str, size_t start, intptr_t stop); +STC_API cbitset_t cbitset_clone(cbitset_t other); +STC_API void cbitset_resize(cbitset_t* self, size_t size, bool value); +STC_API size_t cbitset_count(cbitset_t set); +STC_API bool cbitset_is_disjoint(cbitset_t set, cbitset_t other); +STC_API bool cbitset_is_subset(cbitset_t set, cbitset_t other); +STC_API bool cbitset_is_superset(cbitset_t set, cbitset_t other); + +STC_INLINE size_t cbitset_size(cbitset_t set) {return set.size;} +#define cbitset_inits {NULL, 0} STC_INLINE cbitset_t cbitset_init() { - cbitset_t bs = {NULL, 0}; return bs; + cbitset_t bs = cbitset_inits; return bs; } STC_INLINE void cbitset_del(cbitset_t* self) { c_free(self->_arr); diff --git a/stc/clist.h b/stc/clist.h index 768de474..f322fd00 100644 --- a/stc/clist.h +++ b/stc/clist.h @@ -36,10 +36,10 @@ int main() { clist_ix list = clist_inits; - crand_rng32_t pcg = crand_rng32_init(12345); + cstc64_t rng = cstc64_init(12345); int n; for (int i=0; i<1000000; ++i) // one million - clist_ix_push_back(&list, crand_i32(&pcg)); + clist_ix_push_back(&list, cstc64_rand(&rng) >> 32); n = 0; c_foreach (i, clist_ix, list) if (++n % 10000 == 0) printf("%8d: %10zd\n", n, i.val->value); diff --git a/stc/cpqueue.h b/stc/cpqueue.h index f7f8150c..2bdf2807 100644 --- a/stc/cpqueue.h +++ b/stc/cpqueue.h @@ -29,13 +29,13 @@ using_cpqueue(f, cvec_f, >); // min-heap (increasing values) int main() { - crand_rng32_t gen = crand_rng32_init(1234); - crand_uniform_f32_t dist = crand_uniform_f32_init(10.0f, 100.0f); + cstc64_t gen = cstc64_init(1234); + cstc64_uniformf_t dist = cstc64_uniformf_init(10.0f, 100.0f); cpqueue_f queue = cpqueue_f_init(); // Push ten million random numbers onto the queue. for (int i=0; i<10000000; ++i) - cpqueue_f_push(&queue, crand_uniform_f32(&gen, dist)); + cpqueue_f_push(&queue, cstc64_uniformf(&gen, dist)); // Extract the 100 smallest. for (int i=0; i<100; ++i) { printf("%f ", *cpqueue_f_top(queue)); diff --git a/stc/cqueue.h b/stc/cqueue.h index 7ad87c62..6226b59b 100644 --- a/stc/cqueue.h +++ b/stc/cqueue.h @@ -32,19 +32,19 @@ int main() { int n = 10000000; - crand_rng32_t gen = crand_rng32_init(1234); - crand_uniform_i32_t dist = crand_uniform_i32_init(gen, 0, n); + cstc64_t gen = cstc64_init(1234); + cstc64_uniform_t dist = cstc64_uniform_init(gen, 0, n); cqueue_i queue = cqueue_i_init(); // Push ten million random numbers onto the queue. for (int i=0; i0; --i) { - int r = crand_uniform_i32(&dist); + int r = cstc64_uniform(&dist); if (r & 1) ++n, cqueue_i_push(&queue, r); else diff --git a/stc/crandom.h b/stc/crandom.h index edd31113..fa353345 100644 --- a/stc/crandom.h +++ b/stc/crandom.h @@ -24,176 +24,118 @@ #define CRANDOM__H__ /* - crand_rng32_t rng = crand_rng32_init(seed); - crand_uniform_f32_t fdist = crand_uniform_f32_init(rng, 1.0f, 6.0f); - crand_uniform_i32_t idist = crand_uniform_i32_init(rng, 1, 6); +// cstc64: Pseudo-random number generator +#include "stc/crandom.h" +int main() { + uint64_t seed = 123456789; + cstc64_t rng = cstc64_init(seed); + cstc64_uniform_t dist1 = cstc64_uniform_init(1, 6); + cstc64_uniformf_t dist2 = cstc64_uniformf_init(1.0, 10.0); + cstc64_normalf_t dist3 = cstc64_normalf_init(1.0, 10.0); - uint32_t i = crand_i32(&rng); - int j = crand_uniform_i32(&idist); - float r = crand_uniform_f32(&fdist); + uint64_t i = cstc64_rand(&rng); + int64_t iu = cstc64_uniform(&rng, &dist1); + double xu = cstc64_uniformf(&rng, &dist2); + double xn = cstc64_normalf(&rng, &dist3); +} */ #include "ccommon.h" #include #include -/* 32-BIT RANDOM NUMBER GENERATOR */ - -typedef struct {uint64_t state[2];} crand_rng32_t; -typedef struct {int32_t offset; uint32_t range;} crand_uniform_i32_t; -typedef struct {float offset, range;} crand_uniform_f32_t; - -/* engine initializers */ -STC_API crand_rng32_t crand_rng32_with_seq(uint64_t seed, uint64_t seq); -STC_INLINE crand_rng32_t crand_rng32_init(uint64_t seed) { - return crand_rng32_with_seq(seed, seed); -} - -/* int random number generator, range [0, 2^32) */ -STC_API uint32_t crand_i32(crand_rng32_t* rng); - -STC_INLINE float crand_f32(crand_rng32_t* rng) { - union {uint32_t i; float f;} u = {0x3F800000u | (crand_i32(rng) >> 9)}; - return u.f - 1.0f; -} - -/* int random number generator in range [low, high] */ -STC_INLINE crand_uniform_i32_t crand_uniform_i32_init(int32_t low, int32_t high) { - crand_uniform_i32_t dist = {low, (uint32_t) (high - low + 1)}; return dist; -} -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_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(float low, float high) { - crand_uniform_f32_t dist = {low, high - low}; return dist; -} -STC_INLINE float crand_uniform_f32(crand_rng32_t* rng, crand_uniform_f32_t* dist) { - return dist->offset + crand_f32(rng) * dist->range; -} - +typedef struct {uint64_t state[4];} cstc64_t; +typedef struct {int64_t offset; uint64_t range, threshold;} cstc64_uniform_t; +typedef struct {double offset, range;} cstc64_uniformf_t; +typedef struct {double mean, stddev, next; bool has_next;} cstc64_normalf_t; -/* 64 BIT RANDOM NUMBER GENERATOR */ -typedef struct {uint64_t state[4];} crand_rng64_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 */ -STC_API crand_rng64_t crand_rng64_with_seq(uint64_t seed, uint64_t seq); -STC_INLINE crand_rng64_t crand_rng64_init(uint64_t seed) { - return crand_rng64_with_seq(seed, 1); -} /* int random number generator, range [0, 2^64). PRNG copyright Tyge Løvset, NORCE Research, 2020 */ -STC_API uint64_t crand_i64(crand_rng64_t* rng); +STC_API cstc64_t cstc64_init(uint64_t seed); +STC_API cstc64_t cstc64_with_seq(uint64_t seed, uint64_t seq); +STC_API uint64_t cstc64_rand(cstc64_t* rng); -/* double random number in range [low, high). 52 bit resolution. */ -STC_INLINE double crand_f64(crand_rng64_t* rng) { - union {uint64_t i; double f;} u = {0x3FF0000000000000ull | (crand_i64(rng) >> 12)}; +/* double random number in range [low, high). */ +STC_INLINE double cstc64_randf(cstc64_t* rng) { + union {uint64_t i; double f;} u = {0x3FF0000000000000ull | (cstc64_rand(rng) >> 12)}; return u.f - 1.0; } -/* int random number generator in range [low, high] */ -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_rng64_t* rng, crand_uniform_i64_t* dist) { - #if defined(__SIZEOF_INT128__) - 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(rng), dist->range, &hi); return dist->offset + hi; - #else - return dist->offset + crand_i64(rng) % dist->range; // slower - #endif -} +/* integer uniform distributed RNG, range [low, high]. */ +STC_API cstc64_uniform_t cstc64_uniform_init(int64_t low, int64_t high); +STC_API int64_t cstc64_uniform(cstc64_t* rng, cstc64_uniform_t* 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; +/* double uniform distributed RNG, range [low, high). */ +STC_INLINE cstc64_uniformf_t cstc64_uniformf_init(double low, double high) { + cstc64_uniformf_t dist = {low, high - low}; return dist; } -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 double cstc64_uniformf(cstc64_t* rng, cstc64_uniformf_t* dist) { + return dist->offset + cstc64_randf(rng) * dist->range; } -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; +/* double normal distributed RNG. */ +STC_INLINE cstc64_normalf_t cstc64_normalf_init(double mean, double stddev) { + cstc64_normalf_t dist = {mean, stddev, 0.0, false}; return dist; } -STC_API double crand_normal_f64(crand_rng64_t* rng, crand_normal_f64_t* dist); +STC_API double cstc64_normalf(cstc64_t* rng, cstc64_normalf_t* dist); #if !defined(STC_HEADER) || defined(STC_IMPLEMENTATION) #if defined(_MSC_VER) #pragma warning(push) - #pragma warning(disable: 4146) + //#pragma warning(disable: 4146) #endif -/* PRNG PCG32 https://www.pcg-random.org/download.html */ -STC_DEF 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); - rng.state[0] += seed; - crand_i32(&rng); - return rng; -} -STC_DEF uint32_t crand_i32(crand_rng32_t* rng) { - uint64_t old = rng->state[0]; - rng->state[0] = old * 6364136223846793005ull + rng->state[1]; - uint32_t xors = (uint32_t) (((old >> 18u) ^ old) >> 27u); - uint32_t rot = old >> 59u; - return (xors >> rot) | (xors << ((-rot) & 31)); -} +/* PRNG stc64: Tyge Løvset, NORCE Research, 2020. + * Extremely fast PRNG suited for parallel usage with Weyl-sequence parameter. + * Faster than sfc64, wyhash64, and almost 50% faster than xoshiro256** on gcc. + * 256bit state, updates only 192bit per rng. + * 2^63 unique threads with a minimum 2^64 period lengths each. + * 2^127 minimum period length for single thread (double loop). + * Passes PractRand tested up to 8TB output, Vigna's Hamming weight test, + * and simple correlation tests, i.e. interleaved streams with one-bit diff state. + */ -/* 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. */ -/* Tested with PractRand to 8 TB output: no issues */ -STC_DEF crand_rng64_t crand_rng64_with_seq(uint64_t seed, uint64_t seq) { - crand_rng64_t rng = {{seed, seed, seed, (seq << 1u) | 1u}}; /* increment must be odd */ - for (int i = 0; i < 12; ++i) crand_i64(&rng); +STC_DEF cstc64_t cstc64_init(uint64_t seed) { + return cstc64_with_seq(seed, 0x3504f333d3aa0b34); +} +STC_DEF cstc64_t cstc64_with_seq(uint64_t seed, uint64_t seq) { + cstc64_t rng = {{seed, seed, seed, (seq << 1u) | 1u}}; + for (int i = 0; i < 8; ++i) cstc64_rand(&rng); return rng; } -STC_DEF uint64_t crand_i64(crand_rng64_t* rng) { +STC_DEF uint64_t cstc64_rand(cstc64_t* rng) { enum {LROT = 24, RSHIFT = 11, LSHIFT = 3}; uint64_t *s = rng->state; - const uint64_t b = s[1], result = s[0] ^ (s[2] += s[3]|1); - s[0] = (b + (b << LSHIFT)) ^ (b >> RSHIFT); - s[1] = ((b << LROT) | (b >> (64 - LROT))) + result; + const uint64_t result = s[0] ^ (s[2] += s[3]|1); + s[0] = (s[1] + (s[1] << LSHIFT)) ^ (s[1] >> RSHIFT); + s[1] = ((s[1] << LROT) | (s[1] >> (64 - LROT))) + result; return result; } -/* Unbiased uniform https://github.com/lemire/fastrange */ -STC_DEF 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(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(rng) * r); +/* unbiased integer random number generator in range [low, high] */ +STC_DEF cstc64_uniform_t cstc64_uniform_init(int64_t low, int64_t high) { + cstc64_uniform_t dist = {low, (uint64_t) (high - low + 1)}; + dist.threshold = (uint64_t)(-dist.range) % dist.range; + return dist; +} +STC_DEF int64_t cstc64_uniform(cstc64_t* rng, cstc64_uniform_t* dist) { + for (;;) { + uint64_t r = cstc64_rand(rng); + if (r >= dist->threshold) return dist->offset + (r % dist->range); } - return dist->offset + (m >> 32); } /* Marsaglia polar method for gaussian distribution. */ -STC_DEF double crand_normal_f64(crand_rng64_t* rng, crand_normal_f64_t* dist) { +STC_DEF double cstc64_normalf(cstc64_t* rng, cstc64_normalf_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(rng) - 1.0; - u2 = 2.0 * crand_f64(rng) - 1.0; + u1 = 2.0 * cstc64_randf(rng) - 1.0; + u2 = 2.0 * cstc64_randf(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