Removed Mersenne Twister and xoroshiro functions. Added PGC32. Refactored carray.

5 files changed, 86 insertions, 220 deletions

diff --git a/examples/benchmark.c b/examples/benchmark.c
index b6e31a22..a8c4fb86 100644
--- a/examples/benchmark.c
+++ b/examples/benchmark.c
@@ -26,12 +26,9 @@ KHASH_MAP_INIT_INT64(ii, uint64_t)
 size_t seed = 1234;

 static const double maxLoadFactor = 0.77;

 

-sfc64_t rng;

+sfc64_random_t rng;

 #define SEED(s) rng = sfc64_seed(seed)

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

-//mt19937_t rng;

-//#define SEED(s) rng = mt19937_seed(s)

-//#define RAND(N) (mt19937_rand(&rng) & ((1 << N) - 1))

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

 

 

 #define CMAP_SETUP(tag, Key, Value) CHash_##tag map = chash_init; \

diff --git a/examples/rngtest.c b/examples/rngtest.c
index f13306b0..d88d1fab 100644
--- a/examples/rngtest.c
+++ b/examples/rngtest.c
@@ -14,45 +14,26 @@ int main(void)
     uint64_t v;
     printf("start\n");
 
-    mt19937_t state = mt19937_default();
-    uint32_t k = mt19937_rand(&state);
-    printf("%u - %g\n", k, c_randToFloat(k));
-
+    pcg32_random_t pcg = pcg32_seed(time(NULL), 1);
     before = clock(); \
     v = 0;
     for (size_t i=0; i<NN; i++) {
-        v += mt19937_rand(&state);
+        v += pcg32_random(&pcg) & 0xffffffff;
     }
     difference = clock() - before;
-    printf("my-mt: %.02f, %zu\n", (float) difference / CLOCKS_PER_SEC, v);
-
-#ifdef __cplusplus
-    std::mt19937 mt_rand;
-    before = clock(); \
-    v = 0;
-    for (size_t i=0; i<NN; i++) {
-        v += mt_rand();
-    }
-    difference = clock() - before;
-    printf("c++mt: %.02f, %zu\n", (float) difference / CLOCKS_PER_SEC, v);
-#endif
+    printf("pcg32: %.02f, %zu\n", (float) difference / CLOCKS_PER_SEC, v);
 
-    xoroshiro128ss_t xo = xoroshiro128ss_seed(1234);
+    sfc64_random_t sfc = sfc64_seed(time(NULL));
     before = clock(); \
     v = 0;
     for (size_t i=0; i<NN; i++) {
-        v += xoroshiro128ss_rand(&xo) & 0xffffffff;
-    }
-    difference = clock() - before;
-    printf("xoros: %.02f, %zu\n", (float) difference / CLOCKS_PER_SEC, v);
-
-
-    sfc64_t sfc = sfc64_seed(1234);
-    before = clock(); \
-    v = 0;
-    for (size_t i=0; i<NN; i++) {
-        v += sfc64_rand(&sfc) & 0xffffffff;
+        v += sfc64_random(&sfc) & 0xffffffff;
     }
     difference = clock() - before;
     printf("sfc64: %.02f, %zu\n", (float) difference / CLOCKS_PER_SEC, v);
+
+    for (int i=0; i<8; ++i) printf("%f ", sfc64_fRandom(&sfc));
+    puts("");
+    for (int i=0; i<8; ++i) printf("%f ", pcg32_fRandom(&pcg));
+    puts("");
 }
\ No newline at end of file
diff --git a/stc/carray.h b/stc/carray.h
index 1e2be0cd..b39bb35c 100644
--- a/stc/carray.h
+++ b/stc/carray.h
@@ -56,8 +56,8 @@ int main()
 #define carray2_size(a) ((a)._yxdim)

 #define carray3_size(a) _carray3_size(&(a)._zdim)

 

-#define _carray_own ((SIZE_MAX>>1)+1)

-#define _carray_sub (SIZE_MAX>>1)

+#define _carray_sub (SIZE_MAX >> 1)

+#define _carray_own (_carray_sub + 1)

 

 static inline size_t _carray_ydim(const size_t* yxdim) {

     return yxdim[0] / (yxdim[-1] & _carray_sub);

@@ -89,22 +89,6 @@ static inline size_t _carray3_size(const size_t* zdim) {
         size_t _xdim, _yxdim, _zdim; \

     } CArray3_##tag; \

  \

-    static inline void \

-    carray1_##tag##_destroy(CArray1_##tag* self) { \

-        size_t n = carray1_size(*self); Value* a = self->data; \

-        if (self->_xdim & _carray_own) {while (n--) valueDestroy(&a[n]); free(a);} \

-    } \

-    static inline void \

-    carray2_##tag##_destroy(CArray2_##tag* self) { \

