Reworked the crandom.h module. Changed API (sorry!), made uniform distr. unbiased. Removed 32-bit pcg32 PRNG.

16 files changed, 182 insertions, 294 deletions

diff --git a/docs/cpqueue_api.md b/docs/cpqueue_api.md
index 39569992..1e92277a 100644
--- a/docs/cpqueue_api.md
+++ b/docs/cpqueue_api.md
@@ -63,13 +63,13 @@ using_cpqueue(i, cvec_i, >); // adaptor type, '>' = min-heap
 int main()
 {
     size_t N = 10000000;
-    crand_rng64_t rng = crand_rng64_init(1234);
-    crand_uniform_i64_t dist = crand_uniform_i64_init(0, N * 10);
+    cstc64_t rng = cstc64_init(1234);
+    cstc64_uniform_t dist = cstc64_uniform_init(0, N * 10);
 
     cpqueue_i heap = cpqueue_i_init();
     // Push ten million random numbers to priority queue, plus some negative ones.
     c_forrange (N)
-        cpqueue_i_push(&heap, crand_uniform_i64(&rng, &dist));
+        cpqueue_i_push(&heap, cstc64_uniform(&rng, &dist));
     c_push_items(&heap, cpqueue_i, {-231, -32, -873, -4, -343});
 
     // Extract and display the fifty smallest.
diff --git a/docs/crandom_api.md b/docs/crandom_api.md
index 9af69b9e..57697067 100644
--- a/docs/crandom_api.md
+++ b/docs/crandom_api.md
@@ -1,19 +1,16 @@
-# Module crand: Pseudo Random Number Generators
+# Module crandom: Pseudo Random Number Generators
 
 This describes the API of module **crand**. It contains *pcg32* created by Melissa O'Neill, and an
 *64-bit PRNG* created by Tyge Løvset. The PRNG can generate uniform and normal distributions.
 
 ## Types
 
-| Name                  | Type definition                             | Used to represent...                 |
-|:----------------------|:--------------------------------------------|:-------------------------------------|
-| `crand_rng32_t`       | `struct {uint64_t state[2];}`               | The crandom type                     |
-| `crand_uniform_i32_t` | `struct {int32_t offset; uint32_t range;}`  | The crandom element type             |
-| `crand_uniform_f32_t` | `struct {float offset, range;}`             | crandom iterator                     |
-| `crand_rng64_t`       | `struct {uint64_t state[4];}`               |                                      |
-| `crand_uniform_i64_t` | `struct {int64_t offset; uint64_t range;}`  |                                      |
-| `crand_uniform_f64_t` | `struct {double offset, range;}`            |                                      |
-| `crand_normal_f64_t`  | `struct {double mean, stddev, ...;}`        |                                      |
+| Name                | Type definition                             | Used to represent...         |
+|:--------------------|:--------------------------------------------|:-----------------------------|
+| `cstc64_t`          | `struct {uint64_t state[4];}`               | The PRNG engine type         |
+| `cstc64_uniform_t`  | `struct {int64_t offset; uint64_t range;}`  | Integer uniform distribution |
+| `cstc64_uniformf_t` | `struct {double offset, range;}`            | Real number uniform distr.   |
+| `cstc64_normalf_t`  | `struct {double mean, stddev;}`             | Normal distribution type     |
 
 ## Header file
 
@@ -25,41 +22,28 @@ All cstr definitions and prototypes may be included in your C source file by inc
 ## Methods
 
 ```c
- 1)     crand_rng32_t           crand_rng32_init(uint64_t seed);
- 2)     crand_rng32_t           crand_rng32_with_seq(uint64_t seed, uint64_t seq);
- 3)     uint32_t                crand_i32(crand_rng32_t* rng);
- 4)     float                   crand_f32(crand_rng32_t* rng);
- 5)     crand_uniform_i32_t     crand_uniform_i32_init(int32_t low, int32_t high);
- 6)     int32_t                 crand_uniform_i32(crand_rng32_t* rng, crand_uniform_i32_t* dist);
- 7)     uint32_t                crand_unbiased_i32(crand_rng32_t* rng, crand_uniform_i32_t* dist);
- 8)     crand_uniform_f32_t     crand_uniform_f32_init(float low, float high); /*  */
- 9)     float                   crand_uniform_f32(crand_rng32_t* rng, crand_uniform_f32_t* dist);
+ 1)     cstc64_t            cstc64_with_seq(uint64_t seed, uint64_t seq);
+ 2)     cstc64_t            cstc64_init(uint64_t seed);
+ 3)     uint64_t            cstc64_rand(cstc64_t* rng);
+ 4)     double              cstc64_randf(cstc64_t* rng);
+ 5)     cstc64_uniform_t    cstc64_uniform_init(int64_t low, int64_t high);
+ 6)     int64_t             cstc64_uniform(cstc64_t* rng, cstc64_uniform_t* dist);
+ 7)     cstc64_uniformf_t   cstc64_uniformf_init(double low, double high);
+ 8)     double              cstc64_uniformf(cstc64_t* rng, cstc64_uniformf_t* dist);
+ 9)     cstc64_normalf_t    cstc64_normalf_init(double mean, double stddev);
+10)     double              cstc64_normalf(cstc64_t* rng, cstc64_normalf_t* dist);
 ```
-`1-2)` PRNG 32-bit engine initializers. `3)` Integer RNG with range \[0, 2^32). `4)` Float RNG with range \[0, 1). 
-`5-6)` Uniform integer RNG with range \[*low*, *high*]. `7)` Unbiased version, see https://github.com/lemire/fastrange.
-`8-9)` Uniform float RNG with range \[*low*, *high*).
-```c
- 1)     crand_rng64_t           crand_rng64_with_seq(uint64_t seed, uint64_t seq);
- 2)     crand_rng64_t           crand_rng64_init(uint64_t seed);
- 3)     uint64_t                crand_i64(crand_rng64_t* rng);
- 4)     double                  crand_f64(crand_rng64_t* rng);
- 5)     crand_uniform_i64_t     crand_uniform_i64_init(int64_t low, int64_t high);
- 6)     int64_t                 crand_uniform_i64(crand_rng64_t* rng, crand_uniform_i64_t* dist);
- 7)     crand_uniform_f64_t     crand_uniform_f64_init(double low, double high);
- 8)     double                  crand_uniform_f64(crand_rng64_t* rng, crand_uniform_f64_t* dist);
- 9)     crand_normal_f64_t      crand_normal_f64_init(double mean, double stddev);
-10)     double                  crand_normal_f64(crand_rng64_t* rng, crand_normal_f64_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 were around 68% of the values are within the range [*mean* - *stddev, *mean* + *stddev*].
-
-The method `crand_i64(crand_rng64_t* rng)` is an extremely fast PRNG suited for parallel usage, featuring
-a Weyl-sequence as part of the state. It is faster than *sfc64*, *wyhash64*, *pcg*, and the *xoroshiro*
-families of RPNGs. It does not require fast multiplication or 128-bit integer operations. The state is
-256-bits, but updates only 192 bit per generated number. It can create 2^63 unique threads with minimum period
-length of 2^64 per thread (expected period lengths 2^127). There is no *jump function*, but incrementing
-the Weyl-increment by 2, achieves the same goal. Passes *PractRand*, tested up to 8TB output.
+`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 are within the range [*mean* - *stddev, *mean* + *stddev*].
+
+The method `cstc64_rand(cstc64_t* rng)` 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 GCC. 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 incrementing
+the Weyl-increment by 2, initializes a new unique 2^64 *minimum* length period. 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.
 
 ## Example
 ```c
@@ -91,13 +75,13 @@ int main()
 
     // Setup random engine with normal distribution.
     uint64_t seed = time(NULL);
-    crand_rng64_t rng = crand_rng64_init(seed);
-    crand_normal_f64_t dist = crand_normal_f64_init(Mean, StdDev);
+    cstc64_t rng = cstc64_init(seed);
+    cstc64_normalf_t dist = cstc64_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_normal_f64(&rng, &dist) );
+        int index = (int) round( cstc64_normalf(&rng, &dist) );
         cmap_i_emplace(&mhist, index, 0).first->second += 1;
     }
 
diff --git a/examples/benchmark.cpp b/examples/benchmark.cpp
index ca744e4d..4830c9ba 100644
--- a/examples/benchmark.cpp
+++ b/examples/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_rng64_t rng;

-#define SEED(s) rng = crand_rng64_init(seed)

-#define RAND(N) (crand_i64(&rng) & ((1 << N) - 1))

+cstc64_t rng;

+#define SEED(s) rng = cstc64_init(seed)

+#define RAND(N) (cstc64_rand(&rng) & ((1 << N) - 1))

 

 

 #define CMAP_SETUP(X, Key, Value) cmap_##X map = cmap_inits \

diff --git a/examples/birthday.c b/examples/birthday.c
index 54ffa30e..fedb2c5c 100644
--- a/examples/birthday.c
+++ b/examples/birthday.c
@@ -14,12 +14,12 @@ const static uint64_t mask = (1ull << 52) - 1;
 

 void repeats(void)

 {

-    crand_rng64_t rng = crand_rng64_init(seed);

+    cstc64_t rng = cstc64_init(seed);

     cmap_ic m = cmap_ic_init();

     cmap_ic_reserve(&m, N);

     clock_t now = clock();

     c_forrange (i, N) {

-        uint64_t k = crand_i64(&rng) & mask;

+        uint64_t k = cstc64_rand(&rng) & mask;

         int v = ++cmap_ic_emplace(&m, k, 0).first->second;

         if (v > 1) printf("%zu: %llx - %d\n", i, k, v);

     }

@@ -33,14 +33,14 @@ using_cvec(x, uint64_t);
 

 void distribution(void)

 {

-    crand_rng32_t rng = crand_rng32_init(seed); // time(NULL), time(NULL));

+    cstc64_t rng = cstc64_init(seed); // time(NULL), time(NULL));

     const size_t N = 1ull << 28, M = 1ull << 9; // 1ull << 10;

     cmap_x map = cmap_x_with_capacity(M);

     clock_t now = clock();

-    crand_uniform_i32_t dist = crand_uniform_i32_init(0, M);

+    cstc64_uniform_t dist = cstc64_uniform_init(0, M);

 

     c_forrange (N) {

-        ++cmap_x_emplace(&map, crand_uniform_i32(&rng, &dist), 0).first->second;

+        ++cmap_x_emplace(&map, cstc64_uniform(&rng, &dist), 0).first->second;

     }

     float diff = (float) (clock() - now) / CLOCKS_PER_SEC;

 

diff --git a/examples/ex_gaussian.c b/examples/ex_gaussian.c
index 3ed30cf2..5303dd0c 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_rng64_t rng = crand_rng64_init(seed);

-    crand_normal_f64_t dist = crand_normal_f64_init(Mean, StdDev);

+    cstc64_t rng = cstc64_init(seed);

+    cstc64_normalf_t dist = cstc64_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_normal_f64(&rng, &dist) );

+        int index = (int) round( cstc64_normalf(&rng, &dist) );

         cmap_i_emplace(&mhist, index, 0).first->second += 1;

     }

 

diff --git a/examples/heap.c b/examples/heap.c
index 873077fe..166dd6ba 100644
--- a/examples/heap.c
+++ b/examples/heap.c
@@ -10,15 +10,15 @@ using_cpqueue(f, cvec_f, >);
 int main()

 {

     uint32_t seed = time(NULL);

-    crand_rng32_t pcg;

+    cstc64_t rng;

     int N = 3000000, M = 100;

 

     cpqueue_f pq = cpqueue_f_init();

 

-    pcg = crand_rng32_init(seed);

+    rng = cstc64_init(seed);

     clock_t start = clock();

     c_forrange (i, int, N)

-        cvec_f_push_back(&pq, (float) crand_f32(&pcg)*100000);

+        cvec_f_push_back(&pq, (float) cstc64_randf(&rng)*100000);

 

     cpqueue_f_make_heap(&pq);

     printf("Built priority queue: %f secs\n", (clock() - start) / (float) CLOCKS_PER_SEC);

@@ -33,10 +33,9 @@ int main()
         cpqueue_f_pop(&pq);

     printf("\n\npopped PQ: %f secs\n", (clock() - start) / (float) CLOCKS_PER_SEC);

 

-    pcg = crand_rng32_init(seed);

     start = clock();

     c_forrange (i, int, N)

-        cpqueue_f_push(&pq, (float) crand_f32(&pcg)*100000);

+        cpqueue_f_push(&pq, (float) cstc64_randf(&rng)*100000);

     printf("pushed PQ: %f secs\n", (clock() - start) / (float) CLOCKS_PER_SEC);

 

     c_forrange (i, int, M) {

diff --git a/examples/list.c b/examples/list.c
index 4758b072..fbfa8830 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_rng64_t eng = crand_rng64_init(1234);

-    crand_uniform_f64_t dist = crand_uniform_f64_init(100.0f, n);

+    cstc64_t eng = cstc64_init(1234);

+    cstc64_uniformf_t dist = cstc64_uniformf_init(100.0f, n);

     int m = 0;

     c_forrange (i, int, n)

-        clist_fx_push_back(&list, crand_uniform_f64(&eng, &dist)), ++m;

+        clist_fx_push_back(&list, cstc64_uniformf(&eng, &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 e4a20da9..e81cd8a3 100644
--- a/examples/priority.c
+++ b/examples/priority.c
@@ -11,19 +11,19 @@ using_cpqueue(i, cvec_i, >); // min-heap (increasing values)
 

 int main() {

     size_t N = 10000000;

-    crand_rng64_t pcg = crand_rng64_init(time(NULL));

-    crand_uniform_i64_t dist = crand_uniform_i64_init(0, N * 10);

+    cstc64_t pcg = cstc64_init(time(NULL));

+    cstc64_uniform_t dist = cstc64_uniform_init(0, N * 10);

     cpqueue_i heap = cpqueue_i_init();

 

     // Push ten million random numbers to priority queue

     c_forrange (N)

-        cpqueue_i_push(&heap, crand_uniform_i64(&pcg, &dist));

+        cpqueue_i_push(&heap, cstc64_uniform(&pcg, &dist));

 

     // push some negative numbers too.

     c_push_items(&heap, cpqueue_i, {-231, -32, -873, -4, -343});

 

     c_forrange (N)

-        cpqueue_i_push(&heap, crand_uniform_i64(&pcg, &dist));

+        cpqueue_i_push(&heap, cstc64_uniform(&pcg, &dist));

 

 

     // Extract the hundred smallest.

diff --git a/examples/queue.c b/examples/queue.c
index c56214a7..0ff29736 100644
--- a/examples/queue.c
+++ b/examples/queue.c
@@ -7,20 +7,20 @@ using_cqueue(i, clist_i); // min-heap (increasing values)
 

 int main() {

     int n = 10000000;

-    crand_uniform_i32_t dist;

-    crand_rng32_t rng = crand_rng32_init(1234);

-    dist = crand_uniform_i32_init(0, n);

+    cstc64_uniform_t dist;

+    cstc64_t rng = cstc64_init(1234);

+    dist = cstc64_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_i32(&rng, &dist));

+        cqueue_i_push(&queue, cstc64_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_i32(&rng, &dist);

+        int r = cstc64_uniform(&rng, &dist);

         if (r & 1)

             ++n, cqueue_i_push(&queue, r);

         else

diff --git a/examples/random.c b/examples/random.c
index 7f1be3f5..4834a974 100644
--- a/examples/random.c
+++ b/examples/random.c
@@ -6,62 +6,38 @@
 int main()

 {

     enum {R = 30};

-    const size_t N = 1000000000;

-    clock_t difference, before;    

+    const size_t N = 100000000;

+    clock_t diff, before;    

     uint64_t sum = 0;

-

     uint64_t seed = time(NULL);

-    crand_rng32_t pcg = crand_rng32_init(seed);

-    uint32_t range = crand_i32(&pcg) & ((1u << 28) - 1);

-    crand_uniform_i32_t dist0 = crand_uniform_i32_init(0, range);

-

-    printf("32 uniform: %u\n", dist0.range);

-    double fsum = 0;

-    before = clock();

-    c_forrange (N) {

-        fsum += (double) crand_uniform_i32(&pcg, &dist0) / dist0.range;

-    }

-    difference = clock() - before;

-    printf("%zu %f: %f secs\n", N, fsum / N, (float) difference / CLOCKS_PER_SEC);

-

-    pcg = crand_rng32_init(seed);

-    dist0 = crand_uniform_i32_init(0, range);

-    puts("32 unbiased");

-    fsum = 0;

-    before = clock();

-    c_forrange (N) {

-        fsum += (double) crand_unbiased_i32(&pcg, &dist0) / dist0.range;

-    }

-    difference = clock() - before;

-    printf("%zu %f: %f secs\n", N, fsum / N, (float) difference / CLOCKS_PER_SEC);

 

-    puts("64 uniform");

-    crand_rng64_t stc = crand_rng64_init(seed);

-    crand_uniform_i64_t dist1 = crand_uniform_i64_init(0, N);

+    cstc64_t stc = cstc64_init(seed);

+    cstc64_uniform_t dist1 = cstc64_uniform_init(0, N);

     sum = 0;

     before = clock();    

     c_forrange (N)  {

-        sum += crand_uniform_i64(&stc, &dist1);

+        sum += cstc64_uniform(&stc, &dist1);

     }

-    difference = clock() - before;

-    printf("%zu %f: %f secs\n", N, (double) sum / N, (float) difference / CLOCKS_PER_SEC);

-

+    diff = clock() - before;

+    printf("uniform: %zu %f: %f secs\n", N, (double) sum / N, (float) diff / CLOCKS_PER_SEC);

 

-    puts("normal distribution");

-    crand_normal_f64_t dist2 = crand_normal_f64_init(R / 2.0, R / 6.0);

+    cstc64_normalf_t dist2 = cstc64_normalf_init(R / 2.0, R / 6.0);

     size_t N2 = 10000000;

     int hist[R] = {0};

     sum = 0;

+    before = clock();

     c_forrange (N2)  {

-        int n = (int) (crand_normal_f64(&stc, &dist2) + 0.5);

+        int n = (int) (cstc64_normalf(&stc, &dist2) + 0.5);

         sum += n;

         if (n >= 0 && n < R) ++hist[n];

     }

-    

-    cstr_t bar = cstr_init();

+    diff = clock() - before;

+    printf("normal : %zu %f: %f secs\n", N, (double) sum / N2, (float) diff / CLOCKS_PER_SEC);

+

+    cstr_t bar = cstr_inits;

     c_forrange (i, int, R)  {

-        cstr_take(&bar, cstr_with_size(hist[i] * 25ull * R / N2, '*'));

-        printf("%2d %s\n", i, bar.str);

+        cstr_resize(&bar, hist[i] * 25ull * R / N2, '*');

+        printf("%3d %s\n", i, bar.str);

     }

     cstr_del(&bar);

 }
\ No newline at end of file
diff --git a/examples/rngtest.c b/examples/rngtest.c
index df3aa41c..76666484 100644
--- a/examples/rngtest.c
+++ b/examples/rngtest.c
@@ -6,49 +6,35 @@
 #endif

 

 

-#define NN 3000000000

+#define NN 1000000000

 

 int main(void)

 {

-    clock_t difference, before;

+    clock_t diff, before;

     uint64_t v;

     

-    crand_rng64_t stc = crand_rng64_init(time(NULL));

-    crand_uniform_i64_t idist = crand_uniform_i64_init(10, 20);

-    crand_uniform_f64_t fdist = crand_uniform_f64_init(10, 20);

+    cstc64_t stc = cstc64_init(time(NULL));

+    cstc64_uniform_t idist = cstc64_uniform_init(10, 20);

+    cstc64_uniformf_t fdist = cstc64_uniformf_init(10, 20);

 

-    c_forrange (30) printf("%02zd ", crand_uniform_i64(&stc, &idist)); 

-    puts("");

-

-    crand_rng32_t pcg = crand_rng32_init(time(NULL));

-    crand_uniform_i32_t i32dist = crand_uniform_i32_init(10, 20);

-    crand_uniform_f32_t f32dist = crand_uniform_f32_init(10, 20);

-

-    before = clock(); \

+    before = clock();

     v = 0;

     c_forrange (NN) {

-        //v += crand_i32(&pcg);

-        v += crand_uniform_i32(&pcg, &i32dist);

+        v += cstc64_rand(&stc);

     }

-    difference = clock() - before;

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

+    diff = clock() - before;

+    printf("stc64_rand: %.02f, %zu\n", (float) diff / CLOCKS_PER_SEC, v);

 

-    before = clock(); \

+    before = clock();

     v = 0;

     c_forrange (NN) {

-        //v += crand_i64(&stc) & 0xffffffff;

-        v += crand_uniform_i64(&stc, &idist);

+        v += cstc64_uniform(&stc, &idist);

     }

-    difference = clock() - before;

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

+    diff = clock() - before;

+    printf("stc64_uniform: %.02f, %zu\n\n", (float) diff / CLOCKS_PER_SEC, v);

 

-    c_forrange (8) printf("%d ", crand_uniform_i32(&pcg, &i32dist));

+    c_forrange (30) printf("%02zd ", cstc64_uniform(&stc, &idist)); 

     puts("");

-

-    

-    c_forrange (8) printf("%f ", crand_uniform_f32(&pcg, &f32dist));

-    puts("");

-    

-    c_forrange (8) printf("%f ", crand_uniform_f64(&stc, &fdist));

+    c_forrange (8) printf("%f ", cstc64_uniformf(&stc, &fdist));

     puts("");

 }
\ No newline at end of file
diff --git a/stc/cbitset.h b/stc/cbitset.h
index 715eddb4..2e1887a2 100644
--- a/stc/cbitset.h
+++ b/stc/cbitset.h
@@ -51,20 +51,21 @@ int main() {
 

 typedef struct cbitset { uint64_t* _arr; size_t size; } cbitset_t;

 

-STC_API cbitset_t cbitset_with_size(size_t size, bool value);

-STC_API cbitset_t cbitset_from_str(const char* str);

-STC_API char*     cbitset_to_str(cbitset_t set, char* str, size_t start, intptr_t stop);

-STC_API cbitset_t cbitset_clone(cbitset_t other);

-STC_API void      cbitset_resize(cbitset_t* self, size_t size, bool value);

-STC_API size_t    cbitset_count(cbitset_t set);

-STC_API bool      cbitset_is_disjoint(cbitset_t set, cbitset_t other);

-STC_API bool      cbitset_is_subset(cbitset_t set, cbitset_t other);

-STC_API bool      cbitset_is_superset(cbitset_t set, cbitset_t other);

-

-STC_INLINE size_t cbitset_size(cbitset_t set) {return set.size;}

+STC_API cbitset_t    cbitset_with_size(size_t size, bool value);

+STC_API cbitset_t    cbitset_from_str(const char* str);

+STC_API char*        cbitset_to_str(cbitset_t set, char* str, size_t start, intptr_t stop);

+STC_API cbitset_t    cbitset_clone(cbitset_t other);

+STC_API void         cbitset_resize(cbitset_t* self, size_t size, bool value);

+STC_API size_t       cbitset_count(cbitset_t set);

+STC_API bool         cbitset_is_disjoint(cbitset_t set, cbitset_t other);

+STC_API bool         cbitset_is_subset(cbitset_t set, cbitset_t other);

+STC_API bool         cbitset_is_superset(cbitset_t set, cbitset_t other);

+

+STC_INLINE size_t    cbitset_size(cbitset_t set) {return set.size;}

+#define              cbitset_inits {NULL, 0}

 

 STC_INLINE cbitset_t cbitset_init() {

-    cbitset_t bs = {NULL, 0}; return bs;

+    cbitset_t bs = cbitset_inits; return bs;

 }

 STC_INLINE void cbitset_del(cbitset_t* self) {

     c_free(self->_arr);

diff --git a/stc/clist.h b/stc/clist.h
index 768de474..f322fd00 100644
--- a/stc/clist.h
+++ b/stc/clist.h
@@ -36,10 +36,10 @@
 

     int main() {

         clist_ix list = clist_inits;

-        crand_rng32_t pcg = crand_rng32_init(12345);

+        cstc64_t rng = cstc64_init(12345);

         int n;

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

-            clist_ix_push_back(&list, crand_i32(&pcg));

+            clist_ix_push_back(&list, cstc64_rand(&rng) >> 32);

         n = 0;

         c_foreach (i, clist_ix, list)

             if (++n % 10000 == 0) printf("%8d: %10zd\n", n, i.val->value);

diff --git a/stc/cpqueue.h b/stc/cpqueue.h
index f7f8150c..2bdf2807 100644
--- a/stc/cpqueue.h
+++ b/stc/cpqueue.h
@@ -29,13 +29,13 @@
     using_cpqueue(f, cvec_f, >); // min-heap (increasing values)

 

     int main() {

-        crand_rng32_t gen = crand_rng32_init(1234);

-        crand_uniform_f32_t dist = crand_uniform_f32_init(10.0f, 100.0f);

+        cstc64_t gen = cstc64_init(1234);

+        cstc64_uniformf_t dist = cstc64_uniformf_init(10.0f, 100.0f);

 

         cpqueue_f queue = cpqueue_f_init();

         // Push ten million random numbers onto the queue.

         for (int i=0; i<10000000; ++i)

-            cpqueue_f_push(&queue, crand_uniform_f32(&gen, dist));

+            cpqueue_f_push(&queue, cstc64_uniformf(&gen, dist));

         // Extract the 100 smallest.

         for (int i=0; i<100; ++i) {

             printf("%f ", *cpqueue_f_top(queue));

diff --git a/stc/cqueue.h b/stc/cqueue.h
index 7ad87c62..6226b59b 100644
--- a/stc/cqueue.h
+++ b/stc/cqueue.h
@@ -32,19 +32,19 @@
 

     int main() {

         int n = 10000000;

-        crand_rng32_t gen = crand_rng32_init(1234);

-        crand_uniform_i32_t dist = crand_uniform_i32_init(gen, 0, n);

+        cstc64_t gen = cstc64_init(1234);

+        cstc64_uniform_t dist = cstc64_uniform_init(gen, 0, n);

 

         cqueue_i queue = cqueue_i_init();

 

         // Push ten million random numbers onto the queue.

         for (int i=0; i<n; ++i)

-            cqueue_i_push(&queue, crand_uniform_i32(&dist));

+            cqueue_i_push(&queue, cstc64_uniform(&dist));

 

         // Push or pop on the queue ten million times

         printf("%d\n", n);

         for (int i=n; i>0; --i) {

-            int r = crand_uniform_i32(&dist);

+            int r = cstc64_uniform(&dist);

             if (r & 1)

                 ++n, cqueue_i_push(&queue, r);

             else

diff --git a/stc/crandom.h b/stc/crandom.h
index edd31113..fa353345 100644
--- a/stc/crandom.h
+++ b/stc/crandom.h
@@ -24,176 +24,118 @@
 #define CRANDOM__H__

 

 /*

-    crand_rng32_t rng = crand_rng32_init(seed);

-    crand_uniform_f32_t fdist = crand_uniform_f32_init(rng, 1.0f, 6.0f);

-    crand_uniform_i32_t idist = crand_uniform_i32_init(rng, 1, 6);

+// cstc64: Pseudo-random number generator

+#include "stc/crandom.h"

+int main() {

+    uint64_t seed = 123456789;

+    cstc64_t rng = cstc64_init(seed);

+    cstc64_uniform_t dist1 = cstc64_uniform_init(1, 6);

+    cstc64_uniformf_t dist2 = cstc64_uniformf_init(1.0, 10.0);

+    cstc64_normalf_t dist3 = cstc64_normalf_init(1.0, 10.0);

     

-    uint32_t i = crand_i32(&rng);

-    int j = crand_uniform_i32(&idist);

-    float r = crand_uniform_f32(&fdist);

+    uint64_t i = cstc64_rand(&rng);

+    int64_t iu = cstc64_uniform(&rng, &dist1);

+    double xu = cstc64_uniformf(&rng, &dist2);

+    double xn = cstc64_normalf(&rng, &dist3);

+}

 */

 #include "ccommon.h"

 #include <string.h>

 #include <math.h>

 

-/* 32-BIT RANDOM NUMBER GENERATOR */

-

-typedef struct {uint64_t state[2];} crand_rng32_t;

-typedef struct {int32_t offset; uint32_t range;} crand_uniform_i32_t;

-typedef struct {float offset, range;} crand_uniform_f32_t;

-

-/* engine initializers */

-STC_API crand_rng32_t crand_rng32_with_seq(uint64_t seed, uint64_t seq);

-STC_INLINE crand_rng32_t crand_rng32_init(uint64_t seed) {

-    return crand_rng32_with_seq(seed, seed);

-}

-

-/* int random number generator, range [0, 2^32) */

-STC_API uint32_t crand_i32(crand_rng32_t* rng);

-

-STC_INLINE float crand_f32(crand_rng32_t* rng) {

-    union {uint32_t i; float f;} u = {0x3F800000u | (crand_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, (uint32_t) (high - low + 1)}; return dist;

-}

-STC_INLINE int32_t crand_uniform_i32(crand_rng32_t* rng, crand_uniform_i32_t* dist) {

-    return dist->offset + (int32_t) (((uint64_t) crand_i32(rng) * dist->range) >> 32);

-}

-/* https://github.com/lemire/fastrange */

-STC_API uint32_t crand_unbiased_i32(crand_rng32_t* rng, crand_uniform_i32_t* dist);

-

-/* 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_rng32_t* rng, crand_uniform_f32_t* dist) {

-    return dist->offset + crand_f32(rng) * dist->range;

-}

-

+typedef struct {uint64_t state[4];}                         cstc64_t;

+typedef struct {int64_t offset; uint64_t range, threshold;} cstc64_uniform_t;

+typedef struct {double offset, range;}                      cstc64_uniformf_t;

+typedef struct {double mean, stddev, next; bool has_next;}  cstc64_normalf_t;

 

-/* 64 BIT RANDOM NUMBER GENERATOR */

 

-typedef struct {uint64_t state[4];} crand_rng64_t;

-typedef struct {int64_t offset; uint64_t range;} crand_uniform_i64_t;

-typedef struct {double offset, range;} crand_uniform_f64_t;

-typedef struct {double mean, stddev, next; bool has_next;} crand_normal_f64_t;

-

-

-/* engine initializers */

-STC_API crand_rng64_t crand_rng64_with_seq(uint64_t seed, uint64_t seq);

-STC_INLINE crand_rng64_t crand_rng64_init(uint64_t seed) {

-    return crand_rng64_with_seq(seed, 1);

-}

 /* int random number generator, range [0, 2^64). PRNG copyright Tyge Løvset, NORCE Research, 2020 */

-STC_API uint64_t crand_i64(crand_rng64_t* rng);

+STC_API cstc64_t cstc64_init(uint64_t seed);

+STC_API cstc64_t cstc64_with_seq(uint64_t seed, uint64_t seq);

+STC_API uint64_t cstc64_rand(cstc64_t* rng);

 

-/* double random number in range [low, high). 52 bit resolution. */

-STC_INLINE double crand_f64(crand_rng64_t* rng) {

-    union {uint64_t i; double f;} u = {0x3FF0000000000000ull | (crand_i64(rng) >> 12)};

+/* double random number in range [low, high). */

+STC_INLINE double cstc64_randf(cstc64_t* rng) {

+    union {uint64_t i; double f;} u = {0x3FF0000000000000ull | (cstc64_rand(rng) >> 12)};

     return u.f - 1.0;

 }

 

-/* int random number generator in range [low, high] */

-STC_INLINE crand_uniform_i64_t crand_uniform_i64_init(int64_t low, int64_t high) {

-    crand_uniform_i64_t dist = {low, (uint64_t) (high - low + 1)}; return dist;

-}

-

-#if defined(_MSC_VER) && defined(_WIN64)

-#include <intrin.h>

-#endif

-

-STC_INLINE int64_t crand_uniform_i64(crand_rng64_t* rng, crand_uniform_i64_t* dist) {

-    #if defined(__SIZEOF_INT128__)

-        return dist->offset + (int64_t) (((__uint128_t) crand_i64(rng) * dist->range) >> 64);

-    #elif defined(_MSC_VER) && defined(_WIN64)

-        uint64_t hi; _umul128(crand_i64(rng), dist->range, &hi); return dist->offset + hi;

-    #else

-        return dist->offset + crand_i64(rng) % dist->range; // slower

-    #endif

-}

+/* integer uniform distributed RNG, range [low, high]. */

+STC_API cstc64_uniform_t cstc64_uniform_init(int64_t low, int64_t high);

+STC_API int64_t          cstc64_uniform(cstc64_t* rng, cstc64_uniform_t* dist);

 

-STC_INLINE crand_uniform_f64_t crand_uniform_f64_init(double low, double high) {

-    crand_uniform_f64_t dist = {low, high - low}; return dist;

+/* double uniform distributed RNG, range [low, high). */

+STC_INLINE cstc64_uniformf_t cstc64_uniformf_init(double low, double high) {

+    cstc64_uniformf_t dist = {low, high - low}; return dist;

 }

-STC_INLINE double crand_uniform_f64(crand_rng64_t* rng, crand_uniform_f64_t* dist) {

-    return dist->offset + crand_f64(rng) * dist->range;

+STC_INLINE double cstc64_uniformf(cstc64_t* rng, cstc64_uniformf_t* dist) {

+    return dist->offset + cstc64_randf(rng) * dist->range;

 }

 

-STC_INLINE crand_normal_f64_t crand_normal_f64_init(double mean, double stddev) {

-    crand_normal_f64_t dist = {mean, stddev, 0.0, false}; return dist;

+/* double normal distributed RNG. */

+STC_INLINE cstc64_normalf_t cstc64_normalf_init(double mean, double stddev) {

+    cstc64_normalf_t dist = {mean, stddev, 0.0, false}; return dist;

 }

-STC_API double crand_normal_f64(crand_rng64_t* rng, crand_normal_f64_t* dist);

+STC_API double cstc64_normalf(cstc64_t* rng, cstc64_normalf_t* dist);

 

 

 #if !defined(STC_HEADER) || defined(STC_IMPLEMENTATION)

 

 #if defined(_MSC_VER)

   #pragma warning(push)

-  #pragma warning(disable: 4146)

+  //#pragma warning(disable: 4146)

 #endif

 

-/* PRNG PCG32 https://www.pcg-random.org/download.html */

-STC_DEF crand_rng32_t crand_rng32_with_seq(uint64_t seed, uint64_t seq) {

-    crand_rng32_t rng = {{0u, (seq << 1u) | 1u}}; /* inc must be odd */

-    crand_i32(&rng);

-    rng.state[0] += seed;

-    crand_i32(&rng);

-    return rng;

-}

-STC_DEF uint32_t crand_i32(crand_rng32_t* rng) {

-    uint64_t old = rng->state[0];

-    rng->state[0] = old * 6364136223846793005ull + rng->state[1];

-    uint32_t xors = (uint32_t) (((old >> 18u) ^ old) >> 27u);

-    uint32_t rot = old >> 59u;

-    return (xors >> rot) | (xors << ((-rot) & 31));

-}

+/* PRNG stc64: 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.

+ */

 

-/* PRNG STC64: copyright Tyge Løvset, NORCE Research, 2020 */

-/* Extremely fast PRNG suited for parallel usage with Weyl-sequence parameter. */

-/* Even faster than sfc64: updates only 192bit state. Better for parallel processing: */

-/* Guarantees 2^63 unique threads with minimum 2^64 period length ~ 2^160 average period. */

-/* Tested with PractRand to 8 TB output: no issues */

-STC_DEF crand_rng64_t crand_rng64_with_seq(uint64_t seed, uint64_t seq) {

-    crand_rng64_t rng = {{seed, seed, seed, (seq << 1u) | 1u}}; /* increment must be odd */

-    for (int i = 0; i < 12; ++i) crand_i64(&rng);

+STC_DEF cstc64_t cstc64_init(uint64_t seed) {

+    return cstc64_with_seq(seed, 0x3504f333d3aa0b34);

+}

+STC_DEF cstc64_t cstc64_with_seq(uint64_t seed, uint64_t seq) {

+    cstc64_t rng = {{seed, seed, seed, (seq << 1u) | 1u}};

+    for (int i = 0; i < 8; ++i) cstc64_rand(&rng);

     return rng;

 }

-STC_DEF uint64_t crand_i64(crand_rng64_t* rng) {

+STC_DEF uint64_t cstc64_rand(cstc64_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;

+    const uint64_t result = s[0] ^ (s[2] += s[3]|1);

+    s[0] = (s[1] + (s[1] << LSHIFT)) ^ (s[1] >> RSHIFT);

+    s[1] = ((s[1] << LROT) | (s[1] >> (64 - LROT))) + result;

     return result;

 }

 

-/* Unbiased uniform https://github.com/lemire/fastrange */

-STC_DEF uint32_t crand_unbiased_i32(crand_rng32_t* rng, crand_uniform_i32_t* dist) {

-    uint32_t r = dist->range;

-    uint64_t m = (uint64_t) crand_i32(rng) * r;

-    uint32_t l = (uint32_t) m;

-    if (l < r) {

-        uint32_t t = -r;

-        if (t >= r) if ((t -= r) >= r) t %= r;

-        while (l < t) l = (uint32_t) (m = (uint64_t) crand_i32(rng) * r);

+/* unbiased integer random number generator in range [low, high] */

+STC_DEF cstc64_uniform_t cstc64_uniform_init(int64_t low, int64_t high) {

+    cstc64_uniform_t dist = {low, (uint64_t) (high - low + 1)};

+    dist.threshold = (uint64_t)(-dist.range) % dist.range;

+    return dist;

+}

+STC_DEF int64_t cstc64_uniform(cstc64_t* rng, cstc64_uniform_t* dist) {

+    for (;;) {

+        uint64_t r = cstc64_rand(rng);

+        if (r >= dist->threshold) return dist->offset + (r % dist->range);

     }

-    return dist->offset + (m >> 32);

 }

 

 /* Marsaglia polar method for gaussian distribution. */

-STC_DEF double crand_normal_f64(crand_rng64_t* rng, crand_normal_f64_t* dist) {

+STC_DEF double cstc64_normalf(cstc64_t* rng, cstc64_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_f64(rng) - 1.0;

-        u2 = 2.0 * crand_f64(rng) - 1.0;

+        u1 = 2.0 * cstc64_randf(rng) - 1.0;

+        u2 = 2.0 * cstc64_randf(rng) - 1.0;

         s = u1*u1 + u2*u2;

     } while (s >= 1.0 || s == 0.0);

     m = sqrt(-2.0 * log(s) / s);