From eb748a6f7778237e77ca1704fd857e3570c1e9c4 Mon Sep 17 00:00:00 2001 From: Tyge Løvset Date: Sat, 1 Aug 2020 23:55:45 +0200 Subject: Renamed files cvec_pq.h --> cpqueue.h and crand.h --> crandom.h --- README.md | 12 ++-- examples/benchmark.c | 8 +-- examples/geek7.c | 2 +- examples/heap.c | 12 ++-- examples/inits.c | 2 +- examples/list.c | 8 +-- examples/priority.c | 10 ++-- examples/rngbirthday.c | 12 ++-- examples/rngtest.c | 24 ++++---- stc/clist.h | 6 +- stc/cpqueue.h | 125 +++++++++++++++++++++++++++++++++++++++++ stc/crand.h | 150 ------------------------------------------------- stc/crandom.h | 150 +++++++++++++++++++++++++++++++++++++++++++++++++ stc/cvec_pq.h | 125 ----------------------------------------- 14 files changed, 323 insertions(+), 323 deletions(-) create mode 100644 stc/cpqueue.h delete mode 100644 stc/crand.h create mode 100644 stc/crandom.h delete mode 100644 stc/cvec_pq.h diff --git a/README.md b/README.md index 52ce7e99..26149a44 100644 --- a/README.md +++ b/README.md @@ -12,9 +12,9 @@ An elegant, fully typesafe, generic, customizable, user-friendly, consistent, an - **stc/cset.h** - A generic **unordered set** implemented in tandem with *unordered map* - **stc/cstr.h** - Compact and powerful **string** class. - **stc/cvec.h** - Dynamic generic **vector** class, works well as a **stack**. -- **stc/cvec_pq.h** - Priority queue adapter for **cvec.h**, as a **heap**. +- **stc/cpqueue.h** - Priority queue adapter for **cvec.h**, as a **heap**. - **stc/copt.h** - Implementation of a **getopt_long()**-like function, *copt_get()*, to parse command line arguments. -- **stc/crand.h** - A few very efficent modern random number generators *pcg32* and my own *64-bit PRNG* inspired by *sfc64*. +- **stc/crandom.h** - A few very efficent modern random number generators *pcg32* and my own *64-bit PRNG* inspired by *sfc64*. - **stc/cdefs.h** - A common include file with some general definitions. The usage of the containers is vert similar to the C++ standard containers, so it should be easy if you are familiar with them. @@ -254,17 +254,17 @@ int main() { #include #include #include -#include +#include declare_clist(fx, double); int main() { clist_fx list = clist_init; - crand_eng64_t eng = crand_eng64_init(time(NULL)); - crand_uniform_f64_t dist = crand_uniform_f64_init(100.0, 1000.0); + crandom_eng64_t eng = crandom_eng64_init(time(NULL)); + crandom_uniform_f64_t dist = crandom_uniform_f64_init(100.0, 1000.0); int k; for (int i = 0; i < 10000000; ++i) - clist_fx_push_back(&list, crand_uniform_f64(&eng, dist)); + clist_fx_push_back(&list, crandom_uniform_f64(&eng, dist)); k = 0; c_foreach (i, clist_fx, list) if (++k <= 100) printf("%8d: %10f\n", k, i.item->value); else break; diff --git a/examples/benchmark.c b/examples/benchmark.c index 07219ec2..7b707260 100644 --- a/examples/benchmark.c +++ b/examples/benchmark.c @@ -1,4 +1,4 @@ -#include +#include #include #include #include "others/khash.h" @@ -26,9 +26,9 @@ KHASH_MAP_INIT_INT64(ii, uint64_t) size_t seed; static const float max_load_factor = 0.77f; -crand_eng64_t rng; -#define SEED(s) rng = crand_eng64_init(seed) -#define RAND(N) (crand_gen_i64(&rng) & ((1 << N) - 1)) +crandom_eng64_t rng; +#define SEED(s) rng = crandom_eng64_init(seed) +#define RAND(N) (crandom_gen_i64(&rng) & ((1 << N) - 1)) #define CMAP_SETUP(tag, Key, Value) cmap_##tag map = cmap_init \ diff --git a/examples/geek7.c b/examples/geek7.c index 58fadcae..9a60c0a9 100644 --- a/examples/geek7.c +++ b/examples/geek7.c @@ -24,7 +24,7 @@ After inserting all the elements excluding the ones which are to be deleted, Pop #include #include #include -#include +#include declare_cmap(ii, int, int); declare_cvec(i, int); diff --git a/examples/heap.c b/examples/heap.c index 72232d58..58895115 100644 --- a/examples/heap.c +++ b/examples/heap.c @@ -1,7 +1,7 @@ #include #include -#include -#include +#include +#include declare_cvec(f, float); declare_cvec_pqueue(f, >); @@ -9,12 +9,12 @@ declare_cvec_pqueue(f, >); int main() { uint32_t seed = time(NULL); - crand_eng32_t pcg = crand_eng32_init(seed); + crandom_eng32_t pcg = crandom_eng32_init(seed); int N = 30000000, M = 100; cvec_f vec = cvec_init; clock_t start = clock(); for (int i=0; i #include #include -#include +#include #include declare_cmap(id, int, cstr_t, cstr_destroy); // Map of int -> cstr_t diff --git a/examples/list.c b/examples/list.c index bf777611..b981b1ec 100644 --- a/examples/list.c +++ b/examples/list.c @@ -1,17 +1,17 @@ #include #include #include -#include +#include declare_clist(fx, double); int main() { int k, n = 100000; clist_fx list = clist_init; - crand_eng64_t eng = crand_eng64_init(time(NULL)); - crand_uniform_f64_t dist = crand_uniform_f64_init(0.0f, n); + crandom_eng64_t eng = crandom_eng64_init(time(NULL)); + crandom_uniform_f64_t dist = crandom_uniform_f64_init(0.0f, n); for (int i = 0; i < 100000; ++i) - clist_fx_push_back(&list, crand_uniform_f64(&eng, dist)); + clist_fx_push_back(&list, crandom_uniform_f64(&eng, dist)); k = 0; c_foreach (i, clist_fx, list) if (++k <= 10) printf("%8d: %10f\n", k, i.item->value); else break; diff --git a/examples/priority.c b/examples/priority.c index 015f24fb..70d89f2a 100644 --- a/examples/priority.c +++ b/examples/priority.c @@ -1,21 +1,21 @@ #include #include -#include +#include #include -#include +#include declare_cvec(i, uint32_t); declare_cvec_pqueue(i, >); // min-heap (increasing values) int main() { - crand_eng32_t pcg = crand_eng32_init(time(NULL)); - crand_uniform_i32_t dist = crand_uniform_i32_init(0, 100000000); + crandom_eng32_t pcg = crandom_eng32_init(time(NULL)); + crandom_uniform_i32_t dist = crandom_uniform_i32_init(0, 100000000); cvec_i heap = cvec_init; // Push ten million random numbers to priority queue for (int i=0; i<10000000; ++i) - cvec_i_pqueue_push(&heap, crand_uniform_i32(&pcg, dist)); + cvec_i_pqueue_push(&heap, crandom_uniform_i32(&pcg, dist)); // Extract the hundred smallest. for (int i=0; i<100; ++i) { diff --git a/examples/rngbirthday.c b/examples/rngbirthday.c index 03c5b783..216a4cb1 100644 --- a/examples/rngbirthday.c +++ b/examples/rngbirthday.c @@ -2,7 +2,7 @@ #include #include -#include +#include #include #include #include @@ -15,12 +15,12 @@ const static uint64_t mask = (1ull << 52) - 1; void repeats(void) { - crand_eng64_t rng = crand_eng64_init(seed); + crandom_eng64_t rng = crandom_eng64_init(seed); cmap_ic m = cmap_init; cmap_ic_reserve(&m, N); clock_t now = clock(); for (size_t i = 0; i < N; ++i) { - uint64_t k = crand_gen_i64(&rng) & mask; + uint64_t k = crandom_gen_i64(&rng) & mask; int v = ++cmap_ic_insert(&m, k, 0)->value; if (v > 1) printf("%zu: %llx - %d\n", i, k, v); } @@ -34,13 +34,13 @@ declare_cvec(x, uint64_t); void distribution(void) { - crand_eng32_t rng = crand_eng32_init(seed); // time(NULL), time(NULL)); + crandom_eng32_t rng = crandom_eng32_init(seed); // time(NULL), time(NULL)); const size_t N = 1ull << 28, M = 1ull << 9; // 1ull << 10; cmap_x map = cmap_x_make(M); clock_t now = clock(); - crand_uniform_i32_t dist = crand_uniform_i32_init(0, M); + crandom_uniform_i32_t dist = crandom_uniform_i32_init(0, M); for (size_t i = 0; i < N; ++i) { - ++cmap_x_insert(&map, crand_uniform_i32(&rng, dist), 0)->value; + ++cmap_x_insert(&map, crandom_uniform_i32(&rng, dist), 0)->value; } float diff = (float) (clock() - now) / CLOCKS_PER_SEC; diff --git a/examples/rngtest.c b/examples/rngtest.c index f8fe44cf..ac464afd 100644 --- a/examples/rngtest.c +++ b/examples/rngtest.c @@ -1,6 +1,6 @@ #include #include -#include +#include #ifdef __cplusplus #include #endif @@ -14,34 +14,34 @@ int main(void) uint64_t v; printf("start\n"); - crand_eng32_t pcg = crand_eng32_init(time(NULL)); + crandom_eng32_t pcg = crandom_eng32_init(time(NULL)); before = clock(); \ v = 0; for (size_t i=0; i #include - #include + #include declare_clist(ix, int64_t); int main() { clist_ix list = clist_init; - crand_eng32_t pcg = crand_eng32_init(12345); + crandom_eng32_t pcg = crandom_eng32_init(12345); int n; for (int i=0; i<1000000; ++i) // one million - clist_ix_push_back(&list, crand_gen_i32(&pcg)); + clist_ix_push_back(&list, crandom_gen_i32(&pcg)); n = 0; c_foreach (i, clist_ix, list) if (++n % 10000 == 0) printf("%8d: %10zd\n", n, i.item->value); diff --git a/stc/cpqueue.h b/stc/cpqueue.h new file mode 100644 index 00000000..771c7b87 --- /dev/null +++ b/stc/cpqueue.h @@ -0,0 +1,125 @@ +/* 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. + */ + +/* Priority Queue using cvec as heap. + + #include + #include + declare_cvec(f, float); + declare_cvec_pqueue(f, >); // min-heap (increasing values) + + int main() { + crandom_eng32_t gen = crandom_eng32_init(1234); + crandom_uniform_f32_t dist = crandom_uniform_f32_init(10.0f, 100.0f); + + cvec_f queue = cvec_init; + // Push ten million random numbers onto the queue. + for (int i=0; i<10000000; ++i) + cvec_f_pqueue_push(&queue, crandom_uniform_f32(&gen, dist)); + // Extract the 100 smallest. + for (int i=0; i<100; ++i) { + printf("%f ", cvec_f_pqueue_top(&queue)); + cvec_f_pqueue_pop(&queue); + } + cvec_f_destroy(&queue); + } +*/ + +#ifndef CPQUEUE__H__ +#define CPQUEUE__H__ + +#include "cvec.h" + +#define declare_cvec_pqueue(tag, cmpOpr) /* < or > */ \ + \ +STC_API void \ +cvec_##tag##_pqueue_build(cvec_##tag* self); \ +STC_API void \ +cvec_##tag##_pqueue_erase(cvec_##tag* self, size_t i); \ +STC_INLINE cvec_##tag##_value_t \ +cvec_##tag##_pqueue_top(cvec_##tag* self) {return self->data[0];} \ +STC_INLINE void \ +cvec_##tag##_pqueue_pop(cvec_##tag* self) {cvec_##tag##_pqueue_erase(self, 0);} \ +STC_API void \ +cvec_##tag##_pqueue_push(cvec_##tag* self, cvec_##tag##_value_t value); \ +STC_API void \ +cvec_##tag##_pqueue_push_n(cvec_##tag *self, const cvec_##tag##_value_t in[], size_t size); \ + \ +implement_cvec_pqueue(tag, cmpOpr) \ +typedef cvec_##tag##_value_t cvec_##tag##_pqueue_input_t + +/* -------------------------- IMPLEMENTATION ------------------------- */ + +#if !defined(STC_HEADER) || defined(STC_IMPLEMENTATION) +#define implement_cvec_pqueue(tag, cmpOpr) \ + \ +STC_INLINE void \ +_cvec_##tag##_pqueue_sift_down(cvec_##tag##_value_t* arr, size_t i, size_t n) { \ + size_t r = i, c = i << 1; \ + while (c <= n) { \ + if (c < n && cvec_##tag##_sort_compare(&arr[c], &arr[c + 1]) cmpOpr 0) \ + ++c; \ + if (cvec_##tag##_sort_compare(&arr[r], &arr[c]) cmpOpr 0) { \ + cvec_##tag##_value_t t = arr[r]; arr[r] = arr[c]; arr[r = c] = t; \ + } else \ + return; \ + c <<= 1; \ + } \ +} \ + \ +STC_API void \ +cvec_##tag##_pqueue_build(cvec_##tag* self) { \ + size_t n = cvec_size(*self); \ + cvec_##tag##_value_t *arr = self->data - 1; \ + for (size_t k = n >> 1; k != 0; --k) \ + _cvec_##tag##_pqueue_sift_down(arr, k, n); \ +} \ + \ +STC_API void \ +cvec_##tag##_pqueue_erase(cvec_##tag* self, size_t i) { \ + size_t n = cvec_size(*self) - 1; \ + self->data[i] = self->data[n]; \ + cvec_##tag##_pop_back(self); \ + _cvec_##tag##_pqueue_sift_down(self->data - 1, i + 1, n); \ +} \ + \ +STC_API void \ +cvec_##tag##_pqueue_push(cvec_##tag* self, cvec_##tag##_value_t value) { \ + cvec_##tag##_push_back(self, value); /* sift-up the value */ \ + size_t n = cvec_size(*self), c = n; \ + cvec_##tag##_value_t *arr = self->data - 1; \ + for (; c > 1 && cvec_##tag##_sort_compare(&arr[c >> 1], &value) cmpOpr 0; c >>= 1) \ + arr[c] = arr[c >> 1]; \ + if (c != n) arr[c] = value; \ +} \ +STC_API void \ +cvec_##tag##_pqueue_push_n(cvec_##tag *self, const cvec_##tag##_value_t in[], size_t size) { \ + cvec_##tag##_reserve(self, cvec_size(*self) + size); \ + for (size_t i=0; i -/* - crand_eng32_t eng = crand_eng32_init(seed); - crand_uniform_f32_t fdist = crand_uniform_f32_init(1.0f, 6.0f); - crand_uniform_i32_t idist = crand_uniform_i32_init(1, 6); - - uint32_t i = crand_gen_i32(&eng); - int j = crand_uniform_i32(&eng, idist); - float r = crand_uniform_f32(&eng, fdist); -*/ - -typedef struct {uint64_t state[2];} crand_eng32_t; -typedef struct {int32_t min, range;} crand_uniform_i32_t; -typedef struct {float min, range;} crand_uniform_f32_t; - -/* 32 bit random number generator engine */ -STC_API crand_eng32_t crand_eng32_with_seq(uint64_t seed, uint64_t seq); -STC_INLINE crand_eng32_t crand_eng32_init(uint64_t seed) { - return crand_eng32_with_seq(seed, seed); -} - -/* int random number generator, range [0, 2^32) */ -STC_API uint32_t crand_gen_i32(crand_eng32_t* rng); - -STC_INLINE float crand_gen_f32(crand_eng32_t* rng) { - union {uint32_t i; float f;} u = {0x3F800000u | (crand_gen_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, high - low + 1}; return dist; -} -STC_INLINE int32_t crand_uniform_i32(crand_eng32_t* rng, crand_uniform_i32_t dist) { - return dist.min + (int32_t) (((uint64_t) crand_gen_i32(rng) * dist.range) >> 32); -} - -/* 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_eng32_t* rng, crand_uniform_f32_t dist) { - return dist.min + crand_gen_f32(rng) * dist.range; -} - - -typedef struct {uint64_t state[4];} crand_eng64_t; -typedef struct {double min, range;} crand_uniform_f64_t; - -/* 64 bit random number generator engine */ -STC_API crand_eng64_t crand_eng64_with_seq(uint64_t seed, uint64_t seq); -STC_INLINE crand_eng64_t crand_eng64_init(uint64_t seed) { - return crand_eng64_with_seq(seed, 1); -} -/* int random number generator, range [0, 2^64) */ -STC_API uint64_t crand_gen_i64(crand_eng64_t* rng); - -STC_INLINE double crand_gen_f64(crand_eng64_t* rng) { - union {uint64_t i; double f;} u = {0x3FF0000000000000ull | (crand_gen_i64(rng) >> 12)}; - return u.f - 1.0; -} - -/* double random number in range [low, high). 52 bit resolution. */ -STC_INLINE crand_uniform_f64_t crand_uniform_f64_init(float low, float high) { - crand_uniform_f64_t dist = {low, high - low}; return dist; -} -STC_INLINE double crand_uniform_f64(crand_eng64_t* rng, crand_uniform_f64_t dist) { - return dist.min + crand_gen_f64(rng) * dist.range; -} - - -#if !defined(STC_HEADER) || defined(STC_IMPLEMENTATION) - -/* PCG32 random number generator: https://www.pcg-random.org/download.html */ - -STC_API crand_eng32_t crand_eng32_with_seq(uint64_t seed, uint64_t seq) { - crand_eng32_t rng = {0u, (seq << 1u) | 1u}; /* inc must be odd */ - crand_gen_i32(&rng); - rng.state[0] += seed; - crand_gen_i32(&rng); - return rng; -} - -STC_API uint32_t crand_gen_i32(crand_eng32_t* rng) { - uint64_t old = rng->state[0]; - rng->state[0] = old * 6364136223846793005ull + rng->state[1]; - uint32_t xors = ((old >> 18u) ^ old) >> 27u; - uint32_t rot = old >> 59u; - return (xors >> rot) | (xors << ((-rot) & 31)); -} - -/* New 64bit PRNG losely based on SFC64: Copyright 2020, Tyge Løvset */ -/* Faster: Updates only 192bit state. Parallel: Ensures unique sequence for each seq (2^63 seqs) */ -/* Minimal period is 2^64 per seq, average ~ 2^127 per seq */ -STC_API crand_eng64_t crand_eng64_with_seq(uint64_t seed, uint64_t seq) { - crand_eng64_t rng = {seed, seed, seed, (seq << 1u) | 1u}; /* increment must be odd */ - for (int i = 0; i < 12; ++i) crand_gen_i64(&rng); - return rng; -} - -STC_API uint64_t crand_gen_i64(crand_eng64_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; -} - -/* // SFC64 random number generator: http://pracrand.sourceforge.net -STC_API uint64_t crand_gen_i64(crand_eng64_t* rng) { - enum {LROT = 24, RSHIFT = 11, LSHIFT = 3}; - uint64_t *s = rng->state; - const uint64_t result = s[0] + s[1] + s[3]++; - s[0] = s[1] ^ (s[1] >> RSHIFT); - s[1] = s[2] + (s[2] << LSHIFT); - s[2] = ((s[2] << LROT) | (s[2] >> (64 - LROT))) + result; - return result; -} -*/ -#endif - -#endif diff --git a/stc/crandom.h b/stc/crandom.h new file mode 100644 index 00000000..e687ce76 --- /dev/null +++ b/stc/crandom.h @@ -0,0 +1,150 @@ +/* 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__ + +#include "cdefs.h" +#include +/* + crandom_eng32_t eng = crandom_eng32_init(seed); + crandom_uniform_f32_t fdist = crandom_uniform_f32_init(1.0f, 6.0f); + crandom_uniform_i32_t idist = crandom_uniform_i32_init(1, 6); + + uint32_t i = crandom_gen_i32(&eng); + int j = crandom_uniform_i32(&eng, idist); + float r = crandom_uniform_f32(&eng, fdist); +*/ + +typedef struct {uint64_t state[2];} crandom_eng32_t; +typedef struct {int32_t min, range;} crandom_uniform_i32_t; +typedef struct {float min, range;} crandom_uniform_f32_t; + +/* 32 bit random number generator engine */ +STC_API crandom_eng32_t crandom_eng32_with_seq(uint64_t seed, uint64_t seq); +STC_INLINE crandom_eng32_t crandom_eng32_init(uint64_t seed) { + return crandom_eng32_with_seq(seed, seed); +} + +/* int random number generator, range [0, 2^32) */ +STC_API uint32_t crandom_gen_i32(crandom_eng32_t* rng); + +STC_INLINE float crandom_gen_f32(crandom_eng32_t* rng) { + union {uint32_t i; float f;} u = {0x3F800000u | (crandom_gen_i32(rng) >> 9)}; + return u.f - 1.0f; +} + +/* int random number generator in range [low, high] */ +STC_INLINE crandom_uniform_i32_t crandom_uniform_i32_init(int32_t low, int32_t high) { + crandom_uniform_i32_t dist = {low, high - low + 1}; return dist; +} +STC_INLINE int32_t crandom_uniform_i32(crandom_eng32_t* rng, crandom_uniform_i32_t dist) { + return dist.min + (int32_t) (((uint64_t) crandom_gen_i32(rng) * dist.range) >> 32); +} + +/* float random number in range [low, high). Note: 23 bit resolution. */ +STC_INLINE crandom_uniform_f32_t crandom_uniform_f32_init(float low, float high) { + crandom_uniform_f32_t dist = {low, high - low}; return dist; +} +STC_INLINE float crandom_uniform_f32(crandom_eng32_t* rng, crandom_uniform_f32_t dist) { + return dist.min + crandom_gen_f32(rng) * dist.range; +} + + +typedef struct {uint64_t state[4];} crandom_eng64_t; +typedef struct {double min, range;} crandom_uniform_f64_t; + +/* 64 bit random number generator engine */ +STC_API crandom_eng64_t crandom_eng64_with_seq(uint64_t seed, uint64_t seq); +STC_INLINE crandom_eng64_t crandom_eng64_init(uint64_t seed) { + return crandom_eng64_with_seq(seed, 1); +} +/* int random number generator, range [0, 2^64) */ +STC_API uint64_t crandom_gen_i64(crandom_eng64_t* rng); + +STC_INLINE double crandom_gen_f64(crandom_eng64_t* rng) { + union {uint64_t i; double f;} u = {0x3FF0000000000000ull | (crandom_gen_i64(rng) >> 12)}; + return u.f - 1.0; +} + +/* double random number in range [low, high). 52 bit resolution. */ +STC_INLINE crandom_uniform_f64_t crandom_uniform_f64_init(float low, float high) { + crandom_uniform_f64_t dist = {low, high - low}; return dist; +} +STC_INLINE double crandom_uniform_f64(crandom_eng64_t* rng, crandom_uniform_f64_t dist) { + return dist.min + crandom_gen_f64(rng) * dist.range; +} + + +#if !defined(STC_HEADER) || defined(STC_IMPLEMENTATION) + +/* PCG32 random number generator: https://www.pcg-random.org/download.html */ + +STC_API crandom_eng32_t crandom_eng32_with_seq(uint64_t seed, uint64_t seq) { + crandom_eng32_t rng = {0u, (seq << 1u) | 1u}; /* inc must be odd */ + crandom_gen_i32(&rng); + rng.state[0] += seed; + crandom_gen_i32(&rng); + return rng; +} + +STC_API uint32_t crandom_gen_i32(crandom_eng32_t* rng) { + uint64_t old = rng->state[0]; + rng->state[0] = old * 6364136223846793005ull + rng->state[1]; + uint32_t xors = ((old >> 18u) ^ old) >> 27u; + uint32_t rot = old >> 59u; + return (xors >> rot) | (xors << ((-rot) & 31)); +} + +/* New 64bit PRNG losely based on SFC64: Copyright 2020, Tyge Løvset */ +/* Faster: Updates only 192bit state. Parallel: Ensures unique sequence for each seq (2^63 seqs) */ +/* Minimal period is 2^64 per seq, average ~ 2^127 per seq */ +STC_API crandom_eng64_t crandom_eng64_with_seq(uint64_t seed, uint64_t seq) { + crandom_eng64_t rng = {seed, seed, seed, (seq << 1u) | 1u}; /* increment must be odd */ + for (int i = 0; i < 12; ++i) crandom_gen_i64(&rng); + return rng; +} + +STC_API uint64_t crandom_gen_i64(crandom_eng64_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; +} + +/* // SFC64 random number generator: http://pracrand.sourceforge.net +STC_API uint64_t crandom_gen_i64(crandom_eng64_t* rng) { + enum {LROT = 24, RSHIFT = 11, LSHIFT = 3}; + uint64_t *s = rng->state; + const uint64_t result = s[0] + s[1] + s[3]++; + s[0] = s[1] ^ (s[1] >> RSHIFT); + s[1] = s[2] + (s[2] << LSHIFT); + s[2] = ((s[2] << LROT) | (s[2] >> (64 - LROT))) + result; + return result; +} +*/ +#endif + +#endif diff --git a/stc/cvec_pq.h b/stc/cvec_pq.h deleted file mode 100644 index c6a896d3..00000000 --- a/stc/cvec_pq.h +++ /dev/null @@ -1,125 +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. - */ - -/* Priority Queue using cvec as heap. - - #include - #include - declare_cvec(f, float); - declare_cvec_pqueue(f, >); // min-heap (increasing values) - - int main() { - crand_eng32_t gen = crand_eng32_init(1234); - crand_uniform_f32_t dist = crand_uniform_f32_init(10.0f, 100.0f); - - cvec_f queue = cvec_init; - // Push ten million random numbers onto the queue. - for (int i=0; i<10000000; ++i) - cvec_f_pqueue_push(&queue, crand_uniform_f32(&gen, dist)); - // Extract the 100 smallest. - for (int i=0; i<100; ++i) { - printf("%f ", cvec_f_pqueue_top(&queue)); - cvec_f_pqueue_pop(&queue); - } - cvec_f_destroy(&queue); - } -*/ - -#ifndef CVEC_PQ__H__ -#define CVEC_PQ__H__ - -#include "cvec.h" - -#define declare_cvec_pqueue(tag, cmpOpr) /* < or > */ \ - \ -STC_API void \ -cvec_##tag##_pqueue_build(cvec_##tag* self); \ -STC_API void \ -cvec_##tag##_pqueue_erase(cvec_##tag* self, size_t i); \ -STC_INLINE cvec_##tag##_value_t \ -cvec_##tag##_pqueue_top(cvec_##tag* self) {return self->data[0];} \ -STC_INLINE void \ -cvec_##tag##_pqueue_pop(cvec_##tag* self) {cvec_##tag##_pqueue_erase(self, 0);} \ -STC_API void \ -cvec_##tag##_pqueue_push(cvec_##tag* self, cvec_##tag##_value_t value); \ -STC_API void \ -cvec_##tag##_pqueue_push_n(cvec_##tag *self, const cvec_##tag##_value_t in[], size_t size); \ - \ -implement_cvec_pqueue(tag, cmpOpr) \ -typedef cvec_##tag##_value_t cvec_##tag##_pqueue_input_t - -/* -------------------------- IMPLEMENTATION ------------------------- */ - -#if !defined(STC_HEADER) || defined(STC_IMPLEMENTATION) -#define implement_cvec_pqueue(tag, cmpOpr) \ - \ -STC_INLINE void \ -_cvec_##tag##_pqueue_sift_down(cvec_##tag##_value_t* arr, size_t i, size_t n) { \ - size_t r = i, c = i << 1; \ - while (c <= n) { \ - if (c < n && cvec_##tag##_sort_compare(&arr[c], &arr[c + 1]) cmpOpr 0) \ - ++c; \ - if (cvec_##tag##_sort_compare(&arr[r], &arr[c]) cmpOpr 0) { \ - cvec_##tag##_value_t t = arr[r]; arr[r] = arr[c]; arr[r = c] = t; \ - } else \ - return; \ - c <<= 1; \ - } \ -} \ - \ -STC_API void \ -cvec_##tag##_pqueue_build(cvec_##tag* self) { \ - size_t n = cvec_size(*self); \ - cvec_##tag##_value_t *arr = self->data - 1; \ - for (size_t k = n >> 1; k != 0; --k) \ - _cvec_##tag##_pqueue_sift_down(arr, k, n); \ -} \ - \ -STC_API void \ -cvec_##tag##_pqueue_erase(cvec_##tag* self, size_t i) { \ - size_t n = cvec_size(*self) - 1; \ - self->data[i] = self->data[n]; \ - cvec_##tag##_pop_back(self); \ - _cvec_##tag##_pqueue_sift_down(self->data - 1, i + 1, n); \ -} \ - \ -STC_API void \ -cvec_##tag##_pqueue_push(cvec_##tag* self, cvec_##tag##_value_t value) { \ - cvec_##tag##_push_back(self, value); /* sift-up the value */ \ - size_t n = cvec_size(*self), c = n; \ - cvec_##tag##_value_t *arr = self->data - 1; \ - for (; c > 1 && cvec_##tag##_sort_compare(&arr[c >> 1], &value) cmpOpr 0; c >>= 1) \ - arr[c] = arr[c >> 1]; \ - if (c != n) arr[c] = value; \ -} \ -STC_API void \ -cvec_##tag##_pqueue_push_n(cvec_##tag *self, const cvec_##tag##_value_t in[], size_t size) { \ - cvec_##tag##_reserve(self, cvec_size(*self) + size); \ - for (size_t i=0; i