Added rngtest.c and included xoroshiro128ss_rand() and sfc64_rand()

2 files changed, 208 insertions, 6 deletions

diff --git a/rngtest.c b/rngtest.c
new file mode 100644
index 00000000..8987ff5d
--- /dev/null
+++ b/rngtest.c
@@ -0,0 +1,142 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <time.h>
+#include "stc/crandom.h"
+#ifdef __cplusplus
+#include <random>
+#endif
+
+
+
+/* Period parameters */  
+#define N 624
+#define M 397
+#define MATRIX_A 0x9908b0dfUL   /* constant vector a */
+#define UPPER_MASK 0x80000000UL /* most significant w-r bits */
+#define LOWER_MASK 0x7fffffffUL /* least significant r bits */
+
+static unsigned long mt[N]; /* the array for the state vector  */
+static int mti=N+1; /* mti==N+1 means mt[N] is not initialized */
+
+/* initializes mt[N] with a seed */
+void init_genrand(unsigned long s)
+{
+    mt[0]= s & 0xffffffffUL;
+    for (mti=1; mti<N; mti++) {
+        mt[mti] = 
+        (1812433253UL * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti); 
+        /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
+        /* In the previous versions, MSBs of the seed affect   */
+        /* only MSBs of the array mt[].                        */
+        /* 2002/01/09 modified by Makoto Matsumoto             */
+        mt[mti] &= 0xffffffffUL;
+        /* for >32 bit machines */
+    }
+}
+
+/* generates a random number on [0,0xffffffff]-interval */
+unsigned long genrand_int32(void)
+{
+    unsigned long y;
+    static unsigned long mag01[2]={0x0UL, MATRIX_A};
+    /* mag01[x] = x * MATRIX_A  for x=0,1 */
+
+    if (mti >= N) { /* generate N words at one time */
+        int kk;
+
+        if (mti == N+1)   /* if init_genrand() has not been called, */
+            init_genrand(5489UL); /* a default initial seed is used */
+
+        for (kk=0;kk<N-M;kk++) {
+            y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
+            mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1UL];
+        }
+        for (;kk<N-1;kk++) {
+            y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
+            mt[kk] = mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1UL];
+        }
+        y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK);
+        mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1UL];
+
+        mti = 0;
+    }
+  
+    y = mt[mti++];
+
+    /* Tempering */
+    y ^= (y >> 11);
+    y ^= (y << 7) & 0x9d2c5680UL;
+    y ^= (y << 15) & 0xefc60000UL;
+    y ^= (y >> 18);
+
+    return y;
+}
+
+
+
+#define NN 1000000000
+
+int main(void)
+{
+    clock_t difference, before;
+    uint64_t v;
+    printf("start\n");
+
+    before = clock(); \
+    v = 0;
+    for (size_t i=0; i<NN; i++) {
+        v += genrand_int32();
+    }
+    difference = clock() - before;
+    printf("refmt: %.02f, %llu\n", (float) difference / CLOCKS_PER_SEC, v);
+
+
+    mt19937_t state = mt19937_default();
+
+    before = clock(); \
+    v = 0;
+    for (size_t i=0; i<NN; i++) {
+        v += mt19937_rand(&state);
+    }
+    difference = clock() - before;
+    printf("mymt : %.02f, %llu\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, %llu\n", (float) difference / CLOCKS_PER_SEC, v);
+#endif
+
+    xoroshiro128ss_t xo = xoroshiro128ss_seed(1234);
+    before = clock(); \
+    v = 0;
+    for (size_t i=0; i<NN; i++) {
+        v += xoroshiro128ss_rand(&xo) & 0xffffffff;
+    }
+    difference = clock() - before;
+    printf("xoros: %.02f, %llu\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;
+    }
+    difference = clock() - before;
+    printf("sfc64: %.02f, %llu\n", (float) difference / CLOCKS_PER_SEC, v);
+
+
+    before = clock(); \
+    v = 0;
+    for (size_t i=0; i<NN; i++) {
+        v += rand();
+    }
+    difference = clock() - before;
+    printf("rand : %.02f, %llu\n", (float) difference / CLOCKS_PER_SEC, v);
+}
\ No newline at end of file
diff --git a/stc/crandom.h b/stc/crandom.h
index 3b3286aa..064c86eb 100644
--- a/stc/crandom.h
+++ b/stc/crandom.h
@@ -52,12 +52,12 @@ enum { /* period parameters */
 };
 
 typedef struct mt19937 {
-    uint32_t idx;
-    uint32_t arr[mt19937_N]; /* the array for the state vector */
+    uint_fast32_t idx;
+    uint_fast32_t arr[mt19937_N]; /* the array for the state vector */
 } mt19937_t;
 
 /* initializes state with a seed */
