Implemented mt19937_64 pseudo-random generator of 64-bit numbers

sparky-game · Dec 6, 2024 · 35c1724 · 35c1724
1 parent 0e1508b
commit 35c1724
Show file tree

Hide file tree

Showing 2 changed files with 46 additions and 0 deletions.
diff --git a/carbon.h b/carbon.h
@@ -310,6 +310,8 @@ typedef struct {
 CARBON_API void carbon_math_srand(u64 seed);
 CARBON_API int carbon_math_rand(void);
 CARBON_API f32 carbon_math_randf(void);
+CARBON_API void carbon_math_mt19937_64_srand(u64 seed);
+CARBON_API u64 carbon_math_mt19937_64_rand(void);
 CARBON_API f32 carbon_math_abs(f32 x);
 CARBON_API f32 carbon_math_round(f32 x);
 CARBON_API f32 carbon_math_floor(f32 x);

diff --git a/src/carbon_math.c b/src/carbon_math.c
@@ -5,7 +5,15 @@
 #include <carbon.h>
 #endif  // CARBON_IMPLEMENTATION
 
+#define CARBON_MATH_MT19937_64_RAND_NN 312
+#define CARBON_MATH_MT19937_64_RAND_MM 156
+#define CARBON_MATH_MT19937_64_RAND_MATRIX_A 0xB5026F5AA96619E9ULL
+#define CARBON_MATH_MT19937_64_RAND_UM 0xFFFFFFFF80000000ULL
+#define CARBON_MATH_MT19937_64_RAND_LM 0x7FFFFFFFULL
+
 static u64 carbon_math__rand_seed;
+static u64 carbon_math__mt19937_64_rand_state_vec[CARBON_MATH_MT19937_64_RAND_NN];
+static u64 carbon_math__mt19937_64_rand_state_idx = CARBON_MATH_MT19937_64_RAND_NN + 1;
 
 void carbon_math_srand(u64 seed) {
   carbon_math__rand_seed = seed - 1;
@@ -20,6 +28,42 @@ f32 carbon_math_randf(void) {
   return (f32) carbon_math_rand() / (f32) CARBON_RAND_MAX;
 }
 
+void carbon_math_mt19937_64_srand(u64 seed) {
+  u64 *sv = carbon_math__mt19937_64_rand_state_vec;
+  u64 *si = &carbon_math__mt19937_64_rand_state_idx;
+  sv[0] = seed;
+  for (*si = 1; *si < CARBON_MATH_MT19937_64_RAND_NN; ++(*si)) {
+    sv[*si] = (6364136223846793005ULL * (sv[*si - 1] ^ (sv[*si - 1] >> 62)) + *si);
+  }
+}
+
+u64 carbon_math_mt19937_64_rand(void) {
+  static u64 mag[] = {0ULL, CARBON_MATH_MT19937_64_RAND_MATRIX_A};
+  i32 i; u64 x;
+  u64 *sv = carbon_math__mt19937_64_rand_state_vec;
+  u64 *si = &carbon_math__mt19937_64_rand_state_idx;
+  if (*si >= CARBON_MATH_MT19937_64_RAND_NN) {
+    if (*si == CARBON_MATH_MT19937_64_RAND_NN + 1) carbon_math_mt19937_64_srand(5489ULL);
+    for (i = 0; i < CARBON_MATH_MT19937_64_RAND_NN - CARBON_MATH_MT19937_64_RAND_MM; ++i) {
+      x = (sv[i] & CARBON_MATH_MT19937_64_RAND_UM) | (sv[i + 1] & CARBON_MATH_MT19937_64_RAND_LM);
+      sv[i] = sv[i + CARBON_MATH_MT19937_64_RAND_MM] ^ (x >> 1) ^ mag[(i32) (x & 1ULL)];
+    }
+    for (; i < CARBON_MATH_MT19937_64_RAND_NN - 1; ++i) {
+      x = (sv[i] & CARBON_MATH_MT19937_64_RAND_UM) | (sv[i + 1] & CARBON_MATH_MT19937_64_RAND_LM);
+      sv[i] = sv[i + (CARBON_MATH_MT19937_64_RAND_MM - CARBON_MATH_MT19937_64_RAND_NN)] ^ (x >> 1) ^ mag[(i32) (x & 1ULL)];
+    }
+    x = (sv[CARBON_MATH_MT19937_64_RAND_NN - 1] & CARBON_MATH_MT19937_64_RAND_UM) | (sv[0] & CARBON_MATH_MT19937_64_RAND_LM);
+    sv[CARBON_MATH_MT19937_64_RAND_NN - 1] = sv[CARBON_MATH_MT19937_64_RAND_MM - 1] ^ (x >> 1) ^ mag[(i32) (x & 1ULL)];
+    *si = 0;
+  }
+  x = sv[(*si)++];
+  x ^= (x >> 29) & 0x5555555555555555ULL;
+  x ^= (x << 17) & 0x71D67FFFEDA60000ULL;
+  x ^= (x << 37) & 0xFFF7EEE000000000ULL;
+  x ^= (x >> 43);
+  return x;
+}
+
 f32 carbon_math_abs(f32 x) {
   union { f32 f; u32 i; } u = {x};
   u.i &= CARBON_I32_MAX;