1 files changed, 145 insertions, 0 deletions
diff --git a/test/MTwikipedia.c b/test/MTwikipedia.c
new file mode 100644
index 0000000..83eda43
--- /dev/null
+++ b/test/MTwikipedia.c
@@ -0,0 +1,145 @@
+// This file is the C Code from the Mersenne Twister wikipedia article
+//		https://en.wikipedia.org/wiki/Mersenne_Twister#C_code
+// bodged together with the init_by_array and main functions from
+//		https://www.math.sci.hiroshima-u.ac.jp/m-mat/MT/MT2002/CODES/mt19937ar.c
+// so I can test if the C Code on wikipedia generates the correct random numbers
+// as seen in
+//		https://www.math.sci.hiroshima-u.ac.jp/m-mat/MT/MT2002/CODES/mt19937ar.out
+// 
+// I did not check them all, but a few at random. They match, which makes me think
+// they operate the same. Now I need to look for the bug in my Rust,,,
+#include <stdint.h>
+#include <stdio.h>
+
+#define n 624
+#define m 397
+#define w 32
+#define r 31
+#define UMASK (0xffffffffUL << r)
+#define LMASK (0xffffffffUL >> (w-r))
+#define a 0x9908b0dfUL
+#define u 11
+#define s 7
+#define t 15
+#define l 18
+#define b 0x9d2c5680UL
+#define c 0xefc60000UL
+#define f 1812433253UL
+
+typedef struct
+{
+    uint32_t state_array[n];         // the array for the state vector 
+    int state_index;                 // index into state vector array, 0 <= state_index <= n-1   always
+} mt_state;
+
+
+void initialize_state(mt_state* state, uint32_t seed) 
+{
+    uint32_t* state_array = &(state->state_array[0]);
+    
+    state_array[0] = seed;                          // suggested initial seed = 19650218UL
+    
+    for (int i=1; i<n; i++)
+    {
+        seed = f * (seed ^ (seed >> (w-2))) + i;    // Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier.
+        state_array[i] = seed; 
+    }
+    
+    state->state_index = 0;
+}
+
+/* initialize by an array with array-length */
+/* init_key is the array for initializing keys */
+/* key_length is its length */
+/* slight change for C++, 2004/2/26 */
+void init_by_array(mt_state* state, unsigned long init_key[], int key_length)
+{
+    int i, j, k;
+    initialize_state(state, 19650218UL);
+    i=1; j=0;
+
+    k = (n>key_length ? n : key_length);
+
+	uint32_t* mt = &(state->state_array[0]);
+
+    for (; k; k--) {
+        mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525UL))
+          + init_key[j] + j; /* non linear */
+        mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
+        i++; j++;
+        if (i>=n) { mt[0] = mt[n-1]; i=1; }
+        if (j>=key_length) j=0;
+    }
+    for (k=n-1; k; k--) {
+        mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941UL))
+          - i; /* non linear */
+        mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
+        i++;
+        if (i>=n) { mt[0] = mt[n-1]; i=1; }
+    }
+
+    mt[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */ 
+}
+
+
+uint32_t random_uint32(mt_state* state)
+{
+    uint32_t* state_array = &(state->state_array[0]);
+    
+    int k = state->state_index;      // point to current state location
+                                     // 0 <= state_index <= n-1   always
+    
+//  int k = k - n;                   // point to state n iterations before
+//  if (k < 0) k += n;               // modulo n circular indexing
+                                     // the previous 2 lines actually do nothing
+                                     //  for illustration only
+    
+    int j = k - (n-1);               // point to state n-1 iterations before
+    if (j < 0) j += n;               // modulo n circular indexing
+
+    uint32_t x = (state_array[k] & UMASK) | (state_array[j] & LMASK);
+    
+    uint32_t xA = x >> 1;
+    if (x & 0x00000001UL) xA ^= a;
+    
+    j = k - (n-m);                   // point to state n-m iterations before
+    if (j < 0) j += n;               // modulo n circular indexing
+    
+    x = state_array[j] ^ xA;         // compute next value in the state
+    state_array[k++] = x;            // update new state value
+    
+    if (k >= n) k = 0;               // modulo n circular indexing
+    state->state_index = k;
+    
+    uint32_t y = x ^ (x >> u);       // tempering 
+             y = y ^ ((y << s) & b);
+             y = y ^ ((y << t) & c);
+    uint32_t z = y ^ (y >> l);
+    
+    return z; 
+}
+
+int main(void)
+{
+    int i;
+    unsigned long init[4]={0x123, 0x234, 0x345, 0x456}, length=4;
+
+	mt_state state;
+	state.state_index = 0;
+
+    initialize_state(&state, 19650218UL);
+    //init_by_array(&state, init, length);
+
+    /*printf("1000 outputs of genrand_int32()\n");
+    for (i=0; i<1000; i++) {
+      printf("%10lu ", random_uint32(&state));
+      if (i%5==4) printf("\n");
+    }*/
+
+    // Print the internal state so we can check it against the Rust
+    for (i = 0; i < n; ++i) {
+        printf("%10lu\n", state.state_array[i]);
+    }
+
+    return 0;
+}
\ No newline at end of file