summaryrefslogtreecommitdiff
path: root/util/mersenne_twister.c
diff options
context:
space:
mode:
Diffstat (limited to 'util/mersenne_twister.c')
-rw-r--r--util/mersenne_twister.c68
1 files changed, 68 insertions, 0 deletions
diff --git a/util/mersenne_twister.c b/util/mersenne_twister.c
new file mode 100644
index 0000000..2d3059f
--- /dev/null
+++ b/util/mersenne_twister.c
@@ -0,0 +1,68 @@
+/*
+ * base on pseudocode from
+ * https://en.wikipedia.org/wiki/Mersenne_Twister
+ * MT19937 version (32-bit word)
+ */
+
+#include "mersenne_twister.h"
+
+//MT-32 algorithm coefficients
+#define W 32UL
+#define N 624UL
+#define M 397UL
+#define R 31UL
+#define A 0x9908B0DFUL
+#define U 11UL
+#define D 0xFFFFFFFFUL
+#define S 7UL
+#define B 0x9D2C5680UL
+#define T 15UL
+#define C 0xEFC60000UL
+#define L 18UL
+#define F 1812433253UL //initialization constant
+
+//derived constants
+#define LOWER_MASK ((1UL << R) - 1)
+#define UPPER_MASK (~LOWER_MASK)
+
+//global state
+static uint32_t MT[N], idx = N + 1; //essentially unseeded
+
+//generate next 'N' values
+static void twist() {
+ for (uint32_t i = 0; i < N; i++) {
+ uint32_t x = (MT[i] & UPPER_MASK) + (MT[(i + 1) % N] & LOWER_MASK);
+ uint32_t x_a = x >> 1;
+ if ((x % 2) != 0)
+ x_a ^= A;
+ MT[i] = MT[(i + M) % N] ^ x_a;
+ }
+ idx = 0;
+}
+
+void set_state(uint32_t *new_state) {
+ idx = N;
+ for (uint32_t i = 0; i < N; i++)
+ MT[i] = *new_state++;
+}
+
+void seed(uint32_t seed) {
+ idx = N;
+ MT[0] = seed;
+ for (uint32_t i = 1; i < N; i++)
+ MT[i] = F * (MT[i - 1] ^ (MT[i - 1] >> (W - 2))) + i;
+}
+
+uint32_t extract() {
+ if (idx >= N) {
+ if (idx > N)
+ seed(5489); //default seed
+ twist();
+ }
+ uint32_t y = MT[idx++];
+ y ^= (y >> U) & D;
+ y ^= (y << S) & B;
+ y ^= (y << T) & C;
+ y ^= (y >> L);
+ return y;
+}