Mersenne Twister 64

5 files changed, 137 insertions, 2 deletions

diff --git a/source/kit/CMakeLists.txt b/source/kit/CMakeLists.txt
index 9cc1f7e..5891851 100644
--- a/source/kit/CMakeLists.txt
+++ b/source/kit/CMakeLists.txt
@@ -3,6 +3,7 @@ target_sources(
     PRIVATE
       atomic.c input_buffer.c input_stream.c lower_bound.c
       string_ref.c async_function.c allocator.c array_ref.c dynamic_array.c
+      mersenne_twister_64.c
     PUBLIC
       $<BUILD_INTERFACE:${CMAKE_CURRENT_LIST_DIR}/atomic.h>
       $<BUILD_INTERFACE:${CMAKE_CURRENT_LIST_DIR}/allocator.h>
@@ -12,4 +13,5 @@ target_sources(
       $<BUILD_INTERFACE:${CMAKE_CURRENT_LIST_DIR}/input_stream.h>
       $<BUILD_INTERFACE:${CMAKE_CURRENT_LIST_DIR}/input_buffer.h>
       $<BUILD_INTERFACE:${CMAKE_CURRENT_LIST_DIR}/lower_bound.h>
-      $<BUILD_INTERFACE:${CMAKE_CURRENT_LIST_DIR}/array_ref.h>)
+      $<BUILD_INTERFACE:${CMAKE_CURRENT_LIST_DIR}/array_ref.h>
+      $<BUILD_INTERFACE:${CMAKE_CURRENT_LIST_DIR}/mersenne_twister_64.h>)
diff --git a/source/kit/mersenne_twister_64.c b/source/kit/mersenne_twister_64.c
new file mode 100644
index 0000000..4be4313
--- /dev/null
+++ b/source/kit/mersenne_twister_64.c
@@ -0,0 +1,62 @@
+#include "mersenne_twister_64.h"
+
+#include <time.h>
+
+void kit_mt64_init(kit_mt64_state_t *state, uint64_t seed) {
+  state->mt[0] = seed;
+  for (state->index = 1; state->index < KIT_MT64_N; state->index++)
+    state->mt[state->index] = (6364136223846793005ull *
+                                   (state->mt[state->index - 1] ^
+                                    (state->mt[state->index - 1] >>
+                                     62u)) +
+                               state->index);
+}
+
+uint64_t kit_mt64_generate(kit_mt64_state_t *state) {
+  static uint64_t const mag01[2] = { 0ull, 0xB5026F5AA96619E9ull };
+
+  int      i;
+  uint64_t x;
+
+  if (state->index >= KIT_MT64_N) {
+    for (i = 0; i < KIT_MT64_N - 156; i++) {
+      x = (state->mt[i] & 0xFFFFFFFF80000000ull) |
+          (state->mt[i + 1] & 0x7FFFFFFFull);
+      state->mt[i] = state->mt[i + 156] ^ (x >> 1u) ^
+                     mag01[(int) (x & 1ull)];
+    }
+
+    for (; i < KIT_MT64_N - 1; i++) {
+      x = (state->mt[i] & 0xFFFFFFFF80000000ull) |
+          (state->mt[i + 1] & 0x7FFFFFFFull);
+      state->mt[i] = state->mt[i + (156 - KIT_MT64_N)] ^ (x >> 1u) ^
+                     mag01[(int) (x & 1ull)];
+    }
+
+    x = (state->mt[KIT_MT64_N - 1] & 0xFFFFFFFF80000000ull) |
+        (state->mt[0] & 0x7FFFFFFFull);
+    state->mt[KIT_MT64_N - 1] = state->mt[156 - 1] ^ (x >> 1u) ^
+                                mag01[(int) (x & 1ull)];
+
+    state->index = 0;
+  }
+
+  x = state->mt[state->index++];
+
+  x ^= (x >> 29u) & 0x5555555555555555ull;
+  x ^= (x << 17u) & 0x71D67FFFEDA60000ull;
+  x ^= (x << 37u) & 0xFFF7EEE000000000ull;
+  x ^= (x >> 43u);
+
+  return x;
+}
+
+uint64_t kit_mt64_seed() {
+  struct timespec t;
+  timespec_get(&t, TIME_UTC);
+
+  kit_mt64_state_t s;
+  kit_mt64_init(&s, (uint64_t) t.tv_nsec);
+
+  return kit_mt64_generate(&s);
+}
diff --git a/source/kit/mersenne_twister_64.h b/source/kit/mersenne_twister_64.h
new file mode 100644
index 0000000..605a25f
--- /dev/null
+++ b/source/kit/mersenne_twister_64.h
@@ -0,0 +1,37 @@
+#ifndef KIT_MERSENNE_TWISTER_64_H
+#define KIT_MERSENNE_TWISTER_64_H
+
+#include <stddef.h>
+#include <stdint.h>
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+enum {
+  KIT_MT64_N = 312,
+};
+
+typedef struct {
+  uint64_t mt[KIT_MT64_N];
+  uint64_t index;
+} kit_mt64_state_t;
+
+void kit_mt64_init(kit_mt64_state_t *state, uint64_t seed);
+
+uint64_t kit_mt64_generate(kit_mt64_state_t *state);
+
+uint64_t kit_mt64_seed();
+
+#ifndef LAPLACE_DISABLE_SHORT_NAMES
+#  define mt64_state_t kit_mt64_state_t
+#  define mt64_init kit_mt64_init
+#  define mt64_generate kit_mt64_generate
+#  define mt64_seed kit_mt64_seed
+#endif
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif
diff --git a/source/test/unittests/CMakeLists.txt b/source/test/unittests/CMakeLists.txt
index 3dda690..5ec950d 100644
--- a/source/test/unittests/CMakeLists.txt
+++ b/source/test/unittests/CMakeLists.txt
@@ -3,4 +3,5 @@ target_sources(
     PRIVATE
       async_function.test.c test_duration.test.c main.test.c
       string_ref.test.c atomic.test.c array_ref.test.c input_stream.test.c
-      lower_bound.test.c input_buffer.test.c dynamic_array.test.c)
+      lower_bound.test.c mersenne_twister_64.test.c input_buffer.test.c
+      dynamic_array.test.c)
diff --git a/source/test/unittests/mersenne_twister_64.test.c b/source/test/unittests/mersenne_twister_64.test.c
new file mode 100644
index 0000000..37c82b4
--- /dev/null
+++ b/source/test/unittests/mersenne_twister_64.test.c
@@ -0,0 +1,33 @@
+#include "../../kit/mersenne_twister_64.h"
+
+#define TEST_FILE mersenne_twister_64
+#include "../../kit_test/test.h"
+
+enum { SIZE = 1000 };
+
+TEST("mt64 same seeds") {
+  uint64_t seed = mt64_seed();
+
+  mt64_state_t foo, bar;
+  mt64_init(&foo, seed);
+  mt64_init(&bar, seed);
+
+  int ok = 1;
+  for (ptrdiff_t i = 0; i < SIZE; i++)
+    ok = ok && mt64_generate(&foo) == mt64_generate(&bar);
+
+  REQUIRE(ok);
+}
+
+TEST("mt64 different seeds") {
+  mt64_state_t foo, bar;
+  mt64_init(&foo, 42);
+  mt64_init(&bar, 4242424242);
+
+  ptrdiff_t difference_count = 0;
+  for (ptrdiff_t i = 0; i < SIZE; i++)
+    if (mt64_generate(&foo) != mt64_generate(&bar))
+      difference_count++;
+
+  REQUIRE(difference_count > SIZE / 2);
+}