Some refactoring.

1 files changed, 47 insertions, 25 deletions

diff --git a/include/stc/crandom.h b/include/stc/crandom.h
index b6bd0c93..3e911d93 100644
--- a/include/stc/crandom.h
+++ b/include/stc/crandom.h
@@ -49,10 +49,10 @@ typedef struct {double lower, range;}                      stc64_uniformf_t;
 typedef struct {double mean, stddev, next; unsigned has_next;} stc64_normalf_t;

 

 /* PRNG stc64.

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

- * 320bit state, 256bit mutable.

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

+ * 320-bit state, 256 bit is mutable.

  * Noticable faster than xoshiro and pcg, but slighly slower than wyrand64, 

- * which only has a single 64-bit state and is unfit when longer periods or

+ * which only has a single 64-bit state and therefore unfit when

  * multiple threads are required.

  * stc64 supports 2^63 unique threads with a minimum 2^64 period lengths each.

  * Passes PractRand tested up to 8TB output, Vigna's Hamming weight test,

@@ -60,36 +60,37 @@ typedef struct {double mean, stddev, next; unsigned has_next;} stc64_normalf_t;
  * The 16-bit version with LR=6, RS=5, LS=3 passes PractRand to multiple TB input.

  */

 

-STC_API stc64_t stc64_init(uint64_t seed);

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

+/* Global STC64 PRNG */

+STC_API void            stc64_srandom(uint64_t seed);

+STC_API uint64_t        stc64_random(void);

+/* STC64 PRNG state */ 

+STC_API stc64_t         stc64_init(uint64_t seed);

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

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

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

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

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

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

+

 

 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]|1)) + s[1];

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

 }

 

-/* Global random() */

-static stc64_t stc64_global = {{0x26aa069ea2fb1a4d, 0x70c72c95cd592d04, 0x504f333d3aa0b359, 0x26aa069ea2fb1a4d, 0x6a09e667a754166b}};

-STC_INLINE void stc64_srandom(uint64_t seed) { stc64_global = stc64_init(seed); }

-STC_INLINE uint64_t stc64_random(void) { return stc64_rand(&stc64_global); }

-

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

 STC_INLINE double stc64_randf(stc64_t* rng) {

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

     return u.f - 1.0;

 }

 

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

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

-

 /* Float64 uniform distributed 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};

-}

 STC_INLINE double stc64_uniformf(stc64_t* rng, stc64_uniformf_t* dist) {

     return stc64_randf(rng)*dist->range + dist->lower;

 }

@@ -105,32 +106,53 @@ STC_INLINE int64_t stc64_uniform(stc64_t* rng, stc64_uniform_t* d) {
     return d->lower + hi;

 }

 

-/* Normal distributed RNG, Float64. */

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

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

-}

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

-

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

 

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

 

+/* Global random() */

+static stc64_t stc64_global = {{

+    0x26aa069ea2fb1a4d, 0x70c72c95cd592d04,

+    0x504f333d3aa0b359, 0x9e3779b97f4a7c15,

+    0x6a09e667a754166b

+}};

+

+STC_DEF void stc64_srandom(uint64_t seed) { 

+    stc64_global = stc64_init(seed); 

+}

+

+STC_DEF uint64_t stc64_random(void) {

+    return stc64_rand(&stc64_global);

+}

+

 STC_DEF stc64_t stc64_init(uint64_t seed) {

     return stc64_with_seq(seed, seed + 0x3504f333d3aa0b34);

 }

+

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

 STC_DEF stc64_t stc64_with_seq(uint64_t seed, uint64_t seq) {

     stc64_t rng = {{seed, seed, seed, seed, (seq << 1) | 1}};

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

     return rng;

 }

 

-/* Very fast unbiased uniform int RNG with bounds [low, high] */

+/* Init unbiased uniform uint64 RNG with bounds [low, high] */

 STC_DEF stc64_uniform_t stc64_uniform_init(int64_t low, int64_t high) {

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

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

     return dist;

 }

 

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

+/* 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)