From 74a75b3a8d155ea155358674d3bf8251658183a5 Mon Sep 17 00:00:00 2001 From: Tyge Løvset Date: Mon, 28 Dec 2020 18:18:48 +0100 Subject: Reverted crand.h API back to rename crand_* to stc64_*, and crand_next() -> stc64_rand(). --- README.md | 2 +- benchmarks/cdeq_benchmark.cpp | 10 +-- benchmarks/cmap_benchmark.cpp | 6 +- benchmarks/cpque_benchmark.cpp | 8 +-- benchmarks/crand_benchmark.cpp | 138 +++++++++++++++++++++++++--------------- benchmarks/crand_benchmark2.cpp | 68 ++++++++++++++++---- benchmarks/vector_vs_deque.cpp | 24 +++---- docs/cpque_api.md | 6 +- docs/crand_api.md | 49 +++++++------- examples/birthday.c | 8 +-- examples/ex_gaussian.c | 6 +- examples/list.c | 6 +- examples/priority.c | 8 +-- examples/queue.c | 10 +-- examples/random.c | 12 ++-- stc/clist.h | 4 +- stc/cpque.h | 6 +- stc/cqueue.h | 8 +-- stc/crand.h | 76 +++++++++++----------- 19 files changed, 263 insertions(+), 192 deletions(-) diff --git a/README.md b/README.md index 61c814fb..abf081e6 100644 --- a/README.md +++ b/README.md @@ -6,7 +6,7 @@ Introduction An modern, fully typesafe, generic, customizable, user-friendly, consistent, and very fast standard container library for C99. This is a small headers only library with the most used container components, and a few algorithms: - [***cstr*** - Powerful and compact **string** type](docs/cstr_api.md) -- [***crand*** - An extremely efficent modern **random number generator**](docs/crand_api.md) +- [***crand*** - An extremely efficent modern **random number generator**](docs/stc64_api.md) - [***clist*** - Templated **std::forward_list** alike type](docs/clist_api.md) - [***cmap*** - Templated **std::unordered_map** alike type](docs/cmap_api.md) - [***cset*** - Templated **std::unordered_set** alike type](docs/cset_api.md) diff --git a/benchmarks/cdeq_benchmark.cpp b/benchmarks/cdeq_benchmark.cpp index fd01ad8b..6a376d2a 100644 --- a/benchmarks/cdeq_benchmark.cpp +++ b/benchmarks/cdeq_benchmark.cpp @@ -4,15 +4,15 @@ #include #include -enum {N = 200000000, M = 10000, P = 5000, R = 100}; +enum {N = 200000000, M = 10000, P = 5000, R = 50}; using_cdeq(i, int); void test1() { clock_t t1 = clock(), t2, t3; - crand_t rng = crand_init(0); + stc64_t rng = stc64_init(0); std::deque deq; for (size_t i = 1; i < N; i++) { - deq.push_front(crand_next(&rng)); + deq.push_front(stc64_rand(&rng)); if (i % M == 0) for (int j = 0; j < P; j++) deq.pop_back(); @@ -30,10 +30,10 @@ void test1() { void test2() { clock_t t1 = clock(), t2, t3; - crand_t rng = crand_init(0); + stc64_t rng = stc64_init(0); cdeq_i deq = cdeq_inits; for (size_t i = 1; i < N; i++) { - cdeq_i_push_front(&deq, crand_next(&rng)); + cdeq_i_push_front(&deq, stc64_rand(&rng)); if (i % M == 0) for (int j = 0; j < P; j++) cdeq_i_pop_back(&deq); diff --git a/benchmarks/cmap_benchmark.cpp b/benchmarks/cmap_benchmark.cpp index aaf9fc73..d9b6d1ad 100644 --- a/benchmarks/cmap_benchmark.cpp +++ b/benchmarks/cmap_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_t rng; -#define SEED(s) rng = crand_init(seed) -#define RAND(N) (crand_next(&rng) & ((1 << N) - 1)) +stc64_t rng; +#define SEED(s) rng = stc64_init(seed) +#define RAND(N) (stc64_rand(&rng) & ((1 << N) - 1)) #define CMAP_SETUP(X, Key, Value) cmap_##X map = cmap_inits \ diff --git a/benchmarks/cpque_benchmark.cpp b/benchmarks/cpque_benchmark.cpp index 625bb056..396d763d 100644 --- a/benchmarks/cpque_benchmark.cpp +++ b/benchmarks/cpque_benchmark.cpp @@ -10,15 +10,15 @@ using_cpque(f, cvec_f, >); int main() { uint32_t seed = time(NULL); - crand_t rng; + stc64_t rng; int N = 10000000, M = 10; cpque_f pq = cpque_f_init(); - rng = crand_init(seed); + rng = stc64_init(seed); clock_t start = clock(); c_forrange (i, int, N) - cvec_f_push_back(&pq, (float) crand_nextf(&rng)*100000); + cvec_f_push_back(&pq, (float) stc64_randf(&rng)*100000); cpque_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) - cpque_f_push(&pq, (float) crand_nextf(&rng)*100000); + cpque_f_push(&pq, (float) stc64_randf(&rng)*100000); printf("pushed PQ: %f secs\n", (clock() - start) / (float) CLOCKS_PER_SEC); c_forrange (i, int, M) { diff --git a/benchmarks/crand_benchmark.cpp b/benchmarks/crand_benchmark.cpp index 78099905..70a3bb45 100644 --- a/benchmarks/crand_benchmark.cpp +++ b/benchmarks/crand_benchmark.cpp @@ -4,15 +4,15 @@ #include static inline uint64_t rotl64(const uint64_t x, const int k) - { return (x << k) | (x >> (64 - k)); } + { return (x << k) | (x >> (64 - k)); } static uint64_t splitmix64_x = 87213627321ull; /* The state can be seeded with any value. */ uint64_t splitmix64(void) { - uint64_t z = (splitmix64_x += 0x9e3779b97f4a7c15); - z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9; - z = (z ^ (z >> 27)) * 0x94d049bb133111eb; - return z ^ (z >> 31); + uint64_t z = (splitmix64_x += 0x9e3779b97f4a7c15); + z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9; + z = (z ^ (z >> 27)) * 0x94d049bb133111eb; + return z ^ (z >> 31); } static void init_state(uint64_t *rng, uint64_t seed) { @@ -23,26 +23,47 @@ static void init_state(uint64_t *rng, uint64_t seed) { /* jsf64 */ static inline uint64_t jsf64(uint64_t *s) { - uint64_t e = s[0] - rotl64(s[1], 7); - s[0] = s[1] ^ rotl64(s[2], 13); - s[1] = s[2] + rotl64(s[3], 37); - s[2] = s[3] + e; - s[3] = e + s[0]; - return s[3]; + uint64_t e = s[0] - rotl64(s[1], 7); + s[0] = s[1] ^ rotl64(s[2], 13); + s[1] = s[2] + rotl64(s[3], 37); + s[2] = s[3] + e; + s[3] = e + s[0]; + return s[3]; } +/* sfc64 */ + +static inline uint64_t sfc64(uint64_t *s) { + uint64_t result = s[0] + s[1] + s[3]++; + s[0] = s[1] ^ (s[1] >> 11); + s[1] = s[2] + (s[2] << 3); + s[2] = rotl64(s[2], 24) + result; + return result; +} + +/* sfc64 with Weyl increment */ +static uint64_t weyl = 1234566789123ull; +static inline uint64_t sfc64w(uint64_t *s) { + uint64_t result = s[0] + s[1] + (s[3] += weyl|1); + s[0] = s[1] ^ (s[1] >> 11); + s[1] = s[2] + (s[2] << 3); + s[2] = rotl64(s[2], 24) + result; + return result; +} + + /* xoshiro256** */ static inline uint64_t xoshiro256starstar(uint64_t* s) { - const uint64_t result = rotl64(s[1] * 5, 7) * 9; - const uint64_t t = s[1] << 17; - s[2] ^= s[0]; - s[3] ^= s[1]; - s[1] ^= s[2]; - s[0] ^= s[3]; - s[2] ^= t; - s[3] = rotl64(s[3], 45); - return result; + const uint64_t result = rotl64(s[1] * 5, 7) * 9; + const uint64_t t = s[1] << 17; + s[2] ^= s[0]; + s[3] ^= s[1]; + s[1] ^= s[2]; + s[0] ^= s[3]; + s[2] ^= t; + s[3] = rotl64(s[3], 45); + return result; } @@ -63,92 +84,103 @@ using namespace std; int main(void) { - enum {N = 524288000}; + enum {N = 1000000000}; uint64_t* recipient = new uint64_t[N]; - static crand_t rng; + static stc64_t rng; init_state(rng.state, 12345123); + cout << "WARMUP" << endl; + for (size_t i = 0; i < N; i++) + recipient[i] = wyhash64(rng.state); + clock_t beg, end; for (size_t ti = 0; ti < 4; ti++) { + cout << endl << "ROUND " << ti+1 << endl; beg = clock(); for (size_t i = 0; i < N; i++) recipient[i] = wyhash64(rng.state); end = clock(); - cerr << "ROUND " << ti+1 << endl - << "wyhash64:\t" + cout << "wyhash64:\t" << (float(end - beg) / CLOCKS_PER_SEC) - << " s" << endl; - cout << "bogus:" << recipient[312] << endl; + << " s: " << recipient[312] << endl; beg = clock(); for (size_t i = 0; i < N; i++) - recipient[i] = crand_next(&rng); + recipient[i] = sfc64w(rng.state); end = clock(); - cerr << "stc crand:\t" + cout << "sfc64w:\t\t" << (float(end - beg) / CLOCKS_PER_SEC) - << " s" << endl; - cout << "bogus:" << recipient[312] << endl; + << " s: " << recipient[312] << endl; + + beg = clock(); + for (size_t i = 0; i < N; i++) + recipient[i] = stc64_rand(&rng); + end = clock(); + cout << "stc-crand:\t" + << (float(end - beg) / CLOCKS_PER_SEC) + << " s: " << recipient[312] << endl; beg = clock(); for (size_t i = 0; i < N; i++) recipient[i] = xoshiro256starstar(rng.state); end = clock(); - cerr << "xoshiro256**:\t" + cout << "xoshiro256**:\t" << (float(end - beg) / CLOCKS_PER_SEC) - << " s" << endl; - cout << "bogus:" << recipient[312] << endl; + << " s: " << recipient[312] << endl; beg = clock(); for (size_t i = 0; i < N; i++) recipient[i] = lehmer64(rng.state); end = clock(); - cerr << "lehmer64:\t" + cout << "lehmer64:\t" << ((float) end - beg) / CLOCKS_PER_SEC - << " s" << endl; - cout << "bogus:" << recipient[312] << endl; - - + << " s: " << recipient[312] << endl; - cout << endl - << "Next we do random number computations only, doing no work." + cout << "Next we do random number computations only, doing no work." << endl; uint64_t s = 0; beg = clock(); for (size_t i = 0; i < N; i++) s += wyhash64(rng.state); end = clock(); - cerr << "wyhash64:\t" + cout << "wyhash64:\t" + << ((float) end - beg) / CLOCKS_PER_SEC + << " s: " << s << endl; + + s = 0; + beg = clock(); + for (size_t i = 0; i < N; i++) + s += sfc64w(rng.state); + end = clock(); + cout << "sfc64w:\t\t" << ((float) end - beg) / CLOCKS_PER_SEC - << " s" << endl; - cout << "bogus:" << s << endl; + << " s: " << s << endl; s = 0; beg = clock(); for (size_t i = 0; i < N; i++) - s += crand_next(&rng); + s += stc64_rand(&rng); end = clock(); - cerr << "stc crand:\t" + cout << "stc-crand:\t" << ((float) end - beg) / CLOCKS_PER_SEC - << " s" << endl; - cout << "bogus:" << s << endl; + << " s: " << s << endl; s = 0; beg = clock(); for (size_t i = 0; i < N; i++) s += xoshiro256starstar(rng.state); end = clock(); - cerr << "xoshiro256**:\t" + cout << "xoshiro256**:\t" << ((float) end - beg) / CLOCKS_PER_SEC - << " s" << endl; - cout << "bogus:" << s << endl; + << " s: " << s << endl; + s = 0; beg = clock(); for (size_t i = 0; i < N; i++) s += lehmer64(rng.state); end = clock(); - cerr << "lehmer64:\t" + cout << "lehmer64:\t" << ((float) end - beg) / CLOCKS_PER_SEC - << " s" << endl; - cout << "bogus:" << s << endl << endl; + << " s: " << s << endl; } delete[] recipient; return 0; diff --git a/benchmarks/crand_benchmark2.cpp b/benchmarks/crand_benchmark2.cpp index ac7296fc..966a3675 100644 --- a/benchmarks/crand_benchmark2.cpp +++ b/benchmarks/crand_benchmark2.cpp @@ -1,11 +1,13 @@ #include #include -#include #include +#include "stc/crand.h" +#include "others/pcg_random.hpp" + +enum {N = 1000000000}; void test1(void) { - enum {N = 1000000000}; clock_t diff, before; uint64_t sum; @@ -20,7 +22,7 @@ void test1(void) sum += rng(); } diff = clock() - before; - printf("std::random:\t\t%.02f, %zu\n", (float) diff / CLOCKS_PER_SEC, sum); + printf("std::random:\t\t%.02f, %zu, sz:%zu\n", (float) diff / CLOCKS_PER_SEC, sum, sizeof rng); before = clock(); sum = 0; @@ -36,35 +38,72 @@ void test1(void) puts("\n"); } -void test2(void) +void test2() +{ + clock_t diff, before; + uint64_t sum; + + // Seed with a real random value, if available + pcg_extras::seed_seq_from seed_source; + + // Make a random number engine + pcg64 rng(seed_source); + + // Choose a random mean between 1 and 10 + std::uniform_int_distribution idist(1, 10); + std::uniform_real_distribution fdist(1, 10); + + before = clock(); + sum = 0; + c_forrange (N) { + sum += rng(); + } + diff = clock() - before; + printf("pcg64::random:\t\t%.02f, %zu, sz:%zu\n", (float) diff / CLOCKS_PER_SEC, sum, sizeof rng); + + before = clock(); + sum = 0; + c_forrange (N) { + sum += idist(rng); + } + diff = clock() - before; + printf("pcg64::uniform:\t\t%.02f, %zu\n\n", (float) diff / CLOCKS_PER_SEC, sum); + + c_forrange (30) printf("%02d ", idist(rng)); + puts(""); + c_forrange (8) printf("%f ", fdist(rng)); + puts("\n"); +} + + +void test3(void) { - enum {N = 1000000000}; clock_t diff, before; uint64_t sum; - crand_t rng = crand_init(time(NULL)); - crand_uniform_t idist = crand_uniform_init(1, 10); - crand_uniformf_t fdist = crand_uniformf_init(1, 10); + stc64_t rng = stc64_init(time(NULL)); + stc64_uniform_t idist = stc64_uniform_init(1, 10); + stc64_uniformf_t fdist = stc64_uniformf_init(1, 10); before = clock(); sum = 0; c_forrange (N) { - sum += crand_next(&rng); + sum += stc64_rand(&rng); } diff = clock() - before; - printf("crand_next:\t\t%.02f, %zu\n", (float) diff / CLOCKS_PER_SEC, sum); + printf("stc64_random:\t\t%.02f, %zu sz:%zu\n", (float) diff / CLOCKS_PER_SEC, sum, sizeof rng); before = clock(); sum = 0; c_forrange (N) { - sum += crand_uniform(&rng, &idist); + sum += stc64_uniform(&rng, &idist); } diff = clock() - before; - printf("crand_uniform:\t\t%.02f, %zu\n\n", (float) diff / CLOCKS_PER_SEC, sum); + printf("stc64_uniform:\t\t%.02f, %zu\n\n", (float) diff / CLOCKS_PER_SEC, sum); - c_forrange (30) printf("%02zd ", crand_uniform(&rng, &idist)); + c_forrange (30) printf("%02zd ", stc64_uniform(&rng, &idist)); puts(""); - c_forrange (8) printf("%f ", crand_uniformf(&rng, &fdist)); + c_forrange (8) printf("%f ", stc64_uniformf(&rng, &fdist)); puts("\n"); } @@ -72,4 +111,5 @@ int main() { test1(); test2(); + test3(); } \ No newline at end of file diff --git a/benchmarks/vector_vs_deque.cpp b/benchmarks/vector_vs_deque.cpp index ae37f246..6c26a823 100644 --- a/benchmarks/vector_vs_deque.cpp +++ b/benchmarks/vector_vs_deque.cpp @@ -20,7 +20,7 @@ void add(cvec_si* tm, const char* s, int n) { Si si = {s, n}; cvec_si_push_back( void test_vector(const int num_iterations) { std::vector v; - crand_t rng = crand_init(0); + stc64_t rng = stc64_init(0); v.reserve(num_iterations + 2); //Ensure there is enough space reserved. // == PUSH_BACK @@ -29,7 +29,7 @@ void test_vector(const int num_iterations) for (int i=0; i<(num_iterations>>1); ++i) { - v.push_back(crand_next(&rng)); + v.push_back(stc64_rand(&rng)); } clock_t t2 = std::clock(); @@ -44,7 +44,7 @@ void test_vector(const int num_iterations) for (int i=0; i<(num_iterations>>1); ++i) { // Rather add some more elements to back. - v.push_back(crand_next(&rng)); + v.push_back(stc64_rand(&rng)); } clock_t t2 = std::clock(); @@ -96,7 +96,7 @@ void test_vector(const int num_iterations) void test_deque(const int num_iterations) { std::deque d; - crand_t rng = crand_init(0); + stc64_t rng = stc64_init(0); // == PUSH_BACK { @@ -104,7 +104,7 @@ void test_deque(const int num_iterations) for (int i=0; i<(num_iterations>>1); ++i) { - d.push_back(crand_next(&rng)); + d.push_back(stc64_rand(&rng)); } clock_t t2 = std::clock(); @@ -118,7 +118,7 @@ void test_deque(const int num_iterations) for (int i=0; i<(num_iterations>>1); ++i) { - d.push_front(crand_next(&rng)); + d.push_front(stc64_rand(&rng)); } clock_t t2 = std::clock(); @@ -177,7 +177,7 @@ using_cvec(i, int); void test_cvec(const int num_iterations) { cvec_i v = cvec_inits; - crand_t rng = crand_init(0); + stc64_t rng = stc64_init(0); //v.reserve(num_iterations + 2); //Ensure there is enough space reserved. // == PUSH_BACK @@ -186,7 +186,7 @@ void test_cvec(const int num_iterations) for (int i=0; i<(num_iterations>>1); ++i) { - cvec_i_push_back(&v, crand_next(&rng)); + cvec_i_push_back(&v, stc64_rand(&rng)); } clock_t t2 = std::clock(); @@ -201,7 +201,7 @@ void test_cvec(const int num_iterations) for (int i=0; i<(num_iterations>>1); ++i) { // Rather add some more elements to back. - cvec_i_push_back(&v, crand_next(&rng)); + cvec_i_push_back(&v, stc64_rand(&rng)); } clock_t t2 = std::clock(); @@ -257,7 +257,7 @@ using_cdeq(i, int); void test_cdeq(const int num_iterations) { cdeq_i d = cdeq_i_with_capacity(num_iterations + 2); - crand_t rng = crand_init(0); + stc64_t rng = stc64_init(0); // == PUSH_BACK { @@ -265,7 +265,7 @@ void test_cdeq(const int num_iterations) for (int i=0; i<(num_iterations>>1); ++i) { - cdeq_i_push_back(&d, crand_next(&rng)); + cdeq_i_push_back(&d, stc64_rand(&rng)); } clock_t t2 = std::clock(); @@ -279,7 +279,7 @@ void test_cdeq(const int num_iterations) for (int i=0; i<(num_iterations>>1); ++i) { - cdeq_i_push_front(&d, crand_next(&rng)); + cdeq_i_push_front(&d, stc64_rand(&rng)); } clock_t t2 = std::clock(); diff --git a/docs/cpque_api.md b/docs/cpque_api.md index cd6b8a9e..7e0021d3 100644 --- a/docs/cpque_api.md +++ b/docs/cpque_api.md @@ -63,13 +63,13 @@ using_cpque(i, cvec_i, >); // adaptor type, '>' = min-heap int main() { size_t N = 10000000; - crand_t rng = crand_init(1234); - crand_uniform_t dist = crand_uniform_init(0, N * 10); + stc64_t rng = stc64_init(1234); + stc64_uniform_t dist = stc64_uniform_init(0, N * 10); cpque_i heap = cpque_i_init(); // Push ten million random numbers to priority queue, plus some negative ones. c_forrange (N) - cpque_i_push(&heap, crand_uniform(&rng, &dist)); + cpque_i_push(&heap, stc64_uniform(&rng, &dist)); c_push_items(&heap, cpque_i, {-231, -32, -873, -4, -343}); // Extract and display the fifty smallest. diff --git a/docs/crand_api.md b/docs/crand_api.md index ddab9b23..3299b046 100644 --- a/docs/crand_api.md +++ b/docs/crand_api.md @@ -1,30 +1,29 @@ # Module [crand](../stc/crand.h): Pseudo Random Number Generators -This describes the API of module **crand**. It contains a *64-bit PRNG*, and can generate +This describes the API of module **crand**. It contains **stc64**, 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. +**stc64** 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 *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 +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 (state[3]), 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. -**crand** passes *PractRand*, tested up to 8TB output, Vigna's Hamming weight test, and simple +**stc64** 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 | +| `stc64_t` | `struct {uint64_t state[4];}` | The PRNG engine type | +| `stc64_uniform_t` | `struct {int64_t lower; uint64_t range;}` | Integer uniform distribution | +| `stc64_uniformf_t` | `struct {double lower, range;}` | Real number uniform distr. | +| `stc64_normalf_t` | `struct {double mean, stddev;}` | Normal distribution type | ## Header file @@ -36,16 +35,16 @@ All cstr definitions and prototypes may be included in your C source file by inc ## 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) stc64_t stc64_init(uint64_t seed); + 2) stc64_t stc64_with_seq(uint64_t seed, uint64_t seq); + 3) uint64_t stc64_rand(stc64_t* rng); + 4) double stc64_randf(stc64_t* rng); + 5) stc64_uniform_t stc64_uniform_init(int64_t low, int64_t high); + 6) int64_t stc64_uniform(stc64_t* rng, stc64_uniform_t* dist); + 7) stc64_uniformf_t stc64_uniformf_init(double low, double high); + 8) double stc64_uniformf(stc64_t* rng, stc64_uniformf_t* dist); + 9) stc64_normalf_t stc64_normalf_init(double mean, double stddev); +10) double stc64_normalf(stc64_t* rng, stc64_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*]. @@ -82,13 +81,13 @@ int main() // 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); + stc64_t rng = stc64_init(seed); + stc64_normalf_t dist = stc64_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) ); + int index = (int) round( stc64_normalf(&rng, &dist) ); cmap_i_emplace(&mhist, index, 0).first->second += 1; } diff --git a/examples/birthday.c b/examples/birthday.c index a2856a3f..dc94fa53 100644 --- a/examples/birthday.c +++ b/examples/birthday.c @@ -15,11 +15,11 @@ static void test_repeats(void) const static uint64_t mask = (1ull << BITS) - 1; printf("birthday paradox: value range: 2^%d, testing repeats of 2^%d values\n", BITS, BITS_TEST); - crand_t rng = crand_init(seed); + stc64_t rng = stc64_init(seed); cmap_ic m = cmap_ic_init(); cmap_ic_reserve(&m, N); c_forrange (i, N) { - uint64_t k = crand_next(&rng) & mask; + uint64_t k = stc64_rand(&rng) & mask; int v = ++cmap_ic_emplace(&m, k, 0).first->second; if (v > 1) printf("repeated value %llx (%d) at 2^%d\n", k, v, (int) log2(i)); } @@ -32,12 +32,12 @@ void test_distribution(void) { enum {BITS = 26}; printf("distribution test: 2^%d values\n", BITS); - crand_t rng = crand_init(seed); + stc64_t rng = stc64_init(seed); const size_t N = 1ull << BITS ; cmap_x map = cmap_x_init(); c_forrange (N) { - uint64_t k = crand_next(&rng); + uint64_t k = stc64_rand(&rng); ++cmap_x_emplace(&map, k & 0xf, 0).first->second; } diff --git a/examples/ex_gaussian.c b/examples/ex_gaussian.c index 5357a2f8..a5137614 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_t rng = crand_init(seed); - crand_normalf_t dist = crand_normalf_init(Mean, StdDev); + stc64_t rng = stc64_init(seed); + stc64_normalf_t dist = stc64_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) ); + int index = (int) round( stc64_normalf(&rng, &dist) ); cmap_i_emplace(&mhist, index, 0).first->second += 1; } diff --git a/examples/list.c b/examples/list.c index d65da22b..9d6e8d89 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_t rng = crand_init(1234); - crand_uniformf_t dist = crand_uniformf_init(100.0f, n); + stc64_t rng = stc64_init(1234); + stc64_uniformf_t dist = stc64_uniformf_init(100.0f, n); int m = 0; c_forrange (i, int, n) - clist_fx_push_back(&list, crand_uniformf(&rng, &dist)), ++m; + clist_fx_push_back(&list, stc64_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 4eb762ec..4eceb29f 100644 --- a/examples/priority.c +++ b/examples/priority.c @@ -11,19 +11,19 @@ using_cpque(i, cvec_i, >); // min-heap (increasing values) int main() { size_t N = 10000000; - crand_t rng = crand_init(time(NULL)); - crand_uniform_t dist = crand_uniform_init(0, N * 10); + stc64_t rng = stc64_init(time(NULL)); + stc64_uniform_t dist = stc64_uniform_init(0, N * 10); cpque_i heap = cpque_i_init(); // Push ten million random numbers to priority queue c_forrange (N) - cpque_i_push(&heap, crand_uniform(&rng, &dist)); + cpque_i_push(&heap, stc64_uniform(&rng, &dist)); // push some negative numbers too. c_push_items(&heap, cpque_i, {-231, -32, -873, -4, -343}); c_forrange (N) - cpque_i_push(&heap, crand_uniform(&rng, &dist)); + cpque_i_push(&heap, stc64_uniform(&rng, &dist)); // Extract the hundred smallest. diff --git a/examples/queue.c b/examples/queue.c index d1eae6bc..bf287e7c 100644 --- a/examples/queue.c +++ b/examples/queue.c @@ -13,20 +13,20 @@ using_cqueue(i, cdeq_i); int main() { int n = 100000000; - crand_uniform_t dist; - crand_t rng = crand_init(1234); - dist = crand_uniform_init(0, n); + stc64_uniform_t dist; + stc64_t rng = stc64_init(1234); + dist = stc64_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(&rng, &dist)); + cqueue_i_push(&queue, stc64_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(&rng, &dist); + int r = stc64_uniform(&rng, &dist); if (r & 1) ++n, cqueue_i_push(&queue, r); else diff --git a/examples/random.c b/examples/random.c index 94a68607..7325e023 100644 --- a/examples/random.c +++ b/examples/random.c @@ -9,16 +9,16 @@ int main() enum {R = 30}; const size_t N = 1000000000; uint64_t seed = 1234; // time(NULL); - crand_t rng = crand_init(seed); + stc64_t rng = stc64_init(seed); uint64_t sum = 0; - crand_normalf_t dist2 = crand_normalf_init(R / 2.0, R / 6.0); + stc64_normalf_t dist2 = stc64_normalf_init(R / 2.0, R / 6.0); size_t N2 = 10000000; int hist[R] = {0}; sum = 0; c_forrange (N2) { - int n = round((crand_normalf(&rng, &dist2) + 0.5)); + int n = round((stc64_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 += crand_next(&rng); + sum += stc64_rand(&rng); } diff = clock() - before; printf("random : %f secs, %zu %f\n", (float) diff / CLOCKS_PER_SEC, N, (double) sum / N); - crand_uniform_t dist1 = crand_uniform_init(0, 1000); + stc64_uniform_t dist1 = stc64_uniform_init(0, 1000); sum = 0; before = clock(); c_forrange (N) { - sum += crand_uniform(&rng, &dist1); + sum += stc64_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/stc/clist.h b/stc/clist.h index a4636a29..a7561a00 100644 --- a/stc/clist.h +++ b/stc/clist.h @@ -36,10 +36,10 @@ int main() { clist_ix list = clist_inits; - crand_t rng = crand_init(12345); + stc64_t rng = stc64_init(12345); int n; for (int i=0; i<1000000; ++i) // one million - clist_ix_push_back(&list, crand_next(&rng) >> 32); + clist_ix_push_back(&list, stc64_rand(&rng) >> 32); n = 0; c_foreach (i, clist_ix, list) if (++n % 10000 == 0) printf("%8d: %10zd\n", n, i.ref->value); diff --git a/stc/cpque.h b/stc/cpque.h index 394b85c3..6f3dc201 100644 --- a/stc/cpque.h +++ b/stc/cpque.h @@ -29,13 +29,13 @@ using_cpque(f, cvec_f, >); // min-heap (increasing values) int main() { - crand_t rng = crand_init(1234); - crand_uniformf_t dist = crand_uniformf_init(10.0f, 100.0f); + stc64_t rng = stc64_init(1234); + stc64_uniformf_t dist = stc64_uniformf_init(10.0f, 100.0f); cpque_f queue = cpque_f_init(); // Push ten million random numbers onto the queue. for (int i=0; i<10000000; ++i) - cpque_f_push(&queue, crand_uniformf(&rng, dist)); + cpque_f_push(&queue, stc64_uniformf(&rng, dist)); // Extract the 100 smallest. for (int i=0; i<100; ++i) { printf("%f ", *cpque_f_top(queue)); diff --git a/stc/cqueue.h b/stc/cqueue.h index f16e1aeb..cb4b0563 100644 --- a/stc/cqueue.h +++ b/stc/cqueue.h @@ -32,19 +32,19 @@ int main() { int n = 10000000; - crand_t rng = crand_init(1234); - crand_uniform_t dist = crand_uniform_init(rng, 0, n); + stc64_t rng = stc64_init(1234); + stc64_uniform_t dist = stc64_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 = crand_uniform(&dist); + int r = stc64_uniform(&dist); if (r & 1) ++n, cqueue_i_push(&queue, r); else diff --git a/stc/crand.h b/stc/crand.h index ad477ea4..d58c0681 100644 --- a/stc/crand.h +++ b/stc/crand.h @@ -28,32 +28,32 @@ #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); + stc64_t rng = stc64_init(seed); + stc64_uniform_t dist1 = stc64_uniform_init(1, 6); + stc64_uniformf_t dist2 = stc64_uniformf_init(1.0, 10.0); + stc64_normalf_t dist3 = stc64_normalf_init(1.0, 10.0); + + uint64_t i = stc64_rand(&rng); + int64_t iu = stc64_uniform(&rng, &dist1); + double xu = stc64_uniformf(&rng, &dist2); + double xn = stc64_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; +typedef struct {uint64_t state[4];} stc64_t; +typedef struct {int64_t lower; uint64_t range, threshold;} stc64_uniform_t; +typedef struct {double lower, range;} stc64_uniformf_t; +typedef struct {double mean, stddev, next; bool has_next;} stc64_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 stc64_t stc64_init(uint64_t seed); +STC_API stc64_t stc64_with_seq(uint64_t seed, uint64_t seq); -STC_INLINE uint64_t crand_next(crand_t* rng) { +STC_INLINE uint64_t stc64_rand(stc64_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); @@ -63,20 +63,20 @@ STC_INLINE 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)}; +STC_INLINE double stc64_randf(stc64_t* rng) { + union {uint64_t i; double f;} u = {0x3FF0000000000000ull | (stc64_rand(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 stc64_uniform_t stc64_uniform_init(int64_t low, int64_t high); /* 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 stc64_uniformf_t stc64_uniformf_init(double low, double high) { + stc64_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; +STC_INLINE double stc64_uniformf(stc64_t* rng, stc64_uniformf_t* dist) { + return stc64_randf(rng)*dist->range + dist->lower; } #if defined(__SIZEOF_INT128__) @@ -92,17 +92,17 @@ STC_INLINE double crand_uniformf(crand_t* rng, crand_uniformf_t* dist) { : [lhs] "0" (a), [rhs] "rm" (b)) #endif -STC_INLINE int64_t crand_uniform(crand_t* rng, crand_uniform_t* d) { +STC_INLINE int64_t stc64_uniform(stc64_t* rng, stc64_uniform_t* d) { uint64_t lo, hi; - do { cmul128(crand_next(rng), d->range, &lo, &hi); } while (lo < d->threshold); + do { cmul128(stc64_rand(rng), d->range, &lo, &hi); } while (lo < d->threshold); return d->lower + hi; } /* 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_INLINE stc64_normalf_t stc64_normalf_init(double mean, double stddev) { + stc64_normalf_t dist = {mean, stddev, 0.0, false}; return dist; } -STC_API double crand_normalf(crand_t* rng, crand_normalf_t* dist); +STC_API double stc64_normalf(stc64_t* rng, stc64_normalf_t* dist); #if !defined(STC_HEADER) || defined(STC_IMPLEMENTATION) @@ -117,32 +117,32 @@ STC_API double crand_normalf(crand_t* rng, crand_normalf_t* dist); * 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 stc64_t stc64_init(uint64_t seed) { + return stc64_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); +STC_DEF stc64_t stc64_with_seq(uint64_t seed, uint64_t seq) { + stc64_t rng = {{seed, seed, seed, (seq << 1u) | 1u}}; + for (int i = 0; i < 8; ++i) stc64_rand(&rng); return rng; } /* 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)}; +STC_DEF stc64_uniform_t stc64_uniform_init(int64_t low, int64_t high) { + stc64_uniform_t dist = {low, (uint64_t) (high - low + 1)}; dist.threshold = (uint64_t)(-dist.range) % dist.range; return dist; } /* Marsaglia polar method for gaussian/normal distribution. */ -STC_DEF double crand_normalf(crand_t* rng, crand_normalf_t* dist) { +STC_DEF double stc64_normalf(stc64_t* rng, stc64_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; + u1 = 2.0 * stc64_randf(rng) - 1.0; + u2 = 2.0 * stc64_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