Changed rand.h API. Using my own 64-bit random function, not sfc64: Faster and has seq (sequence, jump) for multiple threads. Cleanup.

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<N; ++i)

         cvecpq_f_push(&vec, crand_gen_i32(&pcg));

diff --git a/examples/list.c b/examples/list.c
index 2689d187..a65a8454 100644
--- a/examples/list.c
+++ b/examples/list.c
@@ -2,28 +2,29 @@
 #include <time.h>

 #include <stc/clist.h>

 #include <stc/crand.h>

-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 <stc/cbitset.h>

-

-#if defined(__GNUC__)

-#define cbitset_popcnt64(i) __builtin_popcountll(i)

-#else

-#define cbitset_popcnt64(i) _mm_popcnt_u64(i)

-#endif

-

 #include <stdio.h>

+#include <stc/cbitset.h>

 

 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<NN; i++) {

@@ -23,7 +23,7 @@ int main(void)
     difference = clock() - before;

     printf("pcg32: %.02f, %zu\n", (float) difference / CLOCKS_PER_SEC, v);

 

-    crand_eng64_t sfc = crand_eng64(time(NULL));

+    crand_eng64_t sfc = crand_eng64_init(time(NULL));

     before = clock(); \

     v = 0;

     for (size_t i=0; i<NN; i++) {

@@ -32,10 +32,16 @@ int main(void)
     difference = clock() - before;

     printf("sfc64: %.02f, %zu\n", (float) difference / CLOCKS_PER_SEC, v);

 

-    crand_f64_uniform_t fdist = crand_f64_uniform(10, 20);

-    for (int i=0; i<8; ++i) printf("%f ", crand_gen_f64_uniform(&sfc, fdist));

+    crand_uniform_i32_t i32dist = crand_uniform_i32_init(10, 20);

+    for (int i=0; i<8; ++i) printf("%d ", crand_uniform_i32(&pcg, i32dist));

     puts("");

-    crand_f32_uniform_t f32dist = crand_f32_uniform(10, 20);

-    for (int i=0; i<8; ++i) printf("%f ", crand_gen_f32_uniform(&pcg, f32dist));

+

+    crand_uniform_f32_t f32dist = crand_uniform_f32_init(10, 20);

+    for (int i=0; i<8; ++i) printf("%f ", crand_uniform_f32(&pcg, f32dist));

+    puts("");

+

+    crand_uniform_f64_t fdist = crand_uniform_f64_init(10, 20);

+    for (int i=0; i<8; ++i) printf("%f ", crand_uniform_f64(&sfc, fdist));

+    //for (int i=0; i<8; ++i) printf("%zu ", crand_gen_i64(&sfc));

     puts("");

 }
\ No newline at end of file
diff --git a/stc/clist.h b/stc/clist.h
index 6f0b2b9c..e55880ba 100644
--- a/stc/clist.h
+++ b/stc/clist.h
@@ -39,7 +39,7 @@
  

     int main() {

         clist_ix list = clist_init;

-        crand_eng32_t pcg = crand_eng32(12345);

+        crand_eng32_t pcg = crand_eng32_init(12345);

         int n;

         for (int i=0; i<1000000; ++i) // one million

             clist_ix_push_back(&list, crand_gen_i32(&pcg));

diff --git a/stc/crand.h b/stc/crand.h
index 105275c6..841cb473 100644
--- a/stc/crand.h
+++ b/stc/crand.h
@@ -27,24 +27,23 @@
 #include "cdefs.h"

 #include <string.h>

 /*

-    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)