BENCH: convert bencmarks to asv format

remove files that were part of the origal repo
rework randomgen docs to integrate with numpy and fix some links
remove convenience functions, require explicit call to gen.brng
move code out of numpy.random.randomgen into numpy.random

10 files changed, 1574 insertions, 0 deletions

diff --git a/numpy/random/src/mt19937/LICENSE.md b/numpy/random/src/mt19937/LICENSE.md
new file mode 100644
index 000000000..f65c3d46e
--- /dev/null
+++ b/numpy/random/src/mt19937/LICENSE.md
@@ -0,0 +1,61 @@
+# MT19937
+
+Copyright (c) 2003-2005, Jean-Sebastien Roy (js@jeannot.org)
+
+The rk_random and rk_seed functions algorithms and the original design of
+the Mersenne Twister RNG:
+
+  Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
+  All rights reserved.
+
+  Redistribution and use in source and binary forms, with or without
+  modification, are permitted provided that the following conditions
+  are met:
+
+  1. Redistributions of source code must retain the above copyright
+  notice, this list of conditions and the following disclaimer.
+
+  2. Redistributions in binary form must reproduce the above copyright
+  notice, this list of conditions and the following disclaimer in the
+  documentation and/or other materials provided with the distribution.
+
+  3. The names of its contributors may not be used to endorse or promote
+  products derived from this software without specific prior written
+  permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER
+OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+Original algorithm for the implementation of rk_interval function from
+Richard J. Wagner's implementation of the Mersenne Twister RNG, optimised by
+Magnus Jonsson.
+
+Constants used in the rk_double implementation by Isaku Wada.
+
+Permission is hereby granted, free of charge, to any person obtaining a
+copy of this software and associated documentation files (the
+"Software"), to deal in the Software without restriction, including
+without limitation the rights to use, copy, modify, merge, publish,
+distribute, sublicense, and/or sell copies of the Software, and to
+permit persons to whom the Software is furnished to do so, subject to
+the following conditions:
+
+The above copyright notice and this permission notice shall be included
+in all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
+CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
+TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
+SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
\ No newline at end of file
diff --git a/numpy/random/src/mt19937/mt19937-benchmark.c b/numpy/random/src/mt19937/mt19937-benchmark.c
new file mode 100644
index 000000000..039f8030a
--- /dev/null
+++ b/numpy/random/src/mt19937/mt19937-benchmark.c
@@ -0,0 +1,31 @@
+/*
+ * cl mt19937-benchmark.c mt19937.c /Ox
+ * Measure-Command { .\mt19937-benchmark.exe }
+ *
+ * gcc mt19937-benchmark.c mt19937.c -O3 -o mt19937-benchmark
+ * time ./mt19937-benchmark
+ */
+#include "mt19937.h"
+#include <inttypes.h>
+#include <stdio.h>
+#include <time.h>
+
+#define Q 1000000000
+
+int main() {
+  int i;
+  uint32_t seed = 0x0;
+  uint64_t sum = 0, count = 0;
+  mt19937_state state;
+  mt19937_seed(&state, seed);
+  clock_t begin = clock();
+  for (i = 0; i < Q; i++) {
+    sum += mt19937_next64(&state);
+    count++;
+  }
+  clock_t end = clock();
+  double time_spent = (double)(end - begin) / CLOCKS_PER_SEC;
+  printf("0x%" PRIx64 "\ncount: %" PRIu64 "\n", sum, count);
+  printf("%" PRIu64 " randoms per second\n",
+         (uint64_t)(Q / time_spent) / 1000000 * 1000000);
+}
diff --git a/numpy/random/src/mt19937/mt19937-jump.c b/numpy/random/src/mt19937/mt19937-jump.c
new file mode 100644
index 000000000..46b28cf96
--- /dev/null
+++ b/numpy/random/src/mt19937/mt19937-jump.c
@@ -0,0 +1,224 @@
+#include "mt19937-jump.h"
+#include "mt19937.h"
+
+/* 32-bits function */
+/* return the i-th coefficient of the polynomial pf */
+unsigned long get_coef(unsigned long *pf, unsigned int deg) {
+  if ((pf[deg >> 5] & (LSB << (deg & 0x1ful))) != 0)
+    return (1);
+  else
+    return (0);
+}
+
+/* 32-bit function */
+/* set the coefficient of the polynomial pf with v */
+void set_coef(unsigned long *pf, unsigned int deg, unsigned long v) {
+  if (v != 0)
+    pf[deg >> 5] ^= (LSB << (deg & 0x1ful));
+  else
+    ;
+}
+
+void gray_code(unsigned long *h) {
+  unsigned int i, j = 1, l = 1, term = LL;
+
+  h[0] = 0;
+
+  for (i = 1; i <= QQ; i++) {
+    l = (l << 1);
+    term = (term >> 1);
+    for (; j < l; j++)
+      h[j] = h[l - j - 1] ^ term;
+  }
+}
+
+void copy_state(mt19937_state *target_state, mt19937_state *state) {
+  int i;
+
+  for (i = 0; i < N; i++)
+    target_state->key[i] = state->key[i];
+
+  target_state->pos = state->pos;
+}
+
+/* next state generating function */
+void gen_next(mt19937_state *state) {
+  int num;
+  unsigned long y;
+  static unsigned long mag02[2] = {0x0ul, MATRIX_A};
+
+  num = state->pos;
+  if (num < N - M) {
+    y = (state->key[num] & UPPER_MASK) | (state->key[num + 1] & LOWER_MASK);
+    state->key[num] = state->key[num + M] ^ (y >> 1) ^ mag02[y % 2];
+    state->pos++;
+  } else if (num < N - 1) {
+    y = (state->key[num] & UPPER_MASK) | (state->key[num + 1] & LOWER_MASK);
+    state->key[num] = state->key[num + (M - N)] ^ (y >> 1) ^ mag02[y % 2];
+    state->pos++;
+  } else if (num == N - 1) {
+    y = (state->key[N - 1] & UPPER_MASK) | (state->key[0] & LOWER_MASK);
+    state->key[N - 1] = state->key[M - 1] ^ (y >> 1) ^ mag02[y % 2];
+    state->pos = 0;
+  }
+}
+
+void add_state(mt19937_state *state1, mt19937_state *state2) {
+  int i, pt1 = state1->pos, pt2 = state2->pos;
+
+  if (pt2 - pt1 >= 0) {
+    for (i = 0; i < N - pt2; i++)
+      state1->key[i + pt1] ^= state2->key[i + pt2];
+    for (; i < N - pt1; i++)
+      state1->key[i + pt1] ^= state2->key[i + (pt2 - N)];
+    for (; i < N; i++)
+      state1->key[i + (pt1 - N)] ^= state2->key[i + (pt2 - N)];
+  } else {
+    for (i = 0; i < N - pt1; i++)
+      state1->key[i + pt1] ^= state2->key[i + pt2];
+    for (; i < N - pt2; i++)
+      state1->key[i + (pt1 - N)] ^= state2->key[i + pt2];
+    for (; i < N; i++)
+      state1->key[i + (pt1 - N)] ^= state2->key[i + (pt2 - N)];
+  }
+}
+
+/*
+void gen_vec_h(mt19937_state *state, mt19937_state *vec_h,
+               unsigned long *h) {
+  int i;
+  unsigned long k, g;
+  mt19937_state v;
+
+  gray_code(h);
+
+  copy_state(&vec_h[0], state);
+
+  for (i = 0; i < QQ; i++)
+    gen_next(&vec_h[0]);
+
+  for (i = 1; i < LL; i++) {
+    copy_state(&v, state);
+    g = h[i] ^ h[i - 1];
+    for (k = 1; k < g; k = (k << 1))
+      gen_next(&v);
+    copy_state(&vec_h[h[i]], &vec_h[h[i - 1]]);
+    add_state(&vec_h[h[i]], &v);
+  }
+}
+*/
+
+/* compute pf(ss) using Sliding window algorithm */
+/*
+void calc_state(unsigned long *pf, mt19937_state *state,
+                mt19937_state *vec_h) {
+  mt19937_state *temp1;
+  int i = MEXP - 1, j, digit, skip = 0;
+
+  temp1 = (mt19937_state *)calloc(1, sizeof(mt19937_state));
+
+  while (get_coef(pf, i) == 0)
+    i--;
+
+  for (; i >= QQ; i--) {
+    if (get_coef(pf, i) != 0) {
+      for (j = 0; j < QQ + 1; j++)
+        gen_next(temp1);
+      digit = 0;
+      for (j = 0; j < QQ; j++)
+        digit = (digit << 1) ^ get_coef(pf, i - j - 1);
+      add_state(temp1, &vec_h[digit]);
+      i -= QQ;
+    } else
+      gen_next(temp1);
+  }
+
+  for (; i > -1; i--) {
+    gen_next(temp1);
+    if (get_coef(pf, i) == 1)
+      add_state(temp1, state);
+    else
+      ;
+  }
+
+  copy_state(state, temp1);
+  free(temp1);
+}
+*/
+
+/* compute pf(ss) using standard Horner method */
+void horner1(unsigned long *pf, mt19937_state *state) {
+  int i = MEXP - 1;
+  mt19937_state *temp;
+
+  temp = (mt19937_state *)calloc(1, sizeof(mt19937_state));
+
+  while (get_coef(pf, i) == 0)
+    i--;
+
+  if (i > 0) {
+    copy_state(temp, state);
+    gen_next(temp);
+    i--;
+    for (; i > 0; i--) {
+      if (get_coef(pf, i) != 0)
+        add_state(temp, state);
+      else
+        ;
+      gen_next(temp);
+    }
+    if (get_coef(pf, 0) != 0)
+      add_state(temp, state);
+    else
+      ;
+  } else if (i == 0)
+    copy_state(temp, state);
+  else
+    ;
+
+  copy_state(state, temp);
+  free(temp);
+}
+
+void mt19937_jump_state(mt19937_state *state, const char *jump_str) {
+  unsigned long *pf;
+  int i;
+
+  pf = (unsigned long *)calloc(P_SIZE, sizeof(unsigned long));
+
+  for (i = MEXP - 1; i > -1; i--) {
+    if (jump_str[i] == '1')
+      set_coef(pf, i, 1);
+  }
+  /* TODO: Should generate the next set and start from 0, but doesn't matter ??
+   */
+  if (state->pos >= N) {
+    state->pos = 0;
+  }
+
+  horner1(pf, state);
+
+  free(pf);
+}
+/*
+void mt19937_jump(mt19937_state *state, const char *jump_str)
+{
+    unsigned long h[LL];
+    mt19937_state vec_h[LL];
+    unsigned long *pf;
+    int i;
+
+    pf = (unsigned long *)calloc(P_SIZE, sizeof(unsigned long));
+
+    for (i = MEXP - 1; i > -1; i--)
+    {
+        if (jump_str[i] == '1')
+            set_coef(pf, i, 1);
+    }
+
+    gen_vec_h(state, &vec_h, &h);
+    calc_state(pf, state, &vec_h);
+
+    free(pf);
+}
+*/
\ No newline at end of file
diff --git a/numpy/random/src/mt19937/mt19937-jump.h b/numpy/random/src/mt19937/mt19937-jump.h
new file mode 100644
index 000000000..394c150a0
--- /dev/null
+++ b/numpy/random/src/mt19937/mt19937-jump.h
@@ -0,0 +1,15 @@
+#pragma once
+#include "mt19937.h"
+#include <stdlib.h>
+
+/* parameters for computing Jump */
+#define W_SIZE 32 /* size of unsigned long */
+#define MEXP 19937
+#define P_SIZE ((MEXP / W_SIZE) + 1)
+#define LSB 0x00000001UL
+#define QQ 7
+#define LL 128 /* LL = 2^(QQ) */
+
+void mt19937_jump_state(mt19937_state *state, const char *jump_str);
+
+void set_coef(unsigned long *pf, unsigned int deg, unsigned long v);
\ No newline at end of file
diff --git a/numpy/random/src/mt19937/mt19937-poly.h b/numpy/random/src/mt19937/mt19937-poly.h
new file mode 100644
index 000000000..b03747881
--- /dev/null
+++ b/numpy/random/src/mt19937/mt19937-poly.h
@@ -0,0 +1,207 @@
+static const char * poly =
+"0001000111110111011100100010101111000000010100100101000001110111100010101000110100101001011001010"
+"1110101101100101011100101101011001110011100011110100001000001011100101100010100000010011101110011"
+"0100001001111010000100100101001011100111101101001100000111001000011101100100010000001111110100010"
+"0000111101000101000101101111001011000011001001001011010011001000001000011100100010110101111111101"
+"0010001001100010011011101111101110111010111000010000011010110011111101100000100100101001010000001"
+"1001111000011010011101001101011000111001110010110000011000110101111010110011011000001110110010001"
+"1001101011011101000011001011111111100011001010111100000001111011111101000101000011000011111100101"
+"0100001111101010101100000110100110010010101011011100110011000101100101011110010101110000101011100"
+"0001010100010110100000111001100000011101011001101000001000101101010100010101100000100011110110011"
+"0101100110111101010111100010100110100011111011100111000001110110010000000100000110101010111001111"
+"0011110010000110101101010001110010100111111111100100101010010011101111011000010111101001110110110"
+"1011101101101100110111000100101100111001011111110101001000011111010011000111110011100100001101111"
+"1001010110110001000100001001000010000000001011011100101010010100011000110101001000010101100111101"
+"0011110101100110111100000111001011011001100101111011000101001011011111110110100010001100101001100"
+"1111110011111111110111011011100011000100110011011011011001101011100110010001111100001111100100001"
+"1000100011001010100101010100111110001100111111011111100100011110011101101000110100101110010111111"
+"1001010110000101001110010110001011011010101111111001110001100100011001000010111001011011000111100"
+"1101001011110111111010011000110100001010000000101010101001111101111110101111110101110101010010100"
+"1100100101010110011111001101100110001011000101010001000110011011111101111110001100000010110110101"
+"1111110100001011101011101110111101100001111000011100000110110100100100100101011000111000100110001"
+"0110110001001000111110101111000000100100010100100101101111100011010100111101110010000001011111111"
+"1101010000011001010101111001111110001111100010100010100001011001110001010010100001011111110110111"
+"1100100100001111000111110111000100010101010110100111100001011001101001111101001110010110110011010"
+"1000010011000110000000110110110000111010010000111001100010100101010101111100010111000000011101110"
+"1100011010110001101100110000001010001100111101101011100111110111000110010011011011001101001111100"
+"1011111001100011010110101111100110111101011100000011000010001010001101001011000001111000101000100"
+"0110001011001010110000001101100000011000011110010000101000011010011110001101111111010010101100100"
+"1111010100000011011001111111011011111001101110101010110111110110101000100001011110111010100111100"
+"0000001001111100111111111000100000100100010001011001100001111100100000001111011101100010011000111"
+"0011110110100011011001110011100011011000010000000101101101001010111000010000010101111110000000100"
+"1011010100001001000011001100011000000111100111100101010100000111000000110111011101011111100010101"
+"0011001100110000010101111001000111001001010100011000110010011011101001001100101100000000111000111"
+"0111111000010010010100000101010010000100101011111111111001100101101010011010100010111001011100011"
+"1001001011010000110000111100010110110100000100110010000010010000001000110010101000110101101100100"
+"0001100001100011110110010000100000100010011001010010110111100011011000101011001100001111110110110"
+"0001100110010100011001101000100001110011011111101001101011110011011011111110111110101110010011001"
+"1000000101100000101100100000100000001011000100100001100100101101010111101010111101010001001010110"
+"0011111011001101001110110010100100000011001001111010001001100101110000000010111101000111111101010"
+"0110101110101110001001110000111110100000101101100110010001111101111011001000101110111010110111110"
+"0011001101011010001011000010000111111111101001011100110101011000000001111000101100011101011011100"
+"1111101110000000000110001110011001101100111111010001110000111110100011000100001100110010000110111"
+"1001011011001111011100000000011011000100000011000010010111000111101000011001001100011010001111000"
+"0011110010100010001101011101010011001100000010101001001101111101000111001110110000000010111101001"
+"1110110011101110111010011100101001010101100100011111100110001111011111110010100000011100110110001"
+"1011100000101000010100011101000010111100101111101100110001010001010000101110000000110100010110011"
+"1111110100101010011010100001100110110110011111110010000100001010011110010110001000000100000111000"
+"0111001010011001000010111001100110100110110101111011110111001001000101010010010011000111110010101"
+"1100110001100101001000010001101010011001110011001110001110010100010000000000000110111001010101000"
+"0111111011011101000111011001011011000101110100010001111100101110000100001011111101111101010011001"
+"0010001100011011101100010010101011001000001001010101100110001111001110011100110111111010110010001"
+"1111111101111001001101101001001010011001110000101000110010111110010110111111000100101000101011010"
+"0000101101101100000110101000101000010001111000100000111110011111111110010010001010001111011001100"
+"0011110111000000111111000100001111101110100010101011001010110110011001010010001011100001010110101"
+"0100000010101101000011001101110010000010110011000101100100000111111100011001110011010011001110000"
+"1110011110000000001001001010100000111001010110001110011100011010010010001110010011001010111100000"
+"1110000101101001011010001001010000111000010011010100001010110000101101110110011000011100111100001"
+"1001000011010001110110111001100100001111110010110010011111000010100000001101110100000000101101000"
+"0011000000100011000111110001000011100111110110000110101111101100011110100111111000000011011110110"
+"1101011010111010010001001101000110110010000010101000000001100100100000001111011001001010110100011"
+"1011000010101111010111000001001100111110000010110010011011110011111001000101111011010011010100001"
+"0110011111100001011111101010010100110001001001001000100010101011011000011100111000110101110000001"
+"1100001111100011110010000101011000010101111010001101010101100001100101100000100100000101011001100"
+"0011001000101010101010100111000100100010101000111111101010000000101010101001000101010100100111001"
+"1001100001010001100110111101010001111010011100000001001110100010010011110100001000011111100010001"
+"0010001000100110101011001110100110101110110110100101111000110101101101001000001110011010110011001"
+"0111111101011011101001111001011100001010110111000001100010110110100011010111011000111010100011000"
+"1111010110001001010000110001000101101100010100000000100001111100000010111001000011000101010100001"
+"0001101100011100010100101110010100000010011011010100000111110110000110101011011010010001110000111"
+"0110101000110101110010011100010010100111001101110110010001101001101101010100001010001110111011011"
+"1010011001010111101001011000100111001110011000000001101000001111001100001100000011001110100110011"
+"0011000110001001010111111111110110111111000111100010010101110000101100101000001010001011010100010"
+"1010010100010011101111100111010010010001110101011110110100001000001001000111001110010001001100100"
+"1100100010001010011011110100000101101011101010110110100100010001110000111010111001111011111001011"
+"0000000000011000100100100111001000101111000000110001011110101111110111100000000100101011000111011"
+"1011010011101000001011001001110001111010000100101101010111001010001000100001000111011010000110111"
+"1010110001001110001100001110011000101100000101100000000110101000000110101100100101110001100100100"
+"0110000110101011100001010001010000011101111011111011011000100100101011110101111000001011110010110"
+"0111011011100111101010110001111011010011111000010111110100001001010001011001000110111100000101011"
+"0010111111010100000110111101001100000100001011101010100011010010000001101100100101001000100011000"
+"0101010111111100100000111011101111100000011011111111010001100011001100101101011110101011101100001"
+"0100010011101111111011000111111101001000101101111001111000101110010111001010101011010111000000101"
+"0110010000010010101111100010111110000000011101001000011111001011111100111100100101100101111010110"
+"1010101001110011111100111110100000111100100000111111000010100001111011111110110010001001000000000"
+"1110100110010111100101111111001010001111001101100001011000111011100010100001000010100000011001000"
+"0000111000110111001001100010111010100111111001111101100101000011001001110011100110101110001101110"
+"1110000010110110010110000111001110110000011011100111000101100101000000001110011011001001111001111"
+"0000101100001000000111100110110000110111111001101001111111010000001011110011011011100100110000110"
+"1001011111101100100111111000000010001110111011010011011101001100000011001010000010101111111010110"
+"0001000100101110101101100001001010100110010000110110100110011001000111011110110011001110111110101"
+"0000011111011011001111010010101011000010011101001011100001010001111001000110000010000101010011111"
+"0110011000001111101001110001101011111111001010010110100001101000000011101000101011101000110101111"
+"0000101110011010010000110100000101100011000100101111100011001111011101001010100111001110100001101"
+"0000110111011000000110011001101011110000101100110110000101100000110110100001001001110001110001001"
+"1100110111111100101001100010010110011011110001000111111111001101111110010000011001011010111101001"
+"1101111110101110110100101100110001101101001010111101101000000011111111100101000101110001000011001"
+"1000111110111011010010101011110110110001010001001001100111111010011101111000000111011000011010011"
+"0111010101001110010100101101000110000110001100010101001110101011010100000110110111111111110011110"
+"0100011110100011001000110101111010000001011011110101001100111100010100101100010000010110011001111"
+"0011011110001110010010100100011111110000110011011100010110110101001110011010101111011001010101011"
+"1001001111001000001100100111000001000110110101100111000101011000000100001000100010011000001110011"
+"0000111100000111001101011111010000010001100000010101101000111100001000010011110000001011001001100"
+"0011011011111011100000111101001011101000010010001001111110010101111010110101101110110111010000101"
+"1100011000000000110110100011010100100010001101010101101110110111111011010110011101011010110101011"
+"1101000000010010011111000000101000110001000011100001101111010101100000100000100111111111100000000"
+"0011100011100101110010111100010111110010101110101000011000111111001110111111000001101101011011111"
+"1100110101001000011111001111000000001010001001010101101000001100111010101100010111001001111100000"
+"1110101101110001011100011101101100001001001011100111100110011101111000100010010001111100001010010"
+"1011001001010100101100010010000110010000101010111111001000011100000000101101110010001101110101001"
+"1110000011100101010000011110000010001000001010110001010000100111001100110001111000100100011100110"
+"1100010011110111001001100000100111001010000000000011100011111111101110010101111010100010000100001"
+"0101101001010111111110000110110010100000001011110100010111110111010000001011110110111000000110010"
+"0001100100111110001100010101000010011111100000100010000101110000111001101100100000011111111100010"
+"1001101101001000001111000100100001010110111011110110001001010001110001001100011001001100000000101"
+"1100011110101101011001100001010110001010000111100000011011011001000010101100010101110011001101110"
+"0000101011010001010011111001011000010101010100110110111110101000111110001000010100000000100010100"
+"1000111111000110110010001111000010101011101101111101011110101111100111111100111101000101000010011"
+"0010111010100010011001000000010111100010000101001011001101100011100001001111010100100110101111111"
+"1000010011110101001010011111111011101001110100001001100010000100001001100101101111011100100011001"
+"1111010001011001111101011110101101000111110101001010011101010010010101001000000000011001100110001"
+"0001000010101010101000010100111000001110000111001110001101111111000010101010111001011101001001011"
+"0011001111011010101110101111110001001100100111010001011000010100000100000001001100000011000011101"
+"1100000110000001011001110000101001010111101000110101000011000000111011100101010000111000010010101"
+"1010100101100001011011011110110011000100100101010011111101000000100001001101000011000101010111101"
+"1110111111100010111000111000010110111010010110000000000100101001000111101101100000000110111011001"
+"0100000000100100011110111011101101101101010110001110100001100001001011000000111111110100011110011"
+"0010000010000000010100110011110000000010000011111000111101011110000000000010101101001100000010010"
+"1011001001101110110011100001100011101001101011110011010001011101000100011111001010100000011111111"
+"1010101100000010001000110000110000101000110100110011100000110010110100011111010001000011100001001"
+"1000101000010111111011100010111000111001010100110000000010011011101010101111000110001000110111011"
+"1011100001100011010001101011010100110110011100000010111001011111110010100110100010001100000011100"
+"0001011001011000101011010000001010011010001011000111000011000011110011111001111010001101011010010"
+"0010010001001001101000101001011011101110001100010001010100010111111001100100000010001111100010111"
+"0100001111001100101001011101010010110010100010001100011010100110000100011010111110001011011001000"
+"1001001111011010010011101110100001111100000110101001010111110001101100110010111010111001011111010"
+"1110111011111110000001110010000010011111000111011011000011000010011110011111111101100101111011100"
+"0101101100000110101110000111111111111010110101010100111000011111011001100000100000101011000101110"
+"1011010010100000000100100000010111101110111001000011111011111110100011010010000110001101111101100"
+"1100010111001011011001011001010100100110100101001000111011011001100011001010010101111001100100110"
+"1000110000111011100101110101101000011001010010100011000001111001110110101101010010110110001100100"
+"0100001011101100111001010001111011010110010010110010110111110001001001111001111010010001010001101"
+"1110100110101100011110100100110111000111010110011000100100110110001101111100111110100001000110000"
+"1110011011001101100101100000001010100011101000010100111011111100011010000110000001011100010000101"
+"0100101000010001110010001100010110011111111101111000011001110111011100110010010100100010001000010"
+"0100001110010000011000110001101011101001110100100011011001000111010101110100110011010111001100001"
+"0100001001101010010111110101110111000000010100111101011010101001000001001000001000101101111000000"
+"0110000101110100001111001101110111011110010111101000100101110111010101001101100001110001101101101"
+"0010101100100101000100100100110111000111000111100111000001100001000111101011000110111110001010000"
+"0100110010001101100011010111000111010111000111110000110000101111101110010110111001011000111010001"
+"1011000010010101010010011001000011010110111011010001001010100111001000010110110110101110000110000"
+"1110110010011001011011000100011101001001000111011100100000000000100001101101000101000100000111001"
+"0011100001100110101011011101110111101111000100100011100001010001011001110010101010001110101101110"
+"1011001110111111111010101101000010111111011011011100011100101010001011011100011111011100101011000"
+"1000110100101000011111010011110000000101101110010000101100001000100000000010010110000000000110011"
+"1000000000001111001001000100000111001110111111001111101100001100111000101100011000100111111110011"
+"1110010101011010111100110010110001010000101111111101001010100010001001111010111000010000010010001"
+"1111111101100100001101011011100001010101000111110111111101011010011111111101000111011001011011000"
+"0000101011100011101110110011101111011110011110010000011001111001110111011011111010011011001110111"
+"0101100111110100000100010110010010101001010100010111000101111001011011001001110010100011101111110"
+"1101011110010101101011010010011111110000011010011101000000010000111010100100111110111000001101010"
+"0101100001111001111010101011110001001010000011010110010100011100100100111110100110000010011111001"
+"0100010011001001010101110111111010011101101100000101011110111010011110001111110100111011110011010"
+"0111001010110101010110000011001010000000101101010101001101011000011011010110101010101111101101100"
+"1100100000111101010111011011011110011001100010010000010100101000111111101011100111010101011000111"
+"1100110010101100010011111100000110011111101011100100001110001100001010101001001100010011001000100"
+"1101101000101101110010000001101001001110101111000110111000011101111110100100110111000000101011110"
+"0001100100001010101001101111001000001100000011010000100101100000001110100010010000110110101010111"
+"1100010100000110011100101010111110010110111100000010110011011001011110111001010011011110010001110"
+"1101110000001011101101011111101011111110110110000111110011101100110100010000100000110100010010110"
+"0011000011000110101001110100111010110000100010110101110111100010110001000111100111001011011110010"
+"0001001110101001101101011010111001001101100011101001011011001110011010001010110100111001111100101"
+"1000111001010010000010111010101110001100110111111000011101001000001010010011101000111001100111110"
+"1110100100100110010111111101010011101111011011111011011010011110100101100001011000001001001010010"
+"1100001000000110110011011101010001011110010001001110110100100001101101001011101010001110111111010"
+"1100011100101000011110111110110011111111100010110010110111010010001111101110011011010110000001000"
+"0010110100010101110100001000010011100110001110001110010100010010010110011100100110010100001110011"
+"1100001011010000001101011011011110100000001110100111001000101000001000001001000010000111010000100"
+"0111100000101010110010111010010101100000001100110101001001000110001110111011110001010010010011000"
+"1100001111101101100001111000101100110010001000111001101101011110100110100011101000011111011010101"
+"0101000011111010010110001001100110110111000100100011011101000010001010110001111001111101110001111"
+"0100100000010111010011111110000101001001011110100100010011101110011010100101100001010000001110100"
+"0011111101111000100110011000011001100100001010110011111100111010100011110100010101011110011001000"
+"0000110000100100001011101110111010001001011110010101111100001111101101111011011110001010000100010"
+"1001100100100100110010010101100110000000100000000111110011100111101001010000010000000000101011100"
+"0011101011100110000001100101010101011111111011010011110010011011001010011101010010100010001011010"
+"1100010011101011010111110100001010100011000011001001011011101111110011001110010001100101011001101"
+"0100010001111111100000101000001011010100011100111011010111001100110110001100110101000011010001010"
+"1011100001001010011110001010100100001101110011101011100100101010001100110011110010001100100001000"
+"0110001001110110010111101011101101010111001010011010101110000010100010000111011000010110011000001"
+"1000110010100001110001100010010000001101111110000010010110100000000000001111110010001110111100001"
+"0100111101000011101110010101011011000101011010111100111111001011110001110011110011011010010111101"
+"1010111011101101000001110111001010011001110010100100100100001010001100101010111001110100000110111"
+"1010000111000011101101100101101001100000011100100111100110010110011100101000111110111000110111110"
+"1101100101011101100111011111111001111000011110111110101100000111000101100100110111000010100101000"
+"0110000011011101111101111000110101011000010111010000111011000000100011101010100111001111101010111"
+"0001110100001000100001011101001010001110100000101100001011101111100111101011111001111100101101111"
+"0101100001110011111110110100110010000011011111101101110110000110110011100110111000111101000010111"
+"0111101011100100000000011101111011000100001000111000000111011010101010110000111111101010110001111"
+"0000110100111101111011001010101110000011001101001101000010011001101011111110111101010111010011100"
+"0101010011001111101111001100101000101000111110111001011111100000001101111011000001001100111111111"
+"1010111101000001111011110010001001001110100111110010000011110000011000000101001100011110110011001"
+"1010101001000010001010110000010011110101011110010111010001010111101100001001100011101001111101001"
+"0110110100111001110011100011111010010010100010111000001100001011010010000100100110101010111001001"
+"0110000101011011011100110111111001010000001001011010101010010001011010111100111010101101000101101"
+"0100100001011101110111111001111111110110111011000101010000010000011111001000100101100100100110110"
+"1100000111110010110011010100000100011111110001110010110001000001001111001101110110110101101010111"
+"0000100111101100010001110010110111100011100101100011";
diff --git a/numpy/random/src/mt19937/mt19937-test-data-gen.c b/numpy/random/src/mt19937/mt19937-test-data-gen.c
new file mode 100644
index 000000000..4f4ec1d64
--- /dev/null
+++ b/numpy/random/src/mt19937/mt19937-test-data-gen.c
@@ -0,0 +1,59 @@
+/*
+ * Generate testing csv files
+ *
+ * cl mt19937-test-data-gen.c randomkit.c
+ *   -IC:\Anaconda\Lib\site-packages\numpy\core\include -IC:\Anaconda\include
+ *   Advapi32.lib Kernel32.lib C:\Anaconda\libs\python36.lib  -DRK_NO_WINCRYPT=1
+ *
+ */
+#include "randomkit.h"
+#include <inttypes.h>
+#include <stdio.h>
+
+#define N 1000
+
+int main() {
+  uint64_t sum = 0;
+  uint32_t seed = 0xDEADBEAF;
+  int i;
+  rk_state state;
+  rk_seed(seed, &state);
+  uint64_t store[N];
+  for (i = 0; i < N; i++) {
+    store[i] = (uint64_t)rk_random(&state);
+  }
+
+  FILE *fp;
+  fp = fopen("mt19937-testset-1.csv", "w");
+  if (fp == NULL) {
+    printf("Couldn't open file\n");
+    return -1;
+  }
+  fprintf(fp, "seed, 0x%" PRIx32 "\n", seed);
+  for (i = 0; i < N; i++) {
+    fprintf(fp, "%d, 0x%" PRIx64 "\n", i, store[i]);
+    if (i == 999) {
+      printf("%d, 0x%" PRIx64 "\n", i, store[i]);
+    }
+  }
+  fclose(fp);
+
+  seed = 0;
+  rk_seed(seed, &state);
+  for (i = 0; i < N; i++) {
+    store[i] = (uint64_t)rk_random(&state);
+  }
+  fp = fopen("mt19937-testset-2.csv", "w");
+  if (fp == NULL) {
+    printf("Couldn't open file\n");
+    return -1;
+  }
+  fprintf(fp, "seed, 0x%" PRIx32 "\n", seed);
+  for (i = 0; i < N; i++) {
+    fprintf(fp, "%d, 0x%" PRIx64 "\n", i, store[i]);
+    if (i == 999) {
+      printf("%d, 0x%" PRIx64 "\n", i, store[i]);
+    }
+  }
+  fclose(fp);
+}
diff --git a/numpy/random/src/mt19937/mt19937.c b/numpy/random/src/mt19937/mt19937.c
new file mode 100644
index 000000000..e5ca9e0cf
--- /dev/null
+++ b/numpy/random/src/mt19937/mt19937.c
@@ -0,0 +1,107 @@
+#include "mt19937.h"
+#include "mt19937-jump.h"
+#include "mt19937-poly.h"
+
+void mt19937_seed(mt19937_state *state, uint32_t seed) {
+  int pos;
+  seed &= 0xffffffffUL;
+
+  /* Knuth's PRNG as used in the Mersenne Twister reference implementation */
+  for (pos = 0; pos < RK_STATE_LEN; pos++) {
+    state->key[pos] = seed;
+    seed = (1812433253UL * (seed ^ (seed >> 30)) + pos + 1) & 0xffffffffUL;
+  }
+  state->pos = RK_STATE_LEN;
+}
+
+/* initializes mt[RK_STATE_LEN] with a seed */
+static void init_genrand(mt19937_state *state, uint32_t s) {
+  int mti;
+  uint32_t *mt = state->key;
+
+  mt[0] = s & 0xffffffffUL;
+  for (mti = 1; mti < RK_STATE_LEN; mti++) {
+    /*
+     * See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier.
+     * In the previous versions, MSBs of the seed affect
+     * only MSBs of the array mt[].
+     * 2002/01/09 modified by Makoto Matsumoto
+     */
+    mt[mti] = (1812433253UL * (mt[mti - 1] ^ (mt[mti - 1] >> 30)) + mti);
+    /* for > 32 bit machines */
+    mt[mti] &= 0xffffffffUL;
+  }
+  state->pos = mti;
+  return;
+}
+
+/*
+ * initialize by an array with array-length
+ * init_key is the array for initializing keys
+ * key_length is its length
+ */
+void mt19937_init_by_array(mt19937_state *state, uint32_t *init_key,
+                           int key_length) {
+  /* was signed in the original code. RDH 12/16/2002 */
+  int i = 1;
+  int j = 0;
+  uint32_t *mt = state->key;
+  int k;
+
+  init_genrand(state, 19650218UL);
+  k = (RK_STATE_LEN > key_length ? RK_STATE_LEN : key_length);
+  for (; k; k--) {
+    /* non linear */
+    mt[i] = (mt[i] ^ ((mt[i - 1] ^ (mt[i - 1] >> 30)) * 1664525UL)) +
+            init_key[j] + j;
+    /* for > 32 bit machines */
+    mt[i] &= 0xffffffffUL;
+    i++;
+    j++;
+    if (i >= RK_STATE_LEN) {
+      mt[0] = mt[RK_STATE_LEN - 1];
+      i = 1;
+    }
+    if (j >= key_length) {
+      j = 0;
+    }
+  }
+  for (k = RK_STATE_LEN - 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 >= RK_STATE_LEN) {
+      mt[0] = mt[RK_STATE_LEN - 1];
+      i = 1;
+    }
+  }
+
+  mt[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */
+}
+
+void mt19937_gen(mt19937_state *state) {
+  uint32_t y;
+  int i;
+
+  for (i = 0; i < N - M; i++) {
+    y = (state->key[i] & UPPER_MASK) | (state->key[i + 1] & LOWER_MASK);
+    state->key[i] = state->key[i + M] ^ (y >> 1) ^ (-(y & 1) & MATRIX_A);
+  }
+  for (; i < N - 1; i++) {
+    y = (state->key[i] & UPPER_MASK) | (state->key[i + 1] & LOWER_MASK);
+    state->key[i] = state->key[i + (M - N)] ^ (y >> 1) ^ (-(y & 1) & MATRIX_A);
+  }
+  y = (state->key[N - 1] & UPPER_MASK) | (state->key[0] & LOWER_MASK);
+  state->key[N - 1] = state->key[M - 1] ^ (y >> 1) ^ (-(y & 1) & MATRIX_A);
+
+  state->pos = 0;
+}
+
+extern inline uint64_t mt19937_next64(mt19937_state *state);
+
+extern inline uint32_t mt19937_next32(mt19937_state *state);
+
+extern inline double mt19937_next_double(mt19937_state *state);
+
+void mt19937_jump(mt19937_state *state) { mt19937_jump_state(state, poly); }
diff --git a/numpy/random/src/mt19937/mt19937.h b/numpy/random/src/mt19937/mt19937.h
new file mode 100644
index 000000000..8105329ec
--- /dev/null
+++ b/numpy/random/src/mt19937/mt19937.h
@@ -0,0 +1,69 @@
+#pragma once
+#include <math.h>
+#ifdef _WIN32
+#if _MSC_VER == 1500
+#include "../common/stdint.h"
+#else
+#include <stdint.h>
+#endif
+#else
+#include <stdint.h>
+#endif
+
+#ifdef _WIN32
+#define inline __forceinline
+#endif
+
+#define RK_STATE_LEN 624
+
+#define N 624
+#define M 397
+#define MATRIX_A 0x9908b0dfUL
+#define UPPER_MASK 0x80000000UL
+#define LOWER_MASK 0x7fffffffUL
+
+typedef struct s_mt19937_state {
+  uint32_t key[RK_STATE_LEN];
+  int pos;
+} mt19937_state;
+
+extern void mt19937_seed(mt19937_state *state, uint32_t seed);
+
+extern void mt19937_gen(mt19937_state *state);
+
+/* Slightly optimized reference implementation of the Mersenne Twister */
+static inline uint32_t mt19937_next(mt19937_state *state) {
+  uint32_t y;
+
+  if (state->pos == RK_STATE_LEN) {
+    // Move to function to help inlining
+    mt19937_gen(state);
+  }
+  y = state->key[state->pos++];
+
+  /* Tempering */
+  y ^= (y >> 11);
+  y ^= (y << 7) & 0x9d2c5680UL;
+  y ^= (y << 15) & 0xefc60000UL;
+  y ^= (y >> 18);
+
+  return y;
+}
+
+extern void mt19937_init_by_array(mt19937_state *state, uint32_t *init_key,
+                                  int key_length);
+
+static inline uint64_t mt19937_next64(mt19937_state *state) {
+  return (uint64_t)mt19937_next(state) << 32 | mt19937_next(state);
+}
+
+static inline uint32_t mt19937_next32(mt19937_state *state) {
+  return mt19937_next(state);
+}
+
+static inline double mt19937_next_double(mt19937_state *state) {
+  int32_t a = mt19937_next(state) >> 5, b = mt19937_next(state) >> 6;
+  return (a * 67108864.0 + b) / 9007199254740992.0;
+}
+
+void mt19937_jump(mt19937_state *state);
diff --git a/numpy/random/src/mt19937/randomkit.c b/numpy/random/src/mt19937/randomkit.c
new file mode 100644
index 000000000..f8ed4b49e
--- /dev/null
+++ b/numpy/random/src/mt19937/randomkit.c
@@ -0,0 +1,578 @@
+/* Random kit 1.3 */
+
+/*
+ * Copyright (c) 2003-2005, Jean-Sebastien Roy (js@jeannot.org)
+ *
+ * The rk_random and rk_seed functions algorithms and the original design of
+ * the Mersenne Twister RNG:
+ *
+ *   Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
+ *   All rights reserved.
+ *
+ *   Redistribution and use in source and binary forms, with or without
+ *   modification, are permitted provided that the following conditions
+ *   are met:
+ *
+ *   1. Redistributions of source code must retain the above copyright
+ *   notice, this list of conditions and the following disclaimer.
+ *
+ *   2. Redistributions in binary form must reproduce the above copyright
+ *   notice, this list of conditions and the following disclaimer in the
+ *   documentation and/or other materials provided with the distribution.
+ *
+ *   3. The names of its contributors may not be used to endorse or promote
+ *   products derived from this software without specific prior written
+ *   permission.
+ *
+ *   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+ *   "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ *   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+ *   A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER
+ * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+ *   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+ *   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+ *   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+ *   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+ *   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ *   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * Original algorithm for the implementation of rk_interval function from
+ * Richard J. Wagner's implementation of the Mersenne Twister RNG, optimised by
+ * Magnus Jonsson.
+ *
+ * Constants used in the rk_double implementation by Isaku Wada.
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the
+ * "Software"), to deal in the Software without restriction, including
+ * without limitation the rights to use, copy, modify, merge, publish,
+ * distribute, sublicense, and/or sell copies of the Software, and to
+ * permit persons to whom the Software is furnished to do so, subject to
+ * the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+ * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
+ * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
+ * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
+ * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ */
+
+/* static char const rcsid[] =
+  "@(#) $Jeannot: randomkit.c,v 1.28 2005/07/21 22:14:09 js Exp $"; */
+
+#ifdef _WIN32
+/*
+ * Windows
+ * XXX: we have to use this ugly defined(__GNUC__) because it is not easy to
+ * detect the compiler used in distutils itself
+ */
+#if (defined(__GNUC__) && defined(NPY_NEEDS_MINGW_TIME_WORKAROUND))
+
+/*
+ * FIXME: ideally, we should set this to the real version of MSVCRT. We need
+ * something higher than 0x601 to enable _ftime64 and co
+ */
+#define __MSVCRT_VERSION__ 0x0700
+#include <sys/timeb.h>
+#include <time.h>
+
+/*
+ * mingw msvcr lib import wrongly export _ftime, which does not exist in the
+ * actual msvc runtime for version >= 8; we make it an alias to _ftime64, which
+ * is available in those versions of the runtime
+ */
+#define _FTIME(x) _ftime64((x))
+#else
+#include <sys/timeb.h>
+#include <time.h>
+
+#define _FTIME(x) _ftime((x))
+#endif
+
+#ifndef RK_NO_WINCRYPT
+/* Windows crypto */
+#ifndef _WIN32_WINNT
+#define _WIN32_WINNT 0x0400
+#endif
+#include <wincrypt.h>
+#include <windows.h>
+
+#endif
+
+/*
+ * Do not move this include. randomkit.h must be included
+ * after windows timeb.h is included.
+ */
+#include "randomkit.h"
+
+#else
+/* Unix */
+#include "randomkit.h"
+#include <sys/time.h>
+#include <time.h>
+#include <unistd.h>
+
+#endif
+
+#include <assert.h>
+#include <errno.h>
+#include <limits.h>
+#include <math.h>
+#include <stddef.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+#ifndef RK_DEV_URANDOM
+#define RK_DEV_URANDOM "/dev/urandom"
+#endif
+
+#ifndef RK_DEV_RANDOM
+#define RK_DEV_RANDOM "/dev/random"
+#endif
+
+char *rk_strerror[RK_ERR_MAX] = {"no error", "random device unvavailable"};
+
+/* static functions */
+static unsigned long rk_hash(unsigned long key);
+
+void rk_seed(unsigned long seed, rk_state *state) {
+  int pos;
+  seed &= 0xffffffffUL;
+
+  /* Knuth's PRNG as used in the Mersenne Twister reference implementation */
+  for (pos = 0; pos < RK_STATE_LEN; pos++) {
+    state->key[pos] = seed;
+    seed = (1812433253UL * (seed ^ (seed >> 30)) + pos + 1) & 0xffffffffUL;
+  }
+  state->pos = RK_STATE_LEN;
+  state->gauss = 0;
+  state->has_gauss = 0;
+  state->has_binomial = 0;
+}
+
+/* Thomas Wang 32 bits integer hash function */
+unsigned long rk_hash(unsigned long key) {
+  key += ~(key << 15);
+  key ^= (key >> 10);
+  key += (key << 3);
+  key ^= (key >> 6);
+  key += ~(key << 11);
+  key ^= (key >> 16);
+  return key;
+}
+
+rk_error rk_randomseed(rk_state *state) {
+#ifndef _WIN32
+  struct timeval tv;
+#else
+  struct _timeb tv;
+#endif
+  int i;
+
+  if (rk_devfill(state->key, sizeof(state->key), 0) == RK_NOERR) {
+    /* ensures non-zero key */
+    state->key[0] |= 0x80000000UL;
+    state->pos = RK_STATE_LEN;
+    state->gauss = 0;
+    state->has_gauss = 0;
+    state->has_binomial = 0;
+
+    for (i = 0; i < 624; i++) {
+      state->key[i] &= 0xffffffffUL;
+    }
+    return RK_NOERR;
+  }
+
+#ifndef _WIN32
+  gettimeofday(&tv, NULL);
+  rk_seed(rk_hash(getpid()) ^ rk_hash(tv.tv_sec) ^ rk_hash(tv.tv_usec) ^
+              rk_hash(clock()),
+          state);
+#else
+  _FTIME(&tv);
+  rk_seed(rk_hash(tv.time) ^ rk_hash(tv.millitm) ^ rk_hash(clock()), state);
+#endif
+
+  return RK_ENODEV;
+}
+
+/* Magic Mersenne Twister constants */
+#define N 624
+#define M 397
+#define MATRIX_A 0x9908b0dfUL
+#define UPPER_MASK 0x80000000UL
+#define LOWER_MASK 0x7fffffffUL
+
+/*
+ * Slightly optimised reference implementation of the Mersenne Twister
+ * Note that regardless of the precision of long, only 32 bit random
+ * integers are produced
+ */
+unsigned long rk_random(rk_state *state) {
+  unsigned long y;
+
+  if (state->pos == RK_STATE_LEN) {
+    int i;
+
+    for (i = 0; i < N - M; i++) {
+      y = (state->key[i] & UPPER_MASK) | (state->key[i + 1] & LOWER_MASK);
+      state->key[i] = state->key[i + M] ^ (y >> 1) ^ (-(y & 1) & MATRIX_A);
+    }
+    for (; i < N - 1; i++) {
+      y = (state->key[i] & UPPER_MASK) | (state->key[i + 1] & LOWER_MASK);
+      state->key[i] =
+          state->key[i + (M - N)] ^ (y >> 1) ^ (-(y & 1) & MATRIX_A);
+    }
+    y = (state->key[N - 1] & UPPER_MASK) | (state->key[0] & LOWER_MASK);
+    state->key[N - 1] = state->key[M - 1] ^ (y >> 1) ^ (-(y & 1) & MATRIX_A);
+
+    state->pos = 0;
+  }
+  y = state->key[state->pos++];
+
+  /* Tempering */
+  y ^= (y >> 11);
+  y ^= (y << 7) & 0x9d2c5680UL;
+  y ^= (y << 15) & 0xefc60000UL;
+  y ^= (y >> 18);
+
+  return y;
+}
+
+/*
+ * Returns an unsigned 64 bit random integer.
+ */
+NPY_INLINE static npy_uint64 rk_uint64(rk_state *state) {
+  npy_uint64 upper = (npy_uint64)rk_random(state) << 32;
+  npy_uint64 lower = (npy_uint64)rk_random(state);
+  return upper | lower;
+}
+
+/*
+ * Returns an unsigned 32 bit random integer.
+ */
+NPY_INLINE static npy_uint32 rk_uint32(rk_state *state) {
+  return (npy_uint32)rk_random(state);
+}
+
+/*
+ * Fills an array with cnt random npy_uint64 between off and off + rng
+ * inclusive. The numbers wrap if rng is sufficiently large.
+ */
+void rk_random_uint64(npy_uint64 off, npy_uint64 rng, npy_intp cnt,
+                      npy_uint64 *out, rk_state *state) {
+  npy_uint64 val, mask = rng;
+  npy_intp i;
+
+  if (rng == 0) {
+    for (i = 0; i < cnt; i++) {
+      out[i] = off;
+    }
+    return;
+  }
+
+  /* Smallest bit mask >= max */
+  mask |= mask >> 1;
+  mask |= mask >> 2;
+  mask |= mask >> 4;
+  mask |= mask >> 8;
+  mask |= mask >> 16;
+  mask |= mask >> 32;
+
+  for (i = 0; i < cnt; i++) {
+    if (rng <= 0xffffffffUL) {
+      while ((val = (rk_uint32(state) & mask)) > rng)
+        ;
+    } else {
+      while ((val = (rk_uint64(state) & mask)) > rng)
+        ;
+    }
+    out[i] = off + val;
+  }
+}
+
+/*
+ * Fills an array with cnt random npy_uint32 between off and off + rng
+ * inclusive. The numbers wrap if rng is sufficiently large.
+ */
+void rk_random_uint32(npy_uint32 off, npy_uint32 rng, npy_intp cnt,
+                      npy_uint32 *out, rk_state *state) {
+  npy_uint32 val, mask = rng;
+  npy_intp i;
+
+  if (rng == 0) {
+    for (i = 0; i < cnt; i++) {
+      out[i] = off;
+    }
+    return;
+  }
+
+  /* Smallest bit mask >= max */
+  mask |= mask >> 1;
+  mask |= mask >> 2;
+  mask |= mask >> 4;
+  mask |= mask >> 8;
+  mask |= mask >> 16;
+
+  for (i = 0; i < cnt; i++) {
+    while ((val = (rk_uint32(state) & mask)) > rng)
+      ;
+    out[i] = off + val;
+  }
+}
+
+/*
+ * Fills an array with cnt random npy_uint16 between off and off + rng
+ * inclusive. The numbers wrap if rng is sufficiently large.
+ */
+void rk_random_uint16(npy_uint16 off, npy_uint16 rng, npy_intp cnt,
+                      npy_uint16 *out, rk_state *state) {
+  npy_uint16 val, mask = rng;
+  npy_intp i;
+  npy_uint32 buf;
+  int bcnt = 0;
+
+  if (rng == 0) {
+    for (i = 0; i < cnt; i++) {
+      out[i] = off;
+    }
+    return;
+  }
+
+  /* Smallest bit mask >= max */
+  mask |= mask >> 1;
+  mask |= mask >> 2;
+  mask |= mask >> 4;
+  mask |= mask >> 8;
+
+  for (i = 0; i < cnt; i++) {
+    do {
+      if (!bcnt) {
+        buf = rk_uint32(state);
+        bcnt = 1;
+      } else {
+        buf >>= 16;
+        bcnt--;
+      }
+      val = (npy_uint16)buf & mask;
+    } while (val > rng);
+    out[i] = off + val;
+  }
+}
+
+/*
+ * Fills an array with cnt random npy_uint8 between off and off + rng
+ * inclusive. The numbers wrap if rng is sufficiently large.
+ */
+void rk_random_uint8(npy_uint8 off, npy_uint8 rng, npy_intp cnt, npy_uint8 *out,
+                     rk_state *state) {
+  npy_uint8 val, mask = rng;
+  npy_intp i;
+  npy_uint32 buf;
+  int bcnt = 0;
+
+  if (rng == 0) {
+    for (i = 0; i < cnt; i++) {
+      out[i] = off;
+    }
+    return;
+  }
+
+  /* Smallest bit mask >= max */
+  mask |= mask >> 1;
+  mask |= mask >> 2;
+  mask |= mask >> 4;
+
+  for (i = 0; i < cnt; i++) {
+    do {
+      if (!bcnt) {
+        buf = rk_uint32(state);
+        bcnt = 3;
+      } else {
+        buf >>= 8;
+        bcnt--;
+      }
+      val = (npy_uint8)buf & mask;
+    } while (val > rng);
+    out[i] = off + val;
+  }
+}
+
+/*
+ * Fills an array with cnt random npy_bool between off and off + rng
+ * inclusive.
+ */
+void rk_random_bool(npy_bool off, npy_bool rng, npy_intp cnt, npy_bool *out,
+                    rk_state *state) {
+  npy_intp i;
+  npy_uint32 buf;
+  int bcnt = 0;
+
+  if (rng == 0) {
+    for (i = 0; i < cnt; i++) {
+      out[i] = off;
+    }
+    return;
+  }
+
+  /* If we reach here rng and mask are one and off is zero */
+  assert(rng == 1 && off == 0);
+  for (i = 0; i < cnt; i++) {
+    if (!bcnt) {
+      buf = rk_uint32(state);
+      bcnt = 31;
+    } else {
+      buf >>= 1;
+      bcnt--;
+    }
+    out[i] = (buf & 0x00000001) != 0;
+  }
+}
+
+long rk_long(rk_state *state) { return rk_ulong(state) >> 1; }
+
+unsigned long rk_ulong(rk_state *state) {
+#if ULONG_MAX <= 0xffffffffUL
+  return rk_random(state);
+#else
+  return (rk_random(state) << 32) | (rk_random(state));
+#endif
+}
+
+unsigned long rk_interval(unsigned long max, rk_state *state) {
+  unsigned long mask = max, value;
+
+  if (max == 0) {
+    return 0;
+  }
+  /* Smallest bit mask >= max */
+  mask |= mask >> 1;
+  mask |= mask >> 2;
+  mask |= mask >> 4;
+  mask |= mask >> 8;
+  mask |= mask >> 16;
+#if ULONG_MAX > 0xffffffffUL
+  mask |= mask >> 32;
+#endif
+
+  /* Search a random value in [0..mask] <= max */
+#if ULONG_MAX > 0xffffffffUL
+  if (max <= 0xffffffffUL) {
+    while ((value = (rk_random(state) & mask)) > max)
+      ;
+  } else {
+    while ((value = (rk_ulong(state) & mask)) > max)
+      ;
+  }
+#else
+  while ((value = (rk_ulong(state) & mask)) > max)
+    ;
+#endif
+  return value;
+}
+
+double rk_double(rk_state *state) {
+  /* shifts : 67108864 = 0x4000000, 9007199254740992 = 0x20000000000000 */
+  long a = rk_random(state) >> 5, b = rk_random(state) >> 6;
+  return (a * 67108864.0 + b) / 9007199254740992.0;
+}
+
+void rk_fill(void *buffer, size_t size, rk_state *state) {
+  unsigned long r;
+  unsigned char *buf = buffer;
+
+  for (; size >= 4; size -= 4) {
+    r = rk_random(state);
+    *(buf++) = r & 0xFF;
+    *(buf++) = (r >> 8) & 0xFF;
+    *(buf++) = (r >> 16) & 0xFF;
+    *(buf++) = (r >> 24) & 0xFF;
+  }
+
+  if (!size) {
+    return;
+  }
+  r = rk_random(state);
+  for (; size; r >>= 8, size--) {
+    *(buf++) = (unsigned char)(r & 0xFF);
+  }
+}
+
+rk_error rk_devfill(void *buffer, size_t size, int strong) {
+#ifndef _WIN32
+  FILE *rfile;
+  int done;
+
+  if (strong) {
+    rfile = fopen(RK_DEV_RANDOM, "rb");
+  } else {
+    rfile = fopen(RK_DEV_URANDOM, "rb");
+  }
+  if (rfile == NULL) {
+    return RK_ENODEV;
+  }
+  done = fread(buffer, size, 1, rfile);
+  fclose(rfile);
+  if (done) {
+    return RK_NOERR;
+  }
+#else
+
+#ifndef RK_NO_WINCRYPT
+  HCRYPTPROV hCryptProv;
+  BOOL done;
+
+  if (!CryptAcquireContext(&hCryptProv, NULL, NULL, PROV_RSA_FULL,
+                           CRYPT_VERIFYCONTEXT) ||
+      !hCryptProv) {
+    return RK_ENODEV;
+  }
+  done = CryptGenRandom(hCryptProv, size, (unsigned char *)buffer);
+  CryptReleaseContext(hCryptProv, 0);
+  if (done) {
+    return RK_NOERR;
+  }
+#endif
+
+#endif
+  return RK_ENODEV;
+}
+
+rk_error rk_altfill(void *buffer, size_t size, int strong, rk_state *state) {
+  rk_error err;
+
+  err = rk_devfill(buffer, size, strong);
+  if (err) {
+    rk_fill(buffer, size, state);
+  }
+  return err;
+}
+
+double rk_gauss(rk_state *state) {
+  if (state->has_gauss) {
+    const double tmp = state->gauss;
+    state->gauss = 0;
+    state->has_gauss = 0;
+    return tmp;
+  } else {
+    double f, x1, x2, r2;
+
+    do {
+      x1 = 2.0 * rk_double(state) - 1.0;
+      x2 = 2.0 * rk_double(state) - 1.0;
+      r2 = x1 * x1 + x2 * x2;
+    } while (r2 >= 1.0 || r2 == 0.0);
+
+    /* Polar method, a more efficient version of the Box-Muller approach. */
+    f = sqrt(-2.0 * log(r2) / r2);
+    /* Keep for next call */
+    state->gauss = f * x1;
+    state->has_gauss = 1;
+    return f * x2;
+  }
+}
diff --git a/numpy/random/src/mt19937/randomkit.h b/numpy/random/src/mt19937/randomkit.h
new file mode 100644
index 000000000..abb082cb2
--- /dev/null
+++ b/numpy/random/src/mt19937/randomkit.h
@@ -0,0 +1,223 @@
+/* Random kit 1.3 */
+
+/*
+ * Copyright (c) 2003-2005, Jean-Sebastien Roy (js@jeannot.org)
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the
+ * "Software"), to deal in the Software without restriction, including
+ * without limitation the rights to use, copy, modify, merge, publish,
+ * distribute, sublicense, and/or sell copies of the Software, and to
+ * permit persons to whom the Software is furnished to do so, subject to
+ * the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+ * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
+ * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
+ * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
+ * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ */
+
+/* @(#) $Jeannot: randomkit.h,v 1.24 2005/07/21 22:14:09 js Exp $ */
+
+/*
+ * Typical use:
+ *
+ * {
+ *  rk_state state;
+ *  unsigned long seed = 1, random_value;
+ *
+ *  rk_seed(seed, &state); // Initialize the RNG
+ *  ...
+ *  random_value = rk_random(&state); // Generate random values in [0..RK_MAX]
+ * }
+ *
+ * Instead of rk_seed, you can use rk_randomseed which will get a random seed
+ * from /dev/urandom (or the clock, if /dev/urandom is unavailable):
+ *
+ * {
+ *  rk_state state;
+ *  unsigned long random_value;
+ *
+ *  rk_randomseed(&state); // Initialize the RNG with a random seed
+ *  ...
+ *  random_value = rk_random(&state); // Generate random values in [0..RK_MAX]
+ * }
+ */
+
+/*
+ * Useful macro:
+ *  RK_DEV_RANDOM: the device used for random seeding.
+ *                 defaults to "/dev/urandom"
+ */
+
+#ifndef _RANDOMKIT_
+#define _RANDOMKIT_
+
+#include <numpy/npy_common.h>
+#include <stddef.h>
+
+#define RK_STATE_LEN 624
+
+typedef struct rk_state_ {
+  unsigned long key[RK_STATE_LEN];
+  int pos;
+  int has_gauss; /* !=0: gauss contains a gaussian deviate */
+  double gauss;
+
+  /* The rk_state structure has been extended to store the following
+   * information for the binomial generator. If the input values of n or p
+   * are different than nsave and psave, then the other parameters will be
+   * recomputed. RTK 2005-09-02 */
+
+  int has_binomial; /* !=0: following parameters initialized for
+                            binomial */
+  double psave;
+  long nsave;
+  double r;
+  double q;
+  double fm;
+  long m;
+  double p1;
+  double xm;
+  double xl;
+  double xr;
+  double c;
+  double laml;
+  double lamr;
+  double p2;
+  double p3;
+  double p4;
+
+} rk_state;
+
+typedef enum {
+  RK_NOERR = 0,  /* no error */
+  RK_ENODEV = 1, /* no RK_DEV_RANDOM device */
+  RK_ERR_MAX = 2
+} rk_error;
+
+/* error strings */
+extern char *rk_strerror[RK_ERR_MAX];
+
+/* Maximum generated random value */
+#define RK_MAX 0xFFFFFFFFUL
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+/*
+ * Initialize the RNG state using the given seed.
+ */
+extern void rk_seed(unsigned long seed, rk_state *state);
+
+/*
+ * Initialize the RNG state using a random seed.
+ * Uses /dev/random or, when unavailable, the clock (see randomkit.c).
+ * Returns RK_NOERR when no errors occurs.
+ * Returns RK_ENODEV when the use of RK_DEV_RANDOM failed (for example because
+ * there is no such device). In this case, the RNG was initialized using the
+ * clock.
+ */
+extern rk_error rk_randomseed(rk_state *state);
+
+/*
+ * Returns a random unsigned long between 0 and RK_MAX inclusive
+ */
+extern unsigned long rk_random(rk_state *state);
+
+/*
+ * Returns a random long between 0 and LONG_MAX inclusive
+ */
+extern long rk_long(rk_state *state);
+
+/*
+ * Returns a random unsigned long between 0 and ULONG_MAX inclusive
+ */
+extern unsigned long rk_ulong(rk_state *state);
+
+/*
+ * Returns a random unsigned long between 0 and max inclusive.
+ */
+extern unsigned long rk_interval(unsigned long max, rk_state *state);
+
+/*
+ * Fills an array with cnt random npy_uint64 between off and off + rng
+ * inclusive. The numbers wrap if rng is sufficiently large.
+ */
+extern void rk_random_uint64(npy_uint64 off, npy_uint64 rng, npy_intp cnt,
+                             npy_uint64 *out, rk_state *state);
+
+/*
+ * Fills an array with cnt random npy_uint32 between off and off + rng
+ * inclusive. The numbers wrap if rng is sufficiently large.
+ */
+extern void rk_random_uint32(npy_uint32 off, npy_uint32 rng, npy_intp cnt,
+                             npy_uint32 *out, rk_state *state);
+
+/*
+ * Fills an array with cnt random npy_uint16 between off and off + rng
+ * inclusive. The numbers wrap if rng is sufficiently large.
+ */
+extern void rk_random_uint16(npy_uint16 off, npy_uint16 rng, npy_intp cnt,
+                             npy_uint16 *out, rk_state *state);
+
+/*
+ * Fills an array with cnt random npy_uint8 between off and off + rng
+ * inclusive. The numbers wrap if rng is sufficiently large.
+ */
+extern void rk_random_uint8(npy_uint8 off, npy_uint8 rng, npy_intp cnt,
+                            npy_uint8 *out, rk_state *state);
+
+/*
+ * Fills an array with cnt random npy_bool between off and off + rng
+ * inclusive. It is assumed tha npy_bool as the same size as npy_uint8.
+ */
+extern void rk_random_bool(npy_bool off, npy_bool rng, npy_intp cnt,
+                           npy_bool *out, rk_state *state);
+
+/*
+ * Returns a random double between 0.0 and 1.0, 1.0 excluded.
+ */
+extern double rk_double(rk_state *state);
+
+/*
+ * fill the buffer with size random bytes
+ */
+extern void rk_fill(void *buffer, size_t size, rk_state *state);
+
+/*
+ * fill the buffer with randombytes from the random device
+ * Returns RK_ENODEV if the device is unavailable, or RK_NOERR if it is
+ * On Unix, if strong is defined, RK_DEV_RANDOM is used. If not, RK_DEV_URANDOM
+ * is used instead. This parameter has no effect on Windows.
+ * Warning: on most unixes RK_DEV_RANDOM will wait for enough entropy to answer
+ * which can take a very long time on quiet systems.
+ */
+extern rk_error rk_devfill(void *buffer, size_t size, int strong);
+
+/*
+ * fill the buffer using rk_devfill if the random device is available and using
+ * rk_fill if it is not
+ * parameters have the same meaning as rk_fill and rk_devfill
+ * Returns RK_ENODEV if the device is unavailable, or RK_NOERR if it is
+ */
+extern rk_error rk_altfill(void *buffer, size_t size, int strong,
+                           rk_state *state);
+
+/*
+ * return a random gaussian deviate with variance unity and zero mean.
+ */
+extern double rk_gauss(rk_state *state);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif /* _RANDOMKIT_ */