-        size_t n = carray2_size(*self); Value* a = self->data; \

-        if (self->_xdim & _carray_own) {while (n--) valueDestroy(&a[n]); free(a);} \

-    } \

-    static inline void \

-    carray3_##tag##_destroy(CArray3_##tag* self) { \

-        size_t n = carray3_size(*self); Value* a = self->data; \

-        if (self->_xdim & _carray_own) {while (n--) valueDestroy(&a[n]); free(a);} \

-    } \

- \

     static inline CArray1_##tag \

     carray1_##tag##_make(size_t xdim, Value val) { \

         Value* m = c_new_N(Value, xdim); \

@@ -129,6 +113,28 @@ static inline size_t _carray3_size(const size_t* zdim) {
         return a; \

     } \

  \

+    static inline void \

+    carray1_##tag##_destroy(CArray1_##tag* self) { \

+        if (self->_xdim & _carray_own) { \

+            size_t n = carray1_size(*self); Value* a = self->data; \

+            while (n--) valueDestroy(&a[n]); free(a); \

+        } \

+    } \

+    static inline void \

+    carray2_##tag##_destroy(CArray2_##tag* self) { \

+        if (self->_xdim & _carray_own) { \

+            size_t n = carray2_size(*self); Value* a = self->data; \

+            while (n--) valueDestroy(&a[n]); free(a); \

+        } \

+    } \

+    static inline void \

+    carray3_##tag##_destroy(CArray3_##tag* self) { \

+        if (self->_xdim & _carray_own) { \

+            size_t n = carray3_size(*self); Value* a = self->data; \

+            while (n--) valueDestroy(&a[n]); free(a); \

+        } \

+    } \

+ \

     static inline CArray1_##tag \

     carray2_##tag##_at(CArray2_##tag a, size_t y) { \

         CArray1_##tag sub = {a.data + y*carray_xdim(a), carray_xdim(a)}; \

diff --git a/stc/cdefs.h b/stc/cdefs.h
index 6169229b..7f141a77 100644
--- a/stc/cdefs.h
+++ b/stc/cdefs.h
@@ -27,12 +27,13 @@
 #include <stddef.h>

 #include <stdbool.h>

 

-#if 0 // defined(_MSC_VER)

-#  define STC_INLINE static __forceinline

-#elif 0 // defined(__GNUC__)

-#  define STC_INLINE static inline __attribute((always_inline))

-#else // don't force

-#  define STC_INLINE static inline

+#define STC_INLINE static inline

+#if defined(_MSC_VER)

+#define STC_FORCE_INLINE static __forceinline

+#elif defined(__GNUC__)

+#define STC_FORCE_INLINE static inline __attribute((always_inline))

+#else

+#define STC_FORCE_INLINE static inline

 #endif

 

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

diff --git a/stc/crandom.h b/stc/crandom.h
index 1e8ec6ee..7cd22619 100644
--- a/stc/crandom.h
+++ b/stc/crandom.h
@@ -4,193 +4,74 @@
 #include "cdefs.h"
 #include <string.h>
 
-/*
- * Mersenne Twister random number generator MT19937, 32 bit.
- */
+typedef struct {uint64_t state;  uint64_t inc;} pcg32_random_t;
 