-static inline void mt19937_init(mt19937_t *state, uint32_t seed) {
+static inline void mt19937_init(mt19937_t *state, uint_fast32_t seed) {
     state->idx = mt19937_N;
     state->arr[0] = seed;
     for (int i = 1; i < mt19937_N; ++i) {
@@ -69,7 +69,7 @@ static inline void mt19937_init(mt19937_t *state, uint32_t seed) {
     }
 }
 /* creates a new state from a seed */
-static inline mt19937_t mt19937_seed(uint32_t seed) {
+static inline mt19937_t mt19937_seed(uint_fast32_t seed) {
     mt19937_t state;
     mt19937_init(&state, seed);
     return state;
@@ -82,10 +82,10 @@ static inline mt19937_t mt19937_default(void) {
 }
 
 /* generates a random number on [0, 0xffffffff]-interval */
-static inline uint32_t mt19937_rand(mt19937_t *state) {
+static inline uint_fast32_t mt19937_rand(mt19937_t *state) {
     enum {N = mt19937_N, M = mt19937_M};
 
-    uint32_t y, *arr = state->arr;
+    uint_fast32_t y, *arr = state->arr;
     if (state->idx >= N) { /* generate N words at one time */
         int k = 0;
         for (; k < N-M; ++k) {
@@ -110,4 +110,64 @@ static inline uint32_t mt19937_rand(mt19937_t *state) {
     return y;
 }
 
+/* xoroshiro128** and xoshiro128** with splitmix64 seed initialization */
+
+static inline 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);
+}
+
+static inline uint64_t c_rotl(uint64_t x, int s) {
+    return (x << s) | (x >> (64 - s));
+}
+
+/* xoroshiro128** */
+
+typedef struct xoroshiro128ss {
+    uint64_t s[2];
+} xoroshiro128ss_t;
+
+static inline xoroshiro128ss_t xoroshiro128ss_seed(uint64_t seed) {
+    xoroshiro128ss_t state;
+    state.s[0] = splitmix64(&seed);
+    state.s[1] = splitmix64(&seed);
+    return state;
+}
+
+static inline uint64_t xoroshiro128ss_rand(xoroshiro128ss_t *state) {
+    uint64_t *s = state->s, s1 = s[1];
+    const uint64_t s0 = s[0];
+    const uint64_t result = c_rotl(s0 * 5, 7) * 9;
+    s1 ^= s0;
+    s[0] = c_rotl(s0, 24) ^ s1 ^ (s1 << 16);
+    s[1] = c_rotl(s1, 37);
+    return result;
+}
+
+/* http://pracrand.sourceforge.net */
+
+typedef struct sfc64 {
+    uint64_t s[3], counter;
+} sfc64_t;
+
+
+static inline uint64_t sfc64_rand(sfc64_t* state) {
+    enum {LROT = 24, RSHIFT = 11, LSHIFT = 3};
+    uint64_t *s = state->s;
+
+    uint64_t result = s[0] + s[1] + ++state->counter;
+    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;
+}
+
+static inline sfc64_t sfc64_seed(uint64_t seed) {
+    sfc64_t state = {{seed, seed, seed}, 0};
+    for (int i = 0; i < 12; ++i) sfc64_rand(&state);
+    return state;
+}
+
 #endif