Reverted back API to easier naming scheme of crandom.h types and methods.

Highly optimized generation of unbiased uniform bounded random numbers. Only 33% overhead.

15 files changed, 144 insertions, 128 deletions

diff --git a/docs/cpqueue_api.md b/docs/cpqueue_api.md
index 1e92277a..ac33e68f 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;
-    cstc64_t rng = cstc64_init(1234);
-    cstc64_uniform_t dist = cstc64_uniform_init(0, N * 10);
+    crand_t rng = crand_init(1234);
+    crand_uniform_t dist = crand_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, cstc64_uniform(&rng, &dist));
+        cpqueue_i_push(&heap, crand_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 57697067..53e8a0d1 100644
--- a/docs/crandom_api.md
+++ b/docs/crandom_api.md
@@ -5,12 +5,12 @@ This describes the API of module **crand**. It contains *pcg32* created by Melis
 
 ## Types
 
-| 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     |
+| Name               | Type definition                            | Used to represent...         |
+|:-------------------|:-------------------------------------------|:-----------------------------|
+| `crand_t`          | `struct {uint64_t state[4];}`              | The PRNG engine type         |
+| `crand_uniform_t`  | `struct {int64_t offset; uint64_t range;}` | Integer uniform distribution |
+| `crand_uniformf_t` | `struct {double offset, range;}`           | Real number uniform distr.   |
+| `crand_normalf_t`  | `struct {double mean, stddev;}`            | Normal distribution type     |
 
 ## Header file
 
@@ -22,28 +22,33 @@ All cstr definitions and prototypes may be included in your C source file by inc
 ## Methods
 
 ```c
- 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)     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-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*].
+RNG, around 68% of the values fall 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.
+**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\*\** compiled with 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 by incrementing the Weyl-increment by 2, it initializes
+a new unique 2^64 *minimum* length period. Note that for each Weyl-increment, the expected period
+length is about 2^126. For a single thread, a minimum period of 2^127 is generated if the
+Weyl-increment is incremented by 2 every 2^64 number output.
+
+*crand* 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
@@ -75,13 +80,13 @@ int main()
 
     // Setup random engine with normal distribution.
     uint64_t seed = time(NULL);
-    cstc64_t rng = cstc64_init(seed);
-    cstc64_normalf_t dist = cstc64_normalf_init(Mean, StdDev);
+    crand_t rng = crand_init(seed);
+    crand_normalf_t dist = crand_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( cstc64_normalf(&rng, &dist) );
+        int index = (int) round( crand_normalf(&rng, &dist) );
         cmap_i_emplace(&mhist, index, 0).first->second += 1;
     }
 
diff --git a/examples/benchmark.cpp b/examples/benchmark.cpp
index 4830c9ba..206da5ce 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;

 

-cstc64_t rng;

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

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

+crand_t rng;

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

+#define RAND(N) (crand_next(&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 fedb2c5c..ac05a5f5 100644
--- a/examples/birthday.c
+++ b/examples/birthday.c
@@ -14,12 +14,12 @@ const static uint64_t mask = (1ull << 52) - 1;
 

 void repeats(void)

 {

-    cstc64_t rng = cstc64_init(seed);

+    crand_t rng = crand_init(seed);

     cmap_ic m = cmap_ic_init();

     cmap_ic_reserve(&m, N);

     clock_t now = clock();

     c_forrange (i, N) {

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

+        uint64_t k = crand_next(&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)

 {

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

+    crand_t rng = crand_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();

-    cstc64_uniform_t dist = cstc64_uniform_init(0, M);

+    crand_uniform_t dist = crand_uniform_init(0, M);

 

     c_forrange (N) {

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

+        ++cmap_x_emplace(&map, crand_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 5303dd0c..a3842aae 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);

-    cstc64_t rng = cstc64_init(seed);

-    cstc64_normalf_t dist = cstc64_normalf_init(Mean, StdDev);

+    crand_t rng = crand_init(seed);

+    crand_normalf_t dist = crand_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( cstc64_normalf(&rng, &dist) );

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

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

     }

 

diff --git a/examples/heap.c b/examples/heap.c
index 166dd6ba..83f9d97f 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);

-    cstc64_t rng;

+    crand_t rng;

     int N = 3000000, M = 100;

 

     cpqueue_f pq = cpqueue_f_init();

 

-    rng = cstc64_init(seed);

+    rng = crand_init(seed);

     clock_t start = clock();

     c_forrange (i, int, N)

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

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

 

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

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

+        cpqueue_f_push(&pq, (float) crand_nextf(&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 fbfa8830..4f6d16e1 100644
--- a/examples/list.c
+++ b/examples/list.c
@@ -9,11 +9,11 @@ int main() {
     const int n = 2000000;

 

     clist_fx list = clist_inits;

-    cstc64_t eng = cstc64_init(1234);

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

+    crand_t rng = crand_init(1234);

+    crand_uniformf_t dist = crand_uniformf_init(100.0f, n);

     int m = 0;

     c_forrange (i, int, n)

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

+        clist_fx_push_back(&list, crand_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 e81cd8a3..f6b4ecf0 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;

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

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

+    crand_t rng = crand_init(time(NULL));

+    crand_uniform_t dist = crand_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, cstc64_uniform(&pcg, &dist));

+        cpqueue_i_push(&heap, crand_uniform(&rng, &dist));

 

     // push some negative numbers too.

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

 

     c_forrange (N)

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

+        cpqueue_i_push(&heap, crand_uniform(&rng, &dist));

 

 

     // Extract the hundred smallest.

diff --git a/examples/queue.c b/examples/queue.c
index 0ff29736..57362493 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;

-    cstc64_uniform_t dist;

-    cstc64_t rng = cstc64_init(1234);

-    dist = cstc64_uniform_init(0, n);

+    crand_uniform_t dist;

+    crand_t rng = crand_init(1234);

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

+        cqueue_i_push(&queue, crand_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 = cstc64_uniform(&rng, &dist);

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

         if (r & 1)

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

         else

diff --git a/examples/random.c b/examples/random.c
index 0cf9da6f..de1fe778 100644
--- a/examples/random.c
+++ b/examples/random.c
@@ -7,18 +7,18 @@
 int main()

 {

     enum {R = 30};

-    const size_t N = 500000000;

-    uint64_t seed = time(NULL);

-    cstc64_t stc = cstc64_init(seed);

+    const size_t N = 1000000000;

+    uint64_t seed = 1234; // time(NULL);

+    crand_t rng = crand_init(seed);

 

     uint64_t sum = 0;

 

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

+    crand_normalf_t dist2 = crand_normalf_init(R / 2.0, R / 6.0);

     size_t N2 = 10000000;

     int hist[R] = {0};

     sum = 0;

     c_forrange (N2)  {

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

+        int n = round((crand_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 += cstc64_rand(&stc);

+        sum += crand_next(&rng);

     }

     diff = clock() - before;

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

 

-    cstc64_uniform_t dist1 = cstc64_uniform_init(0, 1000);

+    crand_uniform_t dist1 = crand_uniform_init(0, 1000);

     sum = 0;

     before = clock();    

     c_forrange (N)  {

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

+        sum += crand_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/examples/rngtest.c b/examples/rngtest.c
index 76666484..db7a8093 100644
--- a/examples/rngtest.c
+++ b/examples/rngtest.c
@@ -13,14 +13,14 @@ int main(void)
     clock_t diff, before;

     uint64_t v;

     

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

+    crand_t rng = crand_init(time(NULL));

+    crand_uniform_t idist = crand_uniform_init(10, 20);

+    crand_uniformf_t fdist = crand_uniformf_init(10, 20);

 

     before = clock();

     v = 0;

     c_forrange (NN) {

-        v += cstc64_rand(&stc);

+        v += crand_next(&rng);

     }

     diff = clock() - before;

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

@@ -28,13 +28,13 @@ int main(void)
     before = clock();

     v = 0;

     c_forrange (NN) {

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

+        v += crand_uniform(&rng, &idist);

     }

     diff = clock() - before;

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

 

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

+    c_forrange (30) printf("%02zd ", crand_uniform(&rng, &idist)); 

     puts("");

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

+    c_forrange (8) printf("%f ", crand_uniformf(&rng, &fdist));

     puts("");

 }
\ No newline at end of file
diff --git a/stc/clist.h b/stc/clist.h
index f322fd00..6594aa8c 100644
--- a/stc/clist.h
+++ b/stc/clist.h
@@ -36,10 +36,10 @@
 

     int main() {

         clist_ix list = clist_inits;

-        cstc64_t rng = cstc64_init(12345);

+        crand_t rng = crand_init(12345);

         int n;

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

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

+            clist_ix_push_back(&list, crand_next(&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 2bdf2807..c1dfde17 100644
--- a/stc/cpqueue.h
+++ b/stc/cpqueue.h
@@ -29,13 +29,13 @@
     using_cpqueue(f, cvec_f, >); // min-heap (increasing values)

 

     int main() {

-        cstc64_t gen = cstc64_init(1234);

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

+        crand_t rng = crand_init(1234);

+        crand_uniformf_t dist = crand_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, cstc64_uniformf(&gen, dist));

+            cpqueue_f_push(&queue, crand_uniformf(&rng, 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 6226b59b..f385aa4a 100644
--- a/stc/cqueue.h
+++ b/stc/cqueue.h
@@ -32,19 +32,19 @@
 

     int main() {

         int n = 10000000;

-        cstc64_t gen = cstc64_init(1234);

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

+        crand_t rng = crand_init(1234);

+        crand_uniform_t dist = crand_uniform_init(rng, 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, cstc64_uniform(&dist));

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

 

         // Push or pop on the queue ten million times

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

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

-            int r = cstc64_uniform(&dist);

+            int r = crand_uniform(&dist);

             if (r & 1)

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

             else

diff --git a/stc/crandom.h b/stc/crandom.h
index f0413abd..5db3c160 100644
--- a/stc/crandom.h
+++ b/stc/crandom.h
@@ -24,64 +24,64 @@
 #define CRANDOM__H__

 

 /*

-// cstc64: Pseudo-random number generator

+// crand: 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);

+    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 = cstc64_rand(&rng);

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

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

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

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

 }

 */

 #include "ccommon.h"

 #include <string.h>

 #include <math.h>

 

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

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

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

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

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

 

 

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

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

+STC_API crand_t crand_init(uint64_t seed);

+STC_API crand_t crand_with_seq(uint64_t seed, uint64_t seq);

+STC_API uint64_t crand_next(crand_t* rng);

 

 /* 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)};

+STC_INLINE double crand_nextf(crand_t* rng) {

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

     return u.f - 1.0;

 }

 

 /* 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_API crand_uniform_t crand_uniform_init(int64_t low, int64_t high);

+STC_API int64_t         crand_uniform(crand_t* rng, crand_uniform_t* 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 crand_uniformf_t crand_uniformf_init(double low, double high) {

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

 }

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

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

+STC_INLINE double crand_uniformf(crand_t* rng, crand_uniformf_t* dist) {

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

 }

 

 /* 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_INLINE crand_normalf_t crand_normalf_init(double mean, double stddev) {

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

 }

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

+STC_API double crand_normalf(crand_t* rng, crand_normalf_t* dist);

 

 

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

 

-/* PRNG stc64: Tyge Løvset, NORCE Research, 2020.

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

@@ -91,46 +91,58 @@ STC_API double cstc64_normalf(cstc64_t* rng, cstc64_normalf_t* dist);
  * and simple correlation tests, i.e. interleaved streams with one-bit diff state.

  */

 

-STC_DEF cstc64_t cstc64_init(uint64_t seed) {

-    return cstc64_with_seq(seed, 0x3504f333d3aa0b34);

+STC_DEF crand_t crand_init(uint64_t seed) {

+    return crand_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);

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

     return rng;

 }

-STC_DEF uint64_t cstc64_rand(cstc64_t* rng) {

+STC_DEF uint64_t crand_next(crand_t* rng) {

     enum {LROT = 24, RSHIFT = 11, LSHIFT = 3};

     uint64_t *s = rng->state;

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

+    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;

 }

 

 /* 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)};

+

+STC_DEF crand_uniform_t crand_uniform_init(int64_t low, int64_t high) {

+    crand_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);

-    }

+

+#if defined(__SIZEOF_INT128__)

+    #define cmul128(a, b, m) *(__uint128_t *)m = (__uint128_t)a * b

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

+    #include <intrin.h>

+    #define cmul128(a, b, m) m[0] = _umul128(a, b, &m[1])

+#elif defined(__x86_64__)

+    #define cmul128(a, b, m) \

+        asm("mulq %[rhs]" : "=a" (m[0]), "=d" (m[1]) \

+                          : [lhs] "0" (a), [rhs] "rm" (b))

+#endif

+

+STC_DEF int64_t crand_uniform(crand_t* rng, crand_uniform_t* dist) {

+    uint64_t m[2]; 

+    do { cmul128(crand_next(rng), dist->range, m); } while (m[0] < dist->threshold);

+    return dist->offset + m[1];

 }

 

 /* Marsaglia polar method for gaussian distribution. */

-STC_DEF double cstc64_normalf(cstc64_t* rng, cstc64_normalf_t* dist) {

+STC_DEF double crand_normalf(crand_t* rng, crand_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 * cstc64_randf(rng) - 1.0;

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

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

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

         s = u1*u1 + u2*u2;

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

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

@@ -138,6 +150,5 @@ STC_DEF double cstc64_normalf(cstc64_t* rng, cstc64_normalf_t* dist) {
     return (u1 * m) * dist->stddev + dist->mean;

 }

 

-

 #endif

 #endif