Refactored crandom. Removed 32-bit version. Made normal dist inline - does not require math lib if not used. Added stc64_randomf()

2 files changed, 43 insertions, 82 deletions

diff --git a/docs/crandom_api.md b/docs/crandom_api.md
index b099eeaf..896f1c83 100644
--- a/docs/crandom_api.md
+++ b/docs/crandom_api.md
@@ -43,6 +43,7 @@ All crandom definitions and prototypes are available by including a single heade
 ```c
 void                stc64_srandom(uint64_t seed);                            // seed global rng
 uint64_t            stc64_random(void);                                      // range [0, 2^64 - 1]
+double              stc64_randomf(void);                                     // range [0.0, 1.0)
 
 stc64_t             stc64_init(uint64_t seed);
 stc64_t             stc64_with_seq(uint64_t seed, uint64_t seq);             // init with stream
diff --git a/include/stc/crandom.h b/include/stc/crandom.h
index 7543f7b9..97e0275d 100644
--- a/include/stc/crandom.h
+++ b/include/stc/crandom.h
@@ -44,9 +44,7 @@ int main() {
 #include <math.h>

 

 typedef struct stc64 { uint64_t state[5]; } stc64_t;

-typedef struct stc32 { uint32_t state[5]; } stc32_t;

 typedef struct stc64_uniform { int64_t lower; uint64_t range, threshold; } stc64_uniform_t;

-typedef struct stc32_uniform { int32_t lower; uint32_t range, threshold; } stc32_uniform_t;

 typedef struct stc64_uniformf { double lower, range; } stc64_uniformf_t;

 typedef struct stc64_normalf { double mean, stddev, next; unsigned has_next; } stc64_normalf_t;

 

@@ -61,37 +59,28 @@ typedef struct stc64_normalf { double mean, stddev, next; unsigned has_next; } s
  * PractRand to multiple TB input.

  */

 

-/* Global STC64 PRNG */

+/* Global STC64 PRNGs */

 STC_API void     stc64_srandom(uint64_t seed);

 STC_API uint64_t stc64_random(void);

-

-STC_API uint64_t stc64_rand(stc64_t* rng);

-STC_API uint32_t stc32_rand(stc32_t* rng);

+STC_API double   stc64_randomf(void);

 

 /* Init with sequence number */

-STC_API stc64_t  stc64_with_seq(uint64_t seed, uint64_t seq);

-STC_API stc32_t  stc32_with_seq(uint32_t seed, uint32_t seq);

-

-/* Int uniform distributed RNG, range [low, high]. */

-STC_API stc64_uniform_t stc64_uniform_init(int64_t low, int64_t high);

-STC_API stc32_uniform_t stc32_uniform_init(int32_t low, int32_t high);

-

-/* Float64 uniform distributed RNG, range [low, high). */

-STC_API stc64_uniformf_t stc64_uniformf_init(double low, double high);

+STC_API stc64_t stc64_with_seq(uint64_t seed, uint64_t seq);

 

-/* Normal distribution (double) */

-STC_API stc64_normalf_t stc64_normalf_init(double mean, double stddev);

-STC_API double          stc64_normalf(stc64_t* rng, stc64_normalf_t* dist);

-

-/* Unbiased bounded uniform distribution. */

-STC_API int64_t stc64_uniform(stc64_t* rng, stc64_uniform_t* d);

+/* Unbiased bounded uniform distribution. range [low, high] */

+STC_API int64_t stc64_uniform(stc64_t* rng, stc64_uniform_t* dist);

 

 STC_INLINE stc64_t stc64_init(uint64_t seed) {

     return stc64_with_seq(seed, seed + 0x3504f333d3aa0b37);

 }

 

-STC_INLINE stc32_t stc32_init(uint32_t seed) {

-    return stc32_with_seq(seed, seed + 0xd3aa0b37);

+STC_INLINE uint64_t stc64_rand(stc64_t* rng) {

+    uint64_t *s = rng->state; enum {LR=24, RS=11, LS=3};

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

+    s[0] = s[1] ^ (s[1] >> RS);

+    s[1] = s[2] + (s[2] << LS);

+    s[2] = ((s[2] << LR) | (s[2] >> (64 - LR))) + result;

+    return result;

 }

 

 /* Float64 random number in range [0.0, 1.0). */

@@ -105,10 +94,29 @@ STC_INLINE double stc64_uniformf(stc64_t* rng, stc64_uniformf_t* dist) {
     return stc64_randf(rng)*dist->range + dist->lower;

 }

 

-STC_INLINE int32_t stc32_uniform(stc32_t* rng, stc32_uniform_t* d) {

-    uint64_t val;

-    do { val = stc32_rand(rng) * (uint64_t)d->range; } while ((uint32_t)val < d->threshold);

-    return d->lower + (val >> 32);

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

+STC_INLINE stc64_uniformf_t stc64_uniformf_init(double low, double high) {

+    return c_make(stc64_uniformf_t){low, high - low};

+}

+

+/* Marsaglia polar method for gaussian/normal distribution, float64. */

+STC_INLINE stc64_normalf_t stc64_normalf_init(double mean, double stddev) {

+    return c_make(stc64_normalf_t){mean, stddev, 0.0, 0};

+}

+

+/* Normal distribution PRNG */

+STC_INLINE double stc64_normalf(stc64_t* rng, stc64_normalf_t* dist) {

+    double u1, u2, s, m;

+    if (dist->has_next++ & 1)

+        return dist->next * dist->stddev + dist->mean;

+    do {

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

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

+        s = u1*u1 + u2*u2;

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

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

+    dist->next = u2 * m;

+    return (u1 * m) * dist->stddev + dist->mean;

 }

 

 /* -------------------------- IMPLEMENTATION ------------------------- */

@@ -129,20 +137,18 @@ STC_DEF uint64_t stc64_random(void) {
     return stc64_rand(&stc64_global);

 }

 

+STC_DEF double stc64_randomf(void) {

+    return stc64_randf(&stc64_global);

+}

+

 /* rng.state[4] must be odd */

 STC_DEF stc64_t stc64_with_seq(uint64_t seed, uint64_t seq) {

-    stc64_t rng = {{seed, seed+0x26aa069ea2fb1a4d, seed+0x70c72c95cd592d04,

-                          seed+0x504f333d3aa0b359, (seq << 1) | 1}};

+    stc64_t rng = {{seed+0x26aa069ea2fb1a4d, seed+0x70c72c95cd592d04,

+                    seed+0x504f333d3aa0b359, seed, seed<<1 | 1}};

     for (int i = 0; i < 6; ++i) stc64_rand(&rng);

     return rng;

 }

 

-STC_DEF stc32_t stc32_with_seq(uint32_t seed, uint32_t seq) {

-    stc32_t rng = {{seed, seed+0x26aa069e, seed+0xa2fb1a4d,

-                          seed+0x70c72c95, (seq << 1) | 1}};

-    for (int i = 0; i < 6; ++i) stc32_rand(&rng);

-    return rng;

-}

 #ifdef _MSC_VER

 #pragma warning(disable: 4146) // unary minus operator applied to unsigned type

 #endif

@@ -152,34 +158,11 @@ STC_DEF stc64_uniform_t stc64_uniform_init(int64_t low, int64_t high) {
     dist.threshold = (uint64_t)(-dist.range) % dist.range;

     return dist;

 }

-

-STC_DEF stc32_uniform_t stc32_uniform_init(int32_t low, int32_t high) {

-    stc32_uniform_t dist = {low, (uint32_t) (high - low + 1)};

-    dist.threshold = (uint32_t)(-dist.range) % dist.range;

-    return dist;

-}

 #ifdef _MSC_VER

 #pragma warning(default: 4146)

 #endif

 

-STC_DEF uint64_t stc64_rand(stc64_t* rng) {

-    uint64_t *s = rng->state; enum {LR=24, RS=11, LS=3};

-    const uint64_t result = (s[0] ^ (s[3] += s[4])) + s[1];

-    s[0] = s[1] ^ (s[1] >> RS);

-    s[1] = s[2] + (s[2] << LS);

-    s[2] = ((s[2] << LR) | (s[2] >> (64 - LR))) + result;

-    return result;

-}

-

-STC_DEF uint32_t stc32_rand(stc32_t* rng) {

-    uint32_t *s = rng->state; enum {LR=21, RS=9, LS=3};

-    const uint32_t result = (s[0] ^ (s[3] += s[4])) + s[1];

-    s[0] = s[1] ^ (s[1] >> RS);

-    s[1] = s[2] + (s[2] << LS);

-    s[2] = ((s[2] << LR) | (s[2] >> (32 - LR))) + result;

-    return result;

-}

-

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

 STC_DEF int64_t stc64_uniform(stc64_t* rng, stc64_uniform_t* d) {

 #ifdef c_umul128

     uint64_t lo, hi;

@@ -195,29 +178,6 @@ STC_DEF int64_t stc64_uniform(stc64_t* rng, stc64_uniform_t* d) {
 #endif

 }

 

-/* Init uniform distributed float64 RNG, range [low, high). */

-STC_DEF stc64_uniformf_t stc64_uniformf_init(double low, double high) {

-    return c_make(stc64_uniformf_t){low, high - low};

-}

-

-/* Marsaglia polar method for gaussian/normal distribution, float64. */

-STC_DEF stc64_normalf_t stc64_normalf_init(double mean, double stddev) {

-    return c_make(stc64_normalf_t){mean, stddev, 0.0, 0};

-}

-

-STC_DEF double stc64_normalf(stc64_t* rng, stc64_normalf_t* dist) {

-    double u1, u2, s, m;

-    if (dist->has_next++ & 1)

-        return dist->next * dist->stddev + dist->mean;

-    do {

-        u1 = 2.0 * stc64_randf(rng) - 1.0;

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

-        s = u1*u1 + u2*u2;

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

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

-    dist->next = u2 * m;

-    return (u1 * m) * dist->stddev + dist->mean;

-}

 #endif

 #endif

 #undef i_opt