From a147a1be246472d079f72fc22e065080d1588e29 Mon Sep 17 00:00:00 2001 From: Tyge Løvset Date: Wed, 16 Dec 2020 10:58:04 +0100 Subject: Reverted back API to easier naming scheme of crandom.h types and methods. Highly optimized generation of unbiased uniform bounded random numbers. Only 33% overhead. --- docs/cpqueue_api.md | 6 +-- docs/crandom_api.md | 57 ++++++++++++++------------ examples/benchmark.cpp | 6 +-- examples/birthday.c | 10 ++--- examples/ex_gaussian.c | 6 +-- examples/heap.c | 8 ++-- examples/list.c | 6 +-- examples/priority.c | 8 ++-- examples/queue.c | 10 ++--- examples/random.c | 16 ++++---- examples/rngtest.c | 14 +++---- stc/clist.h | 4 +- stc/cpqueue.h | 6 +-- stc/cqueue.h | 8 ++-- stc/crandom.h | 107 +++++++++++++++++++++++++++---------------------- 15 files changed, 144 insertions(+), 128 deletions(-) diff --git a/docs/cpqueue_api.md b/docs/cpqueue_api.md index 1e92277a..ac33e68f 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; - cstc64_t rng = cstc64_init(1234); - cstc64_uniform_t dist = cstc64_uniform_init(0, N * 10); + crand_t rng = crand_init(1234); + crand_uniform_t dist = crand_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, cstc64_uniform(&rng, &dist)); + cpqueue_i_push(&heap, crand_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 57697067..53e8a0d1 100644 --- a/docs/crandom_api.md +++ b/docs/crandom_api.md @@ -5,12 +5,12 @@ This describes the API of module **crand**. It contains *pcg32* created by Melis ## Types -| 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 | +| Name | Type definition | Used to represent... | +|:-------------------|:-------------------------------------------|:-----------------------------| +| `crand_t` | `struct {uint64_t state[4];}` | The PRNG engine type | +| `crand_uniform_t` | `struct {int64_t offset; uint64_t range;}` | Integer uniform distribution | +| `crand_uniformf_t` | `struct {double offset, range;}` | Real number uniform distr. | +| `crand_normalf_t` | `struct {double mean, stddev;}` | Normal distribution type | ## Header file @@ -22,28 +22,33 @@ All cstr definitions and prototypes may be included in your C source file by inc ## Methods ```c - 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) crand_t crand_init(uint64_t seed); + 2) crand_t crand_with_seq(uint64_t seed, uint64_t seq); + 3) uint64_t crand_next(crand_t* rng); + 4) double crand_nextf(crand_t* rng); + 5) crand_uniform_t crand_uniform_init(int64_t low, int64_t high); + 6) int64_t crand_uniform(crand_t* rng, crand_uniform_t* dist); + 7) crand_uniformf_t crand_uniformf_init(double low, double high); + 8) double crand_uniformf(crand_t* rng, crand_uniformf_t* dist); + 9) crand_normalf_t crand_normalf_init(double mean, double stddev); +10) double crand_normalf(crand_t* rng, crand_normalf_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, around 68% of the values are within the range [*mean* - *stddev, *mean* + *stddev*]. +RNG, around 68% of the values fall 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. +**crand** 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\*\** compiled with 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 by incrementing the Weyl-increment by 2, it initializes +a new unique 2^64 *minimum* length period. Note that for each Weyl-increment, the expected period +length is about 2^126. For a single thread, a minimum period of 2^127 is generated if the +Weyl-increment is incremented by 2 every 2^64 number output. + +*crand* 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 @@ -75,13 +80,13 @@ int main() // Setup random engine with normal distribution. uint64_t seed = time(NULL); - cstc64_t rng = cstc64_init(seed); - cstc64_normalf_t dist = cstc64_normalf_init(Mean, StdDev); + crand_t rng = crand_init(seed); + crand_normalf_t dist = crand_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( cstc64_normalf(&rng, &dist) ); + int index = (int) round( crand_normalf(&rng, &dist) ); cmap_i_emplace(&mhist, index, 0).first->second += 1; } diff --git a/examples/benchmark.cpp b/examples/benchmark.cpp index 4830c9ba..206da5ce 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; -cstc64_t rng; -#define SEED(s) rng = cstc64_init(seed) -#define RAND(N) (cstc64_rand(&rng) & ((1 << N) - 1)) +crand_t rng; +#define SEED(s) rng = crand_init(seed) +#define RAND(N) (crand_next(&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 fedb2c5c..ac05a5f5 100644 --- a/examples/birthday.c +++ b/examples/birthday.c @@ -14,12 +14,12 @@ const static uint64_t mask = (1ull << 52) - 1; void repeats(void) { - cstc64_t rng = cstc64_init(seed); + crand_t rng = crand_init(seed); cmap_ic m = cmap_ic_init(); cmap_ic_reserve(&m, N); clock_t now = clock(); c_forrange (i, N) { - uint64_t k = cstc64_rand(&rng) & mask; + uint64_t k = crand_next(&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) { - cstc64_t rng = cstc64_init(seed); // time(NULL), time(NULL)); + crand_t rng = crand_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(); - cstc64_uniform_t dist = cstc64_uniform_init(0, M); + crand_uniform_t dist = crand_uniform_init(0, M); c_forrange (N) { - ++cmap_x_emplace(&map, cstc64_uniform(&rng, &dist), 0).first->second; + ++cmap_x_emplace(&map, crand_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 5303dd0c..a3842aae 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); - cstc64_t rng = cstc64_init(seed); - cstc64_normalf_t dist = cstc64_normalf_init(Mean, StdDev); + crand_t rng = crand_init(seed); + crand_normalf_t dist = crand_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( cstc64_normalf(&rng, &dist) ); + int index = (int) round( crand_normalf(&rng, &dist) ); cmap_i_emplace(&mhist, index, 0).first->second += 1; } diff --git a/examples/heap.c b/examples/heap.c index 166dd6ba..83f9d97f 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); - cstc64_t rng; + crand_t rng; int N = 3000000, M = 100; cpqueue_f pq = cpqueue_f_init(); - rng = cstc64_init(seed); + rng = crand_init(seed); clock_t start = clock(); c_forrange (i, int, N) - cvec_f_push_back(&pq, (float) cstc64_randf(&rng)*100000); + cvec_f_push_back(&pq, (float) crand_nextf(&rng)*100000); cpqueue_f_make_heap(&pq); printf("Built priority queue: %f secs\n", (clock() - start) / (float) CLOCKS_PER_SEC); @@ -35,7 +35,7 @@ int main() start = clock(); c_forrange (i, int, N) - cpqueue_f_push(&pq, (float) cstc64_randf(&rng)*100000); + cpqueue_f_push(&pq, (float) crand_nextf(&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 fbfa8830..4f6d16e1 100644 --- a/examples/list.c +++ b/examples/list.c @@ -9,11 +9,11 @@ int main() { const int n = 2000000; clist_fx list = clist_inits; - cstc64_t eng = cstc64_init(1234); - cstc64_uniformf_t dist = cstc64_uniformf_init(100.0f, n); + crand_t rng = crand_init(1234); + crand_uniformf_t dist = crand_uniformf_init(100.0f, n); int m = 0; c_forrange (i, int, n) - clist_fx_push_back(&list, cstc64_uniformf(&eng, &dist)), ++m; + clist_fx_push_back(&list, crand_uniformf(&rng, &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 e81cd8a3..f6b4ecf0 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; - cstc64_t pcg = cstc64_init(time(NULL)); - cstc64_uniform_t dist = cstc64_uniform_init(0, N * 10); + crand_t rng = crand_init(time(NULL)); + crand_uniform_t dist = crand_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, cstc64_uniform(&pcg, &dist)); + cpqueue_i_push(&heap, crand_uniform(&rng, &dist)); // push some negative numbers too. c_push_items(&heap, cpqueue_i, {-231, -32, -873, -4, -343}); c_forrange (N) - cpqueue_i_push(&heap, cstc64_uniform(&pcg, &dist)); + cpqueue_i_push(&heap, crand_uniform(&rng, &dist)); // Extract the hundred smallest. diff --git a/examples/queue.c b/examples/queue.c index 0ff29736..57362493 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; - cstc64_uniform_t dist; - cstc64_t rng = cstc64_init(1234); - dist = cstc64_uniform_init(0, n); + crand_uniform_t dist; + crand_t rng = crand_init(1234); + dist = crand_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, cstc64_uniform(&rng, &dist)); + cqueue_i_push(&queue, crand_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 = cstc64_uniform(&rng, &dist); + int r = crand_uniform(&rng, &dist); if (r & 1) ++n, cqueue_i_push(&queue, r); else diff --git a/examples/random.c b/examples/random.c index 0cf9da6f..de1fe778 100644 --- a/examples/random.c +++ b/examples/random.c @@ -7,18 +7,18 @@ int main() { enum {R = 30}; - const size_t N = 500000000; - uint64_t seed = time(NULL); - cstc64_t stc = cstc64_init(seed); + const size_t N = 1000000000; + uint64_t seed = 1234; // time(NULL); + crand_t rng = crand_init(seed); uint64_t sum = 0; - cstc64_normalf_t dist2 = cstc64_normalf_init(R / 2.0, R / 6.0); + crand_normalf_t dist2 = crand_normalf_init(R / 2.0, R / 6.0); size_t N2 = 10000000; int hist[R] = {0}; sum = 0; c_forrange (N2) { - int n = round((cstc64_normalf(&stc, &dist2) + 0.5)); + int n = round((crand_normalf(&rng, &dist2) + 0.5)); sum += n; if (n >= 0 && n < R) ++hist[n]; } @@ -33,16 +33,16 @@ int main() sum = 0; before = clock(); c_forrange (N) { - sum += cstc64_rand(&stc); + sum += crand_next(&rng); } diff = clock() - before; printf("random : %f secs, %zu %f\n", (float) diff / CLOCKS_PER_SEC, N, (double) sum / N); - cstc64_uniform_t dist1 = cstc64_uniform_init(0, 1000); + crand_uniform_t dist1 = crand_uniform_init(0, 1000); sum = 0; before = clock(); c_forrange (N) { - sum += cstc64_uniform(&stc, &dist1); + sum += crand_uniform(&rng, &dist1); } diff = clock() - before; printf("uniform: %f secs, %zu %f\n", (float) diff / CLOCKS_PER_SEC, N, (double) sum / N); diff --git a/examples/rngtest.c b/examples/rngtest.c index 76666484..db7a8093 100644 --- a/examples/rngtest.c +++ b/examples/rngtest.c @@ -13,14 +13,14 @@ int main(void) clock_t diff, before; uint64_t v; - 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); + crand_t rng = crand_init(time(NULL)); + crand_uniform_t idist = crand_uniform_init(10, 20); + crand_uniformf_t fdist = crand_uniformf_init(10, 20); before = clock(); v = 0; c_forrange (NN) { - v += cstc64_rand(&stc); + v += crand_next(&rng); } diff = clock() - before; printf("stc64_rand: %.02f, %zu\n", (float) diff / CLOCKS_PER_SEC, v); @@ -28,13 +28,13 @@ int main(void) before = clock(); v = 0; c_forrange (NN) { - v += cstc64_uniform(&stc, &idist); + v += crand_uniform(&rng, &idist); } diff = clock() - before; printf("stc64_uniform: %.02f, %zu\n\n", (float) diff / CLOCKS_PER_SEC, v); - c_forrange (30) printf("%02zd ", cstc64_uniform(&stc, &idist)); + c_forrange (30) printf("%02zd ", crand_uniform(&rng, &idist)); puts(""); - c_forrange (8) printf("%f ", cstc64_uniformf(&stc, &fdist)); + c_forrange (8) printf("%f ", crand_uniformf(&rng, &fdist)); puts(""); } \ No newline at end of file diff --git a/stc/clist.h b/stc/clist.h index f322fd00..6594aa8c 100644 --- a/stc/clist.h +++ b/stc/clist.h @@ -36,10 +36,10 @@ int main() { clist_ix list = clist_inits; - cstc64_t rng = cstc64_init(12345); + crand_t rng = crand_init(12345); int n; for (int i=0; i<1000000; ++i) // one million - clist_ix_push_back(&list, cstc64_rand(&rng) >> 32); + clist_ix_push_back(&list, crand_next(&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 2bdf2807..c1dfde17 100644 --- a/stc/cpqueue.h +++ b/stc/cpqueue.h @@ -29,13 +29,13 @@ using_cpqueue(f, cvec_f, >); // min-heap (increasing values) int main() { - cstc64_t gen = cstc64_init(1234); - cstc64_uniformf_t dist = cstc64_uniformf_init(10.0f, 100.0f); + crand_t rng = crand_init(1234); + crand_uniformf_t dist = crand_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, cstc64_uniformf(&gen, dist)); + cpqueue_f_push(&queue, crand_uniformf(&rng, 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 6226b59b..f385aa4a 100644 --- a/stc/cqueue.h +++ b/stc/cqueue.h @@ -32,19 +32,19 @@ int main() { int n = 10000000; - cstc64_t gen = cstc64_init(1234); - cstc64_uniform_t dist = cstc64_uniform_init(gen, 0, n); + crand_t rng = crand_init(1234); + crand_uniform_t dist = crand_uniform_init(rng, 0, n); cqueue_i queue = cqueue_i_init(); // Push ten million random numbers onto the queue. for (int i=0; i0; --i) { - int r = cstc64_uniform(&dist); + int r = crand_uniform(&dist); if (r & 1) ++n, cqueue_i_push(&queue, r); else diff --git a/stc/crandom.h b/stc/crandom.h index f0413abd..5db3c160 100644 --- a/stc/crandom.h +++ b/stc/crandom.h @@ -24,64 +24,64 @@ #define CRANDOM__H__ /* -// cstc64: Pseudo-random number generator +// crand: 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); + crand_t rng = crand_init(seed); + crand_uniform_t dist1 = crand_uniform_init(1, 6); + crand_uniformf_t dist2 = crand_uniformf_init(1.0, 10.0); + crand_normalf_t dist3 = crand_normalf_init(1.0, 10.0); - 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); + uint64_t i = crand_next(&rng); + int64_t iu = crand_uniform(&rng, &dist1); + double xu = crand_uniformf(&rng, &dist2); + double xn = crand_normalf(&rng, &dist3); } */ #include "ccommon.h" #include #include -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; +typedef struct {uint64_t state[4];} crand_t; +typedef struct {int64_t offset; uint64_t range, threshold;} crand_uniform_t; +typedef struct {double offset, range;} crand_uniformf_t; +typedef struct {double mean, stddev, next; bool has_next;} crand_normalf_t; /* int random number generator, range [0, 2^64). PRNG copyright Tyge Løvset, NORCE Research, 2020 */ -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); +STC_API crand_t crand_init(uint64_t seed); +STC_API crand_t crand_with_seq(uint64_t seed, uint64_t seq); +STC_API uint64_t crand_next(crand_t* rng); /* 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)}; +STC_INLINE double crand_nextf(crand_t* rng) { + union {uint64_t i; double f;} u = {0x3FF0000000000000ull | (crand_next(rng) >> 12)}; return u.f - 1.0; } /* 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_API crand_uniform_t crand_uniform_init(int64_t low, int64_t high); +STC_API int64_t crand_uniform(crand_t* rng, crand_uniform_t* 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 crand_uniformf_t crand_uniformf_init(double low, double high) { + crand_uniformf_t dist = {low, high - low}; return dist; } -STC_INLINE double cstc64_uniformf(cstc64_t* rng, cstc64_uniformf_t* dist) { - return dist->offset + cstc64_randf(rng) * dist->range; +STC_INLINE double crand_uniformf(crand_t* rng, crand_uniformf_t* dist) { + return dist->offset + crand_nextf(rng) * dist->range; } /* 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_INLINE crand_normalf_t crand_normalf_init(double mean, double stddev) { + crand_normalf_t dist = {mean, stddev, 0.0, false}; return dist; } -STC_API double cstc64_normalf(cstc64_t* rng, cstc64_normalf_t* dist); +STC_API double crand_normalf(crand_t* rng, crand_normalf_t* dist); #if !defined(STC_HEADER) || defined(STC_IMPLEMENTATION) -/* PRNG stc64: Tyge Løvset, NORCE Research, 2020. +/* PRNG crand: by 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. @@ -91,46 +91,58 @@ STC_API double cstc64_normalf(cstc64_t* rng, cstc64_normalf_t* dist); * and simple correlation tests, i.e. interleaved streams with one-bit diff state. */ -STC_DEF cstc64_t cstc64_init(uint64_t seed) { - return cstc64_with_seq(seed, 0x3504f333d3aa0b34); +STC_DEF crand_t crand_init(uint64_t seed) { + return crand_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); +STC_DEF crand_t crand_with_seq(uint64_t seed, uint64_t seq) { + crand_t rng = {{seed, seed, seed, (seq << 1u) | 1u}}; + for (int i = 0; i < 8; ++i) crand_next(&rng); return rng; } -STC_DEF uint64_t cstc64_rand(cstc64_t* rng) { +STC_DEF uint64_t crand_next(crand_t* rng) { enum {LROT = 24, RSHIFT = 11, LSHIFT = 3}; uint64_t *s = rng->state; - 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; + 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; return result; } /* 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)}; + +STC_DEF crand_uniform_t crand_uniform_init(int64_t low, int64_t high) { + crand_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); - } + +#if defined(__SIZEOF_INT128__) + #define cmul128(a, b, m) *(__uint128_t *)m = (__uint128_t)a * b +#elif defined(_MSC_VER) && defined(_WIN64) + #include + #define cmul128(a, b, m) m[0] = _umul128(a, b, &m[1]) +#elif defined(__x86_64__) + #define cmul128(a, b, m) \ + asm("mulq %[rhs]" : "=a" (m[0]), "=d" (m[1]) \ + : [lhs] "0" (a), [rhs] "rm" (b)) +#endif + +STC_DEF int64_t crand_uniform(crand_t* rng, crand_uniform_t* dist) { + uint64_t m[2]; + do { cmul128(crand_next(rng), dist->range, m); } while (m[0] < dist->threshold); + return dist->offset + m[1]; } /* Marsaglia polar method for gaussian distribution. */ -STC_DEF double cstc64_normalf(cstc64_t* rng, cstc64_normalf_t* dist) { +STC_DEF double crand_normalf(crand_t* rng, crand_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 * cstc64_randf(rng) - 1.0; - u2 = 2.0 * cstc64_randf(rng) - 1.0; + u1 = 2.0 * crand_nextf(rng) - 1.0; + u2 = 2.0 * crand_nextf(rng) - 1.0; s = u1*u1 + u2*u2; } while (s >= 1.0 || s == 0.0); m = sqrt(-2.0 * log(s) / s); @@ -138,6 +150,5 @@ STC_DEF double cstc64_normalf(cstc64_t* rng, cstc64_normalf_t* dist) { return (u1 * m) * dist->stddev + dist->mean; } - #endif #endif -- cgit v1.2.3