Some API changes in crandom. Added crandom_normal_f64() - normal distributed distribution.

5 files changed, 51 insertions, 41 deletions

diff --git a/examples/list.c b/examples/list.c
index b981b1ec..3018a1f1 100644
--- a/examples/list.c
+++ b/examples/list.c
@@ -8,7 +8,7 @@ int main() {
     int k, n = 100000;

     clist_fx list = clist_init;

     crandom_eng64_t eng = crandom_eng64_init(time(NULL));

-    crandom_uniform_f64_t dist = crandom_uniform_f64_init(0.0f, n);

+    crandom_distrib_f64_t dist = crandom_uniform_f64_init(0.0f, n);

     

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

         clist_fx_push_back(&list, crandom_uniform_f64(&eng, dist));

diff --git a/examples/priority.c b/examples/priority.c
index 70d89f2a..780e5809 100644
--- a/examples/priority.c
+++ b/examples/priority.c
@@ -10,7 +10,7 @@ declare_cvec_pqueue(i, >); // min-heap (increasing values)
 

 int main() {

     crandom_eng32_t pcg = crandom_eng32_init(time(NULL));

-    crandom_uniform_i32_t dist = crandom_uniform_i32_init(0, 100000000);

+    crandom_distrib_i32_t dist = crandom_uniform_i32_init(0, 100000000);

     cvec_i heap = cvec_init;

 

     // Push ten million random numbers to priority queue

diff --git a/examples/rngbirthday.c b/examples/rngbirthday.c
index 595bacdb..22ecdfae 100644
--- a/examples/rngbirthday.c
+++ b/examples/rngbirthday.c
@@ -38,7 +38,7 @@ void distribution(void)
     const size_t N = 1ull << 28, M = 1ull << 9; // 1ull << 10;

     cmap_x map = cmap_x_with_capacity(M);

     clock_t now = clock();

-    crandom_uniform_i32_t dist = crandom_uniform_i32_init(0, M);

+    crandom_distrib_i32_t dist = crandom_uniform_i32_init(0, M);

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

         ++cmap_x_insert(&map, crandom_uniform_i32(&rng, dist), 0)->value;

     }

diff --git a/examples/rngtest.c b/examples/rngtest.c
index e19147c7..ca037292 100644
--- a/examples/rngtest.c
+++ b/examples/rngtest.c
@@ -32,15 +32,15 @@ int main(void)
     difference = clock() - before;

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

 

-    crandom_uniform_i32_t i32dist = crandom_uniform_i32_init(10, 20);

+    crandom_distrib_i32_t i32dist = crandom_uniform_i32_init(10, 20);

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

     puts("");

 

-    crandom_uniform_f32_t f32dist = crandom_uniform_f32_init(10, 20);

+    crandom_distrib_f32_t f32dist = crandom_uniform_f32_init(10, 20);

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

     puts("");

 

-    crandom_uniform_f64_t fdist = crandom_uniform_f64_init(10, 20);

+    crandom_distrib_f64_t fdist = crandom_uniform_f64_init(10, 20);

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

     puts("");

 }
\ No newline at end of file
diff --git a/stc/crandom.h b/stc/crandom.h
index 5d547bed..b73dd11b 100644
--- a/stc/crandom.h
+++ b/stc/crandom.h
@@ -26,10 +26,11 @@
 

 #include "cdefs.h"

 #include <string.h>

+#include <math.h>

 /*

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

+    crandom_distrib_f32_t fdist = crandom_uniform_f32_init(1.0f, 6.0f);

+    crandom_distrib_i32_t idist = crandom_uniform_i32_init(1, 6);

     

     uint32_t i = crandom_i32(&eng);

     int j = crandom_uniform_i32(&eng, idist);

@@ -37,8 +38,8 @@
 */

 

 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;

+typedef struct {int32_t offset, range;} crandom_distrib_i32_t;

+typedef struct {float offset, range;} crandom_distrib_f32_t;

 

 /* 32 bit random number generator engine */

 STC_API crandom_eng32_t crandom_eng32_with_seq(uint64_t seed, uint64_t seq);

