summaryrefslogtreecommitdiff
path: root/tests/rand.c
diff options
context:
space:
mode:
Diffstat (limited to 'tests/rand.c')
-rw-r--r--tests/rand.c83
1 files changed, 83 insertions, 0 deletions
diff --git a/tests/rand.c b/tests/rand.c
new file mode 100644
index 0000000..925d1ba
--- /dev/null
+++ b/tests/rand.c
@@ -0,0 +1,83 @@
+// Taken from
+// http://www.math.sci.hiroshima-u.ac.jp/m-mat/MT/MT2002/CODES/mt19937ar.c
+
+#define _XOPEN_SOURCE 700
+
+#include <assert.h>
+#include <inttypes.h>
+#include <stdbool.h>
+#include <stdio.h>
+
+
+#define N 624U
+static const uint16_t M = 397;
+static const uint32_t MATRIX_A = 0x9908b0dfUL;
+static const uint32_t UPPER_MASK = 0x80000000UL;
+static const uint32_t LOWER_MASK = 0x7fffffffUL;
+
+static uint32_t mt[N];
+static uint16_t mti;
+static bool initialized = false;
+
+static void
+init(uint32_t s) {
+ initialized = true;
+ mt[0] = s & 0xffffffffUL;
+ for (mti = 1; mti < N; mti++) {
+ mt[mti] = (1812433253 * (mt[mti - 1] ^ (mt[mti - 1] >> 30))) + mti;
+ }
+}
+
+/// Generates a random number on the [ 0, 0xffffffff ]-interval.
+static uint32_t
+generate(void) {
+ assert(initialized == true);
+ static uint32_t mag01[2] = { 0, MATRIX_A };
+
+ if (mti >= N) {
+ uint16_t kk;
+
+ for (kk = 0; kk < N - M; kk++) {
+ const uint32_t y0 =
+ (mt[kk] & UPPER_MASK) |
+ (mt[kk + 1] & LOWER_MASK);
+ mt[kk] = mt[kk + M] ^ (y0 >> 1) ^ mag01[y0 & 1U];
+ }
+
+ for (; kk < N - 1; kk++) {
+ const uint32_t y1 =
+ (mt[kk] & UPPER_MASK) |
+ (mt[kk + 1] & LOWER_MASK);
+ mt[kk] = mt[kk + (M - N)] ^ (y1 >> 1) ^ mag01[y1 & 1U];
+ }
+
+ const uint32_t y2 =
+ (mt[N - 1] & UPPER_MASK) |
+ (mt[0] & LOWER_MASK);
+ mt[N - 1] = mt[M - 1] ^ (y2 >> 1) ^ mag01[y2 & 1U];
+ mti = 0;
+ }
+
+ uint32_t y = mt[mti++];
+
+ y ^= (y >> 11);
+ y ^= (y << 7) & 0x9d2c5680UL;
+ y ^= (y << 15) & 0xefc60000UL;
+ y ^= (y >> 18);
+
+ return y;
+}
+
+
+int
+main(void) {
+ init(2222);
+
+ for (uint16_t i = 0; i < 10000; i++) {
+ assert(
+ printf("%" PRIu32 "\n", generate()) > 0
+ );
+ }
+
+ return 0;
+}