From 66c6274103c6501b74e8ae68df088d5b8842141a Mon Sep 17 00:00:00 2001 From: Tyge Løvset Date: Wed, 29 Jul 2020 22:21:19 +0200 Subject: Changed rand.h API. Using my own 64-bit random function, not sfc64: Faster and has seq (sequence, jump) for multiple threads. Cleanup. --- examples/benchmark.c | 2 +- examples/heap.c | 4 +-- examples/list.c | 29 ++++++++-------- examples/prime.c | 11 +----- examples/priority.c | 2 +- examples/rngbirthday.c | 8 ++--- examples/rngtest.c | 18 ++++++---- stc/clist.h | 2 +- stc/crand.h | 92 ++++++++++++++++++++++++++++---------------------- stc/cvecpq.h | 2 +- 10 files changed, 89 insertions(+), 81 deletions(-) diff --git a/examples/benchmark.c b/examples/benchmark.c index ea6c1794..f8536a16 100644 --- a/examples/benchmark.c +++ b/examples/benchmark.c @@ -22,7 +22,7 @@ size_t seed; static const float max_load_factor = 0.77f; crand_eng64_t rng; -#define SEED(s) rng = crand_eng64(seed) +#define SEED(s) rng = crand_eng64_init(seed) #define RAND(N) (crand_gen_i64(&rng) & ((1 << N) - 1)) diff --git a/examples/heap.c b/examples/heap.c index 55aac409..0f4ab4c2 100644 --- a/examples/heap.c +++ b/examples/heap.c @@ -9,7 +9,7 @@ declare_cvec_priority_queue(f, >); int main() { uint32_t seed = time(NULL); - crand_eng32_t pcg = crand_eng32(seed); + crand_eng32_t pcg = crand_eng32_init(seed); int N = 30000000, M = 100; cvec_f vec = cvec_init; clock_t start = clock(); @@ -25,7 +25,7 @@ int main() cvecpq_f_pop(&vec); printf("\n\npopped PQ: %f secs\n", (clock() - start) / (float) CLOCKS_PER_SEC); - pcg = crand_eng32(seed); + pcg = crand_eng32_init(seed); start = clock(); for (int i=0; i #include #include -declare_clist(ix, uint64_t); +declare_clist(fx, double); int main() { - clist_ix list = clist_init; - crand_eng32_t pcg = crand_eng32(time(NULL)); + clist_fx list = clist_init; + crand_eng64_t eng = crand_eng64_init(time(NULL)); + crand_uniform_f64_t dist = crand_uniform_f64_init(1.0, 100.0); int n; - for (int i=0; i<10000000; ++i) // ten million - clist_ix_push_back(&list, crand_gen_i32(&pcg)); + for (int i = 0; i < 10000000; ++i) // ten million + clist_fx_push_back(&list, crand_uniform_f64(&eng, dist)); n = 100; - c_foreach (i, clist_ix, list) - if (n--) printf("%8d: %10zu\n", 100 - n, i.item->value); else break; + c_foreach (i, clist_fx, list) + if (n--) printf("%8d: %10f\n", 100 - n, i.item->value); else break; // Sort them... - clist_ix_sort(&list); // mergesort O(n*log n) + clist_fx_sort(&list); // mergesort O(n*log n) n = 100; puts("sorted"); - c_foreach (i, clist_ix, list) - if (n--) printf("%8d: %10zu\n", 100 - n, i.item->value); else break; + c_foreach (i, clist_fx, list) + if (n--) printf("%8d: %10f\n", 100 - n, i.item->value); else break; - clist_ix_clear(&list); - c_push(&list, clist_ix, c_items(10, 20, 30, 40, 50)); - c_foreach (i, clist_ix, list) printf("%zu ", i.item->value); + clist_fx_clear(&list); + c_push(&list, clist_fx, c_items(10, 20, 30, 40, 50)); + c_foreach (i, clist_fx, list) printf("%f ", i.item->value); puts(""); - clist_ix_destroy(&list); + clist_fx_destroy(&list); } \ No newline at end of file diff --git a/examples/prime.c b/examples/prime.c index 611690ac..2e6c99ee 100644 --- a/examples/prime.c +++ b/examples/prime.c @@ -1,12 +1,5 @@ -#include - -#if defined(__GNUC__) -#define cbitset_popcnt64(i) __builtin_popcountll(i) -#else -#define cbitset_popcnt64(i) _mm_popcnt_u64(i) -#endif - #include +#include static inline void sieveOfEratosthenes(size_t n) { @@ -15,8 +8,6 @@ static inline void sieveOfEratosthenes(size_t n) cbitset_reset(&prime, 0); cbitset_reset(&prime, 1); - uint64_t m = cbitset_popcnt64(123456); - for (size_t i = 2; i <= n; ++i) { // If prime[i] is not changed, then it is a prime if (cbitset_test(prime, i) && i*i <= n) { diff --git a/examples/priority.c b/examples/priority.c index cf1c386c..64449304 100644 --- a/examples/priority.c +++ b/examples/priority.c @@ -9,7 +9,7 @@ declare_cvec(i, uint32_t); declare_cvec_priority_queue(i, >); // min-heap (increasing values) int main() { - crand_eng32_t pcg = crand_eng32(time(NULL)); + crand_eng32_t pcg = crand_eng32_init(time(NULL)); cvec_i heap = cvec_init; // Push ten million random numbers to queue diff --git a/examples/rngbirthday.c b/examples/rngbirthday.c index 43718fe8..c20bac26 100644 --- a/examples/rngbirthday.c +++ b/examples/rngbirthday.c @@ -15,7 +15,7 @@ const static uint64_t mask = (1ull << 52) - 1; void repeats(void) { - crand_eng64_t rng = crand_eng64(seed); + crand_eng64_t rng = crand_eng64_init(seed); cmap_ic m = cmap_init; cmap_ic_reserve(&m, N); clock_t now = clock(); @@ -34,13 +34,13 @@ declare_cvec(x, uint64_t); void distribution(void) { - crand_eng32_t rng = crand_eng32(seed); // time(NULL), time(NULL)); + crand_eng32_t rng = crand_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_i32_uniform_t dist = crand_i32_uniform(0, M); + crand_uniform_i32_t dist = crand_uniform_i32_init(0, M); for (size_t i = 0; i < N; ++i) { - ++cmap_x_insert(&map, crand_gen_i32_uniform(&rng, dist), 0)->value; + ++cmap_x_insert(&map, crand_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 d33daa3b..f8fe44cf 100644 --- a/examples/rngtest.c +++ b/examples/rngtest.c @@ -14,7 +14,7 @@ int main(void) uint64_t v; printf("start\n"); - crand_eng32_t pcg = crand_eng32(time(NULL)); + crand_eng32_t pcg = crand_eng32_init(time(NULL)); before = clock(); \ v = 0; for (size_t i=0; i /* - crand_eng32_t eng = crand_eng32(seed); - crand_f32_uniform_t fdist = crand_f32_uniform(1.0f, 6.0f); - crand_i32_uniform_t idist = crand_i32_uniform(1, 6); + 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_gen_i32_uniform(&eng, idist); - float r = crand_gen_f32_uniform(&eng, fdist); + int j = crand_uniform_i32(&eng, idist); + float r = crand_uniform_f32(&eng, fdist); */ -typedef struct {uint64_t state, inc;} crand_eng32_t; -typedef struct {int32_t min, range;} crand_i32_uniform_t; -typedef struct {float min, range;} crand_f32_uniform_t; +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_id(uint64_t seed, uint64_t seq); - -STC_INLINE crand_eng32_t crand_eng32(uint64_t seed) { - return crand_eng32_with_id(seed, 0); +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, 1); } /* int random number generator, range [0, 2^32) */ @@ -56,28 +55,30 @@ STC_INLINE float crand_gen_f32(crand_eng32_t* rng) { } /* int random number generator in range [low, high] */ -STC_INLINE crand_i32_uniform_t crand_i32_uniform(int32_t low, int32_t high) { - crand_i32_uniform_t dist = {low, high - low}; return dist; +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 uint32_t crand_gen_i32_uniform(crand_eng32_t* rng, crand_i32_uniform_t 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_f32_uniform_t crand_f32_uniform(float low, float high) { - crand_f32_uniform_t dist = {low, high - low}; return dist; +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_gen_f32_uniform(crand_eng32_t* rng, crand_f32_uniform_t 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_f64_uniform_t; +typedef struct {double min, range;} crand_uniform_f64_t; /* 64 bit random number generator engine */ -STC_API crand_eng64_t crand_eng64(const uint64_t seed); - +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); @@ -87,10 +88,10 @@ STC_INLINE double crand_gen_f64(crand_eng64_t* rng) { } /* double random number in range [low, high). 52 bit resolution. */ -STC_INLINE crand_f64_uniform_t crand_f64_uniform(float low, float high) { - crand_f64_uniform_t dist = {low, high - low}; return dist; +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_gen_f64_uniform(crand_eng64_t* rng, crand_f64_uniform_t 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; } @@ -99,40 +100,49 @@ STC_INLINE double crand_gen_f64_uniform(crand_eng64_t* rng, crand_f64_uniform_t /* PCG32 random number generator: https://www.pcg-random.org/index.html */ -STC_API crand_eng32_t crand_eng32_with_id(uint64_t seed, uint64_t seq) { +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 += seed; + 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; - rng->state = old * 6364136223846793005ull + rng->inc; - uint32_t xos = ((old >> 18u) ^ old) >> 27u; + 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 (xos >> rot) | (xos << ((-rot) & 31)); + return (xors >> rot) | (xors << ((-rot) & 31)); } /* SFC64 random number generator: http://pracrand.sourceforge.net */ -STC_API crand_eng64_t crand_eng64(const uint64_t seed) { - crand_eng64_t state = {{seed, seed, seed, 1}}; - for (int i = 0; i < 12; ++i) crand_gen_i64(&state); - return state; +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 {LR=24, RS=11, LS=3}; +#if USE_SFC64 +STC_API uint64_t crand_gen_i64(crand_eng64_t* rng) { /* original sfc64 */ + 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] >> RS); - s[1] = s[2] + (s[2] << LS); - s[2] = (s[2] << LR) | (s[2] >> (64 - LR)); + 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; } - +#else +STC_API uint64_t crand_gen_i64(crand_eng64_t* rng) { /* copyright: Tyge Løvset */ + 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; +} +#endif #endif #endif diff --git a/stc/cvecpq.h b/stc/cvecpq.h index ba81f51b..2e6f2d09 100644 --- a/stc/cvecpq.h +++ b/stc/cvecpq.h @@ -28,7 +28,7 @@ declare_cvec(i, int); declare_cvec_priority_queue(i, >); // min-heap (increasing values) int main() { - crand_eng32_t pcg = crand_eng32(1234); + crand_eng32_t pcg = crand_eng32_init(1234); cvec_i heap = cvec_init; // Push one million random numbers onto the queue. for (int i=0; i<1000000; ++i) -- cgit v1.2.3