From 191b637299bf909d6b7b4db6e5d5e45964574a1c Mon Sep 17 00:00:00 2001 From: Tyge Løvset Date: Wed, 16 Dec 2020 13:31:13 +0100 Subject: Renamed crandom to crand, including renaming of relevant files. --- README.md | 2 +- docs/cpqueue_api.md | 2 +- docs/crand_api.md | 155 ++++++++++++++++++++++++++++++++++++++++++++++++ docs/crandom_api.md | 155 ------------------------------------------------ examples/benchmark.cpp | 2 +- examples/birthday.c | 2 +- examples/ex_gaussian.c | 2 +- examples/heap.c | 2 +- examples/list.c | 2 +- examples/priority.c | 2 +- examples/queue.c | 2 +- examples/random.c | 2 +- examples/rngtest.c | 2 +- stc/clist.h | 2 +- stc/cpqueue.h | 2 +- stc/cqueue.h | 2 +- stc/crand.h | 157 +++++++++++++++++++++++++++++++++++++++++++++++++ stc/crandom.h | 157 ------------------------------------------------- 18 files changed, 326 insertions(+), 326 deletions(-) create mode 100644 docs/crand_api.md delete mode 100644 docs/crandom_api.md create mode 100644 stc/crand.h delete mode 100644 stc/crandom.h diff --git a/README.md b/README.md index e2614973..0d61ace3 100644 --- a/README.md +++ b/README.md @@ -17,7 +17,7 @@ An elegant, fully typesafe, generic, customizable, user-friendly, consistent, an - [***cpqueue*** - A **priority queue** adapter type](docs/cpqueue_api.md) - [***cptr*** - Support for pointers and shared pointers in containers](docs/cptr_api.md) - [***copt*** - Implements *copt_get()*, a **getopt_long**-like function](docs/copt_api.md) -- [***crandom*** - A few very efficent modern **random number generators**](docs/crandom_api.md) +- [***crand*** - A few very efficent modern **random number generators**](docs/crand_api.md) - [***ccommon*** - Collection of general definitions](docs/ccommon_api.md) The usage of the containers is quite similar to the C++ standard containers, so it should be easy if you are familiar with them. diff --git a/docs/cpqueue_api.md b/docs/cpqueue_api.md index ac33e68f..09cca96d 100644 --- a/docs/cpqueue_api.md +++ b/docs/cpqueue_api.md @@ -55,7 +55,7 @@ void cpqueue_X_erase_at(cpqueue_X* self, size_t idx); ```c #include #include "stc/cpqueue.h" -#include "stc/crandom.h" +#include "stc/crand.h" using_cvec(i, int64_t); using_cpqueue(i, cvec_i, >); // adaptor type, '>' = min-heap diff --git a/docs/crand_api.md b/docs/crand_api.md new file mode 100644 index 00000000..1a26519c --- /dev/null +++ b/docs/crand_api.md @@ -0,0 +1,155 @@ +# Module crand: Pseudo Random Number Generators + +This describes the API of module **crand**. It contains a *64-bit PRNG*, and can generate +bounded uniform and normal distributed random numbers. + +**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\*\** on common platforms. 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 starts +a new unique 2^64 *minimum* length period. Note that for each Weyl-increment, the period +length is about 2^126 with a high probability. For a single thread, a minimum period of 2^127 +is generated when the Weyl-increment is incremented by 2 every 2^64 output number. + +**crand** passes *PractRand*, tested up to 8TB output, Vigna's Hamming weight test, and simple +correlation tests, i.e. N interleaved streams with only one-bit differences in initial state. + +## Types + +| Name | Type definition | Used to represent... | +|:-------------------|:------------------------------------------|:-----------------------------| +| `crand_t` | `struct {uint64_t state[4];}` | The PRNG engine type | +| `crand_uniform_t` | `struct {int64_t lower; uint64_t range;}` | Integer uniform distribution | +| `crand_uniformf_t` | `struct {double lower, range;}` | Real number uniform distr. | +| `crand_normalf_t` | `struct {double mean, stddev;}` | Normal distribution type | + +## Header file + +All cstr definitions and prototypes may be included in your C source file by including a single header file. +```c +#include "stc/crand.h" +``` + +## Methods + +```c + 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 fall within the range [*mean* - *stddev*, *mean* + *stddev*]. + +## Example +```c +#include +#include +#include +#include "stc/crand.h" +#include "stc/cstr.h" +#include "stc/cmap.h" +#include "stc/cvec.h" + +// Declare unordered map: int -> int with typetag 'i'. +using_cmap(i, int, size_t); + +// Comparison of map keys. +static int compare(cmap_i_entry_t *a, cmap_i_entry_t *b) { + return c_default_compare(&a->first, &b->first); +} +// Declare vector of map entries, with comparison function. +using_cvec(e, cmap_i_entry_t, c_default_del, compare); + + +int main() +{ + enum {N = 10000000}; + const double Mean = -12.0, StdDev = 6.0, Scale = 74; + + printf("Demo of gaussian / normal distribution of %d random samples\n", N); + + // Setup random engine with normal distribution. + uint64_t seed = time(NULL); + 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( crand_normalf(&rng, &dist) ); + cmap_i_emplace(&mhist, index, 0).first->second += 1; + } + + // Transfer map to vec and sort it by map entry keys. + cvec_e vhist = cvec_e_init(); + c_foreach (i, cmap_i, mhist) + cvec_e_push_back(&vhist, *i.val); + cvec_e_sort(&vhist); + + // Print the gaussian bar chart + cstr_t bar = cstr_init(); + c_foreach (i, cvec_e, vhist) { + size_t n = (size_t) (i.val->second * StdDev * Scale * 2.5 / N); + if (n > 0) { + cstr_resize(&bar, n, '*'); + printf("%4d %s\n", i.val->first, bar.str); + } + } + + // Cleanup + cstr_del(&bar); + cmap_i_del(&mhist); + cvec_e_del(&vhist); +} +``` +Output: +``` +Demo of gaussian / normal distribution of 10000000 random samples + -29 * + -28 ** + -27 *** + -26 **** + -25 ******* + -24 ********* + -23 ************* + -22 ****************** + -21 *********************** + -20 ****************************** + -19 ************************************* + -18 ******************************************** + -17 **************************************************** + -16 *********************************************************** + -15 ***************************************************************** + -14 ********************************************************************* + -13 ************************************************************************ + -12 ************************************************************************* + -11 ************************************************************************ + -10 ********************************************************************* + -9 ***************************************************************** + -8 *********************************************************** + -7 **************************************************** + -6 ******************************************** + -5 ************************************* + -4 ****************************** + -3 *********************** + -2 ****************** + -1 ************* + 0 ********* + 1 ******* + 2 **** + 3 *** + 4 ** + 5 * +``` diff --git a/docs/crandom_api.md b/docs/crandom_api.md deleted file mode 100644 index 5fc89e78..00000000 --- a/docs/crandom_api.md +++ /dev/null @@ -1,155 +0,0 @@ -# Module crandom: Pseudo Random Number Generators - -This describes the API of module **crand**. It contains a *64-bit PRNG*, and can generate -bounded uniform and normal distributed random numbers. - -**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\*\** on common platforms. 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 starts -a new unique 2^64 *minimum* length period. Note that for each Weyl-increment, the period -length is about 2^126 with a high probability. For a single thread, a minimum period of 2^127 -is generated when the Weyl-increment is incremented by 2 every 2^64 output number. - -**crand** passes *PractRand*, tested up to 8TB output, Vigna's Hamming weight test, and simple -correlation tests, i.e. N interleaved streams with only one-bit differences in initial state. - -## Types - -| Name | Type definition | Used to represent... | -|:-------------------|:------------------------------------------|:-----------------------------| -| `crand_t` | `struct {uint64_t state[4];}` | The PRNG engine type | -| `crand_uniform_t` | `struct {int64_t lower; uint64_t range;}` | Integer uniform distribution | -| `crand_uniformf_t` | `struct {double lower, range;}` | Real number uniform distr. | -| `crand_normalf_t` | `struct {double mean, stddev;}` | Normal distribution type | - -## Header file - -All cstr definitions and prototypes may be included in your C source file by including a single header file. -```c -#include "stc/crandom.h" -``` - -## Methods - -```c - 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 fall within the range [*mean* - *stddev*, *mean* + *stddev*]. - -## Example -```c -#include -#include -#include -#include "stc/crandom.h" -#include "stc/cstr.h" -#include "stc/cmap.h" -#include "stc/cvec.h" - -// Declare unordered map: int -> int with typetag 'i'. -using_cmap(i, int, size_t); - -// Comparison of map keys. -static int compare(cmap_i_entry_t *a, cmap_i_entry_t *b) { - return c_default_compare(&a->first, &b->first); -} -// Declare vector of map entries, with comparison function. -using_cvec(e, cmap_i_entry_t, c_default_del, compare); - - -int main() -{ - enum {N = 10000000}; - const double Mean = -12.0, StdDev = 6.0, Scale = 74; - - printf("Demo of gaussian / normal distribution of %d random samples\n", N); - - // Setup random engine with normal distribution. - uint64_t seed = time(NULL); - 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( crand_normalf(&rng, &dist) ); - cmap_i_emplace(&mhist, index, 0).first->second += 1; - } - - // Transfer map to vec and sort it by map entry keys. - cvec_e vhist = cvec_e_init(); - c_foreach (i, cmap_i, mhist) - cvec_e_push_back(&vhist, *i.val); - cvec_e_sort(&vhist); - - // Print the gaussian bar chart - cstr_t bar = cstr_init(); - c_foreach (i, cvec_e, vhist) { - size_t n = (size_t) (i.val->second * StdDev * Scale * 2.5 / N); - if (n > 0) { - cstr_resize(&bar, n, '*'); - printf("%4d %s\n", i.val->first, bar.str); - } - } - - // Cleanup - cstr_del(&bar); - cmap_i_del(&mhist); - cvec_e_del(&vhist); -} -``` -Output: -``` -Demo of gaussian / normal distribution of 10000000 random samples - -29 * - -28 ** - -27 *** - -26 **** - -25 ******* - -24 ********* - -23 ************* - -22 ****************** - -21 *********************** - -20 ****************************** - -19 ************************************* - -18 ******************************************** - -17 **************************************************** - -16 *********************************************************** - -15 ***************************************************************** - -14 ********************************************************************* - -13 ************************************************************************ - -12 ************************************************************************* - -11 ************************************************************************ - -10 ********************************************************************* - -9 ***************************************************************** - -8 *********************************************************** - -7 **************************************************** - -6 ******************************************** - -5 ************************************* - -4 ****************************** - -3 *********************** - -2 ****************** - -1 ************* - 0 ********* - 1 ******* - 2 **** - 3 *** - 4 ** - 5 * -``` diff --git a/examples/benchmark.cpp b/examples/benchmark.cpp index 206da5ce..e919bc3f 100644 --- a/examples/benchmark.cpp +++ b/examples/benchmark.cpp @@ -1,6 +1,6 @@ #include #include -#include +#include #include #include #include "others/khash.h" diff --git a/examples/birthday.c b/examples/birthday.c index ac05a5f5..a13fa8c9 100644 --- a/examples/birthday.c +++ b/examples/birthday.c @@ -2,7 +2,7 @@ #include #include -#include +#include #include #include diff --git a/examples/ex_gaussian.c b/examples/ex_gaussian.c index a3842aae..0fed913f 100644 --- a/examples/ex_gaussian.c +++ b/examples/ex_gaussian.c @@ -1,7 +1,7 @@ #include #include #include -#include "stc/crandom.h" +#include "stc/crand.h" #include "stc/cstr.h" #include "stc/cmap.h" #include "stc/cvec.h" diff --git a/examples/heap.c b/examples/heap.c index 83f9d97f..b45e5ada 100644 --- a/examples/heap.c +++ b/examples/heap.c @@ -1,6 +1,6 @@ #include #include -#include +#include #include #include diff --git a/examples/list.c b/examples/list.c index 4f6d16e1..a838b3bd 100644 --- a/examples/list.c +++ b/examples/list.c @@ -1,7 +1,7 @@ #include #include #include -#include +#include using_clist(fx, double); int main() { diff --git a/examples/priority.c b/examples/priority.c index f6b4ecf0..40fb395e 100644 --- a/examples/priority.c +++ b/examples/priority.c @@ -4,7 +4,7 @@ #include #include #include -#include +#include using_cvec(i, int64_t); using_cpqueue(i, cvec_i, >); // min-heap (increasing values) diff --git a/examples/queue.c b/examples/queue.c index 57362493..2499a9ab 100644 --- a/examples/queue.c +++ b/examples/queue.c @@ -1,4 +1,4 @@ -#include +#include #include #include diff --git a/examples/random.c b/examples/random.c index de1fe778..94a68607 100644 --- a/examples/random.c +++ b/examples/random.c @@ -1,7 +1,7 @@ #include #include #include -#include +#include #include int main() diff --git a/examples/rngtest.c b/examples/rngtest.c index db7a8093..f52099c7 100644 --- a/examples/rngtest.c +++ b/examples/rngtest.c @@ -1,6 +1,6 @@ #include #include -#include +#include #ifdef __cplusplus #include #endif diff --git a/stc/clist.h b/stc/clist.h index 6594aa8c..838edf07 100644 --- a/stc/clist.h +++ b/stc/clist.h @@ -31,7 +31,7 @@ #include #include - #include + #include using_clist(ix, int64_t); int main() { diff --git a/stc/cpqueue.h b/stc/cpqueue.h index c1dfde17..12834a47 100644 --- a/stc/cpqueue.h +++ b/stc/cpqueue.h @@ -23,7 +23,7 @@ /* Priority-Queue adapter (implemented as heap), default uses cvec. - #include + #include #include using_cvec(f, float); using_cpqueue(f, cvec_f, >); // min-heap (increasing values) diff --git a/stc/cqueue.h b/stc/cqueue.h index f385aa4a..f16e1aeb 100644 --- a/stc/cqueue.h +++ b/stc/cqueue.h @@ -25,7 +25,7 @@ /* Queue adapter, default uses clist. - #include + #include #include using_clist(i, int); using_cqueue(i, clist_i); diff --git a/stc/crand.h b/stc/crand.h new file mode 100644 index 00000000..54fbcc10 --- /dev/null +++ b/stc/crand.h @@ -0,0 +1,157 @@ +/* MIT License + * + * Copyright (c) 2020 Tyge Løvset, NORCE, www.norceresearch.no + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in all + * copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + */ +#ifndef CRANDOM__H__ +#define CRANDOM__H__ + +/* +// crand: Pseudo-random number generator +#include "stc/crand.h" +int main() { + uint64_t seed = 123456789; + 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 = 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];} crand_t; +typedef struct {int64_t lower; uint64_t range, threshold;} crand_uniform_t; +typedef struct {double lower, 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 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 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 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 crand_uniformf_t crand_uniformf_init(double low, double high) { + crand_uniformf_t dist = {low, high - low}; return dist; +} +STC_INLINE double crand_uniformf(crand_t* rng, crand_uniformf_t* dist) { + return crand_nextf(rng)*dist->range + dist->lower; +} + +/* double normal distributed RNG. */ +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 crand_normalf(crand_t* rng, crand_normalf_t* dist); + + +#if !defined(STC_HEADER) || defined(STC_IMPLEMENTATION) + +/* 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. + * 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. + */ + +STC_DEF crand_t crand_init(uint64_t seed) { + return crand_with_seq(seed, 0x3504f333d3aa0b34); +} +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 crand_next(crand_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; + return result; +} + +/* Very fast unbiased uniform int RNG with bounds [low, high] */ +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; +} + +#if defined(__SIZEOF_INT128__) + #define cmul128(a, b, lo, hi) \ + do { __uint128_t _z = (__uint128_t)(a) * (b); \ + *(lo) = (uint64_t)_z, *(hi) = _z >> 64; } while(0) +#elif defined(_MSC_VER) && defined(_WIN64) + #include + #define cmul128(a, b, lo, hi) (*(lo) = _umul128(a, b, hi), (void)0) +#elif defined(__x86_64__) + #define cmul128(a, b, lo, hi) \ + asm("mulq %[rhs]" : "=a" (*(lo)), "=d" (*(hi)) \ + : [lhs] "0" (a), [rhs] "rm" (b)) +#endif + +STC_DEF int64_t crand_uniform(crand_t* rng, crand_uniform_t* d) { + uint64_t lo, hi; + do { + cmul128(crand_next(rng), d->range, &lo, &hi); + } while (lo < d->threshold); + return d->lower + hi; +} + +/* Marsaglia polar method for gaussian/normal distribution. */ +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 * 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); + dist->next = u2 * m, dist->has_next = true; + return (u1 * m) * dist->stddev + dist->mean; +} + +#endif +#endif diff --git a/stc/crandom.h b/stc/crandom.h deleted file mode 100644 index 8a67e3ce..00000000 --- a/stc/crandom.h +++ /dev/null @@ -1,157 +0,0 @@ -/* MIT License - * - * Copyright (c) 2020 Tyge Løvset, NORCE, www.norceresearch.no - * - * Permission is hereby granted, free of charge, to any person obtaining a copy - * of this software and associated documentation files (the "Software"), to deal - * in the Software without restriction, including without limitation the rights - * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell - * copies of the Software, and to permit persons to whom the Software is - * furnished to do so, subject to the following conditions: - * - * The above copyright notice and this permission notice shall be included in all - * copies or substantial portions of the Software. - * - * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR - * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, - * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE - * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER - * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, - * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE - * SOFTWARE. - */ -#ifndef CRANDOM__H__ -#define CRANDOM__H__ - -/* -// crand: Pseudo-random number generator -#include "stc/crandom.h" -int main() { - uint64_t seed = 123456789; - 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 = 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];} crand_t; -typedef struct {int64_t lower; uint64_t range, threshold;} crand_uniform_t; -typedef struct {double lower, 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 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 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 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 crand_uniformf_t crand_uniformf_init(double low, double high) { - crand_uniformf_t dist = {low, high - low}; return dist; -} -STC_INLINE double crand_uniformf(crand_t* rng, crand_uniformf_t* dist) { - return crand_nextf(rng)*dist->range + dist->lower; -} - -/* double normal distributed RNG. */ -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 crand_normalf(crand_t* rng, crand_normalf_t* dist); - - -#if !defined(STC_HEADER) || defined(STC_IMPLEMENTATION) - -/* 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. - * 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. - */ - -STC_DEF crand_t crand_init(uint64_t seed) { - return crand_with_seq(seed, 0x3504f333d3aa0b34); -} -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 crand_next(crand_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; - return result; -} - -/* Very fast unbiased uniform int RNG with bounds [low, high] */ -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; -} - -#if defined(__SIZEOF_INT128__) - #define cmul128(a, b, lo, hi) \ - do { __uint128_t _z = (__uint128_t)(a) * (b); \ - *(lo) = (uint64_t)_z, *(hi) = _z >> 64; } while(0) -#elif defined(_MSC_VER) && defined(_WIN64) - #include - #define cmul128(a, b, lo, hi) (*(lo) = _umul128(a, b, hi), (void)0) -#elif defined(__x86_64__) - #define cmul128(a, b, lo, hi) \ - asm("mulq %[rhs]" : "=a" (*(lo)), "=d" (*(hi)) \ - : [lhs] "0" (a), [rhs] "rm" (b)) -#endif - -STC_DEF int64_t crand_uniform(crand_t* rng, crand_uniform_t* d) { - uint64_t lo, hi; - do { - cmul128(crand_next(rng), d->range, &lo, &hi); - } while (lo < d->threshold); - return d->lower + hi; -} - -/* Marsaglia polar method for gaussian/normal distribution. */ -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 * 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); - dist->next = u2 * m, dist->has_next = true; - return (u1 * m) * dist->stddev + dist->mean; -} - -#endif -#endif -- cgit v1.2.3