-enum {mt19937_N = 624, mt19937_M = 397};
-typedef struct {uint32_t idx, arr[mt19937_N];} mt19937_t;
-
-/* initializes state with a seed */
-STC_API void mt19937_init(mt19937_t *state, uint32_t seed) {
-    state->idx = mt19937_N;
-    state->arr[0] = seed;
-    for (int i = 1; i < mt19937_N; ++i) {
-        seed = state->arr[i] = 1812433253 * (seed ^ (seed >> 30)) + i; 
-    }
-}
-/* creates a new state from a seed */
-STC_API mt19937_t mt19937_seed(uint32_t seed) {
-    mt19937_t state;
-    mt19937_init(&state, seed);
-    return state;
-}
-/* creates a new state from default seed */
-STC_API mt19937_t mt19937_default(void) {
-    mt19937_t state;
-    mt19937_init(&state, 5489);
-    return state;
+/* 32 bit random number generator */
+STC_INLINE uint32_t pcg32_random(pcg32_random_t* rng)
+{
+    uint64_t old = rng->state;
+    rng->state = old * 6364136223846793005ull + rng->inc;
+    uint32_t xos = ((old >> 18u) ^ old) >> 27u;
+    uint32_t rot = old >> 59u;
+    return (xos >> rot) | (xos << ((-rot) & 31));
 }
 
-/* generates a random number in [0, 2^32)-interval */
-STC_API uint32_t mt19937_rand(mt19937_t *state) {
-    enum {N = mt19937_N, M = mt19937_M};
+/* float random int number in range [0, 1). NB: 23 bit resolution. */
+STC_INLINE float pcg32_fRandom(pcg32_random_t* rng) {
+    union {uint32_t i; float f;} u = {0x3F800000u | (pcg32_random(rng) >> 9)};
+    return u.f - 1.0f;
+}
 
-    uint32_t y, *arr = state->arr;
-    if (state->idx >= N) { /* generate N words at one time */
-        int k = 0;
-        for (; k < N-M; ++k) {
-            y = (arr[k] & 0x80000000) | (arr[k + 1] & 0x7fffffff);
-            arr[k] = arr[k + M] ^ (y >> 1) ^ ((y & 1) ? 0x9908b0df : 0);
-        }
-        for (; k < N-1; ++k) {
-            y = (arr[k] & 0x80000000) | (arr[k + 1] & 0x7fffffff);
-            arr[k] = arr[k - (N-M)] ^ (y >> 1) ^ ((y & 1) ? 0x9908b0df : 0);
-        }
-        y = (arr[N-1] & 0x80000000) | (arr[0] & 0x7fffffff);
-        arr[N-1] = arr[M-1] ^ (y >> 1) ^ ((y & 1) ? 0x9908b0df : 0);
+/* Uniform random number in range [0, bound) */
+STC_INLINE uint32_t pcg32_bRandom(pcg32_random_t* rng, uint32_t bound) {
+    return (uint32_t) (((uint64_t) pcg32_random(rng) * bound) >> 32);
+}
 
-        state->idx = 0;
-    }
-    /* tempering */
-    y = arr[state->idx++];
-    y ^= (y >> 11);
-    y ^= (y <<  7) & 0x9d2c5680;
-    y ^= (y << 15) & 0xefc60000;
-    y ^= (y >> 18);
-    return y;
+STC_INLINE pcg32_random_t pcg32_seed(uint64_t seed, uint64_t seq) {
+    pcg32_random_t rng = {0u, (seq << 1u) | 1u}; /* inc must be odd */
+    pcg32_random(&rng);
+    rng.state += seed;
+    pcg32_random(&rng);
+    return rng;
 }
 
 /* 
- * rotate bits left uint64
+ * Rotate bits left
  */
 STC_INLINE uint64_t c_rotateLeft64(uint64_t x, int bits) {
     return (x << bits) | (x >> (64 - bits));
 }
 
-/* 
- * xoroshiro128** with splitmix64 seed initialization
- */
-STC_API uint64_t splitmix64(uint64_t *state) {
-    uint64_t z = (*state += 0x9e3779b97f4a7c15);
-    z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9;
-    z = (z ^ (z >> 27)) * 0x94d049bb133111eb;
-    return z ^ (z >> 31);
-}
-
-typedef struct {uint64_t s[2];} xoroshiro128ss_t;
-
-STC_API xoroshiro128ss_t xoroshiro128ss_seed(uint64_t seed) {
-    xoroshiro128ss_t state;
-    state.s[0] = splitmix64(&seed);
-    state.s[1] = splitmix64(&seed);
-    return state;
-}
-
-STC_API uint64_t xoroshiro128ss_rand(xoroshiro128ss_t *state) {
-    uint64_t *s = state->s;
-    const uint64_t xo = s[0] ^ s[1];
-    const uint64_t result = c_rotateLeft64(s[0] * 5, 7) * 9;
-    s[0] = c_rotateLeft64(s[0], 24) ^ xo ^ (xo << 16);
-    s[1] = c_rotateLeft64(xo, 37);
-    return result;
-}
-
-STC_API void xoroshiro128ss_jump(xoroshiro128ss_t *state, bool longJump) {
-    static const uint64_t JUMP2_64[] = {0xdf900294d8f554a5, 0x170865df4b3201fc},
-                          JUMP2_96[] = {0xd2a98b26625eee7b, 0xdddf9b1090aa7ac1};
-    const uint64_t *jump = longJump ? JUMP2_96 : JUMP2_64;
-    uint64_t s0 = 0, s1 = 0, *s = state->s;
-
-    for (int i = 0; i < 2; ++i) for (int b = 0; b < 64; ++b) {
-        if (jump[i] & (1ull << b)) s0 ^= s[0], s1 ^= s[1];
-        xoroshiro128ss_rand(state);
-    }
-    s[0] = s0;
-    s[1] = s1;
-}
-
-/*
- * xoshiro256**:
- */
-typedef struct {uint64_t s[4];} xoshiro256ss_t;
-
-STC_API uint64_t xoshiro256ss_rand(xoshiro256ss_t* state) {
-    uint64_t *s = state->s;
-    const uint64_t result = c_rotateLeft64(s[1] * 5, 7) * 9;
-    const uint64_t t = s[1] << 17;
-    s[2] ^= s[0];
-    s[3] ^= s[1];
-    s[1] ^= s[2];
-    s[0] ^= s[3];
-    s[2] ^= t;
-    s[3] = c_rotateLeft64(s[3], 45);
-    return result;
-}
-
-STC_API xoshiro256ss_t xoshiro256ss_seed(uint64_t seed) {
-    xoshiro256ss_t state;
-    for (int i=0; i<4; ++i) state.s[i] = splitmix64(&seed);
-    return state;
-}
-
-/*
- * sfc32: http://pracrand.sourceforge.net
- */
-typedef struct {uint32_t s[4];} sfc32_t;
-
-STC_API uint32_t sfc32_rand(sfc32_t* state) {
-    enum {LR=21, RS=9, LS=3};
-    uint32_t *s = state->s;
-    const uint32_t out = s[0] + s[1] + s[3]++;
-    s[0] = s[1] ^ (s[1] >> RS);
-    s[1] = s[2] + (s[2] << LS);
-    s[2] = ((s[2] << LR) | (s[2] >> (32 - LR))) + out;
-    return out;
-}
-
-STC_API sfc32_t sfc32_seed(const uint64_t seed) {
-    sfc32_t state = {{0, (uint32_t) seed, (uint32_t) (seed >> 32), 1}};
-    for (int i=0; i<12; ++i) sfc32_rand(&state);
-    return state;
-}
-
 /*
  * sfc64: http://pracrand.sourceforge.net
  */
-typedef struct {uint64_t s[4];} sfc64_t;
+typedef struct {uint64_t state[4];} sfc64_random_t;
 
-STC_API uint64_t sfc64_rand(sfc64_t* state) {
+/* 64 bit random number generator */
+STC_API uint64_t sfc64_random(sfc64_random_t* rng) {
     enum {LR=24, RS=11, LS=3};
-    uint64_t *s = state->s;
-    const uint64_t out = s[0] + s[1] + s[3]++;
+    uint64_t *s = rng->state;
+    const uint64_t result = s[0] + s[1] + s[3]++;
     s[0] = s[1] ^ (s[1] >> RS);
     s[1] = s[2] + (s[2] << LS);
-    s[2] = ((s[2] << LR) | (s[2] >> (64 - LR))) + out;
-    return out;
+    s[2] = c_rotateLeft64(s[2], LR) + result;
+    return result;
 }
 
-STC_API sfc64_t sfc64_seed(const uint64_t seed) {
-    sfc64_t state = {{seed, seed, seed, 1}};
-    for (int i = 0; i < 12; ++i) sfc64_rand(&state);
-    return state;
+/* double random int number in range [0, 1). */
+STC_INLINE double sfc64_fRandom(sfc64_random_t* rng) {
+    union {uint64_t i; double f;} u = {0x3FF0000000000000ull | (sfc64_random(rng) >> 12)};
+    return u.f - 1.0;
 }
 
-/*
- * convert random int number to float in [0, 1) range.
- */
-STC_INLINE float c_randToFloat(uint32_t rnd) {
-    union {uint32_t i; float f;} u = {0x3F800000u | (rnd >> 9)};
-    return u.f - 1.0f;
-}
- 
-STC_INLINE double c_randToDouble(uint64_t rnd) {
-    union {uint64_t i; double f;} u = {0x3FF0000000000000ull | (rnd >> 12)};
-    return u.f - 1.0;
+STC_API sfc64_random_t sfc64_seed(const uint64_t seed) {
+    sfc64_random_t state = {{seed, seed, seed, 1}};
+    for (int i = 0; i < 12; ++i) sfc64_random(&state);
+    return state;
 }
 
 /*
- * siphash implementation.
+ * SipHash implementation.
  */
 #if defined(_WIN32) || (defined(__BYTE_ORDER__) && __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__)
     STC_INLINE uint64_t c_le64ToHost(uint64_t x) { return x; }
@@ -212,7 +93,7 @@ typedef struct siphash_t {
 } siphash_t;
 
 /* c=2, d=4 or c=1, d=3 */
-STC_API void siphash_init_c_d(siphash_t* s, const uint64_t key[2], const int c, const int d) {
+STC_INLINE void siphash_init_c_d(siphash_t* s, const uint64_t key[2], const int c, const int d) {
     s->c = c;
     s->d = d;
     s->length = 0;
@@ -303,7 +184,7 @@ STC_API uint64_t siphash_hash_c_d(const uint64_t key[2], const void* bytes, cons
 }
 
 /* default hash 2-4 */
-STC_API uint64_t siphash_hash(const uint64_t key[2], const void* bytes, const uint64_t size) {
+STC_INLINE uint64_t siphash_hash(const uint64_t key[2], const void* bytes, const uint64_t size) {
     return siphash_hash_c_d(key, bytes, size, 2, 4);
 }