@@ -55,24 +56,24 @@ STC_INLINE float crandom_f32(crandom_eng32_t* rng) {
 }

 

 /* 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 crandom_distrib_i32_t crandom_uniform_i32_init(int32_t low, int32_t high) {

+    crandom_distrib_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_i32(rng) * dist.range) >> 32);

+STC_INLINE int32_t crandom_uniform_i32(crandom_eng32_t* rng, crandom_distrib_i32_t dist) {

+    return dist.offset + (int32_t) (((uint64_t) crandom_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 crandom_distrib_f32_t crandom_uniform_f32_init(float low, float high) {

+    crandom_distrib_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_f32(rng) * dist.range;

+STC_INLINE float crandom_uniform_f32(crandom_eng32_t* rng, crandom_distrib_f32_t dist) {

+    return dist.offset + crandom_f32(rng) * dist.range;

 }

 

 

 typedef struct {uint64_t state[4];} crandom_eng64_t;

-typedef struct {double min, range;} crandom_uniform_f64_t;

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

 

 /* 64 bit random number generator engine */

 STC_API crandom_eng64_t crandom_eng64_with_seq(uint64_t seed, uint64_t seq);

@@ -82,18 +83,23 @@ STC_INLINE crandom_eng64_t crandom_eng64_init(uint64_t seed) {
 /* int random number generator, range [0, 2^64) */

 STC_API uint64_t crandom_i64(crandom_eng64_t* rng);

 

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

 STC_INLINE double crandom_f64(crandom_eng64_t* rng) {

     union {uint64_t i; double f;} u = {0x3FF0000000000000ull | (crandom_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 crandom_distrib_f64_t crandom_uniform_f64_init(double low, double high) {

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

+}

+STC_INLINE double crandom_uniform_f64(crandom_eng64_t* rng, crandom_distrib_f64_t dist) {

+    return dist.offset + crandom_f64(rng) * dist.range;

 }

-STC_INLINE double crandom_uniform_f64(crandom_eng64_t* rng, crandom_uniform_f64_t dist) {

-    return dist.min + crandom_f64(rng) * dist.range;

+

+STC_INLINE crandom_distrib_f64_t crandom_normal_f64_init(double mean, double std_dev) {

+    crandom_distrib_f64_t dist = {mean, std_dev}; return dist;

 }

+STC_API double crandom_normal_f64(crandom_eng64_t* rng, crandom_distrib_f64_t dist);

 

 

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

@@ -122,22 +128,10 @@ STC_API crandom_eng64_t crandom_eng64_with_seq(uint64_t seed, uint64_t seq) {
     return rng;

 }

 

-#ifdef STC_USE_SFC64

-/* SFC64 random number generator: http://pracrand.sourceforge.net */

-STC_API uint64_t crandom_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;

-}

-#else

-/* My own PRNG inspired by SFC64: Faster and has Weyl-sequence parameter. */

-/* Copyright Tyge Løvset, 2020 */

-/* Faster: updates only 192bit state. Parallel: Ensures unique sequence per seq (2^63) */

-/* Minimum period is 2^64 per seq, average ~ 2^127 per seq */

+/* PRNG copyright Tyge Løvset, NORCE Research, 2020 */

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

+/* Updates only 192bit state. Parallel: Ensures unique sequence per seq (2^63) */

+/* Minimum period is 2^64 per seq, but high average per Weyl seq. */

 STC_API uint64_t crandom_i64(crandom_eng64_t* rng) {

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

     uint64_t *s = rng->state;

@@ -146,7 +140,23 @@ STC_API uint64_t crandom_i64(crandom_eng64_t* rng) {
     s[1] = ((b << LROT) | (b >> (64 - LROT))) + result;

     return result;

 }

-#endif

+

+STC_API double crandom_normal_f64(crandom_eng64_t* rng, crandom_distrib_f64_t dist) {

+    static bool spare = false; /* Marsaglia polar method: */

+    static double u2; double u1, s, m;

+    if (spare) { 

+        spare = false;

+        return u2 * dist.range + dist.offset;

+    }

+    do {

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

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

+        s = u1*u1 + u2*u2;

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

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

+    u2 *= m; spare = true;

+    return (u1 * m) * dist.range + dist.offset;

+}

 

 #endif

 #endif