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

5 files changed, 433 insertions, 0 deletions

diff --git a/numpy/random/src/philox/LICENSE.md b/numpy/random/src/philox/LICENSE.md
new file mode 100644
index 000000000..4a9f6bb29
--- /dev/null
+++ b/numpy/random/src/philox/LICENSE.md
@@ -0,0 +1,31 @@
+# THREEFRY
+
+Copyright 2010-2012, D. E. Shaw Research.
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+* Redistributions of source code must retain the above copyright
+  notice, this list of conditions, and the following disclaimer.
+
+* 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.
+
+* Neither the name of D. E. Shaw Research nor the names of its
+  contributors may 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.
diff --git a/numpy/random/src/philox/philox-benchmark.c b/numpy/random/src/philox/philox-benchmark.c
new file mode 100644
index 000000000..0cab04cf5
--- /dev/null
+++ b/numpy/random/src/philox/philox-benchmark.c
@@ -0,0 +1,38 @@
+/*
+ * Simple benchamrk command
+ *
+ *  cl philox-benchmark.c /Ox
+ *
+ *  gcc philox-benchmark.c -O3 -o philox-benchmark
+ *
+ * Requres the Random123 directory containing header files to be located in the
+ * same directory (not included).
+ */
+#include "Random123/philox.h"
+#include <inttypes.h>
+#include <stdio.h>
+#include <time.h>
+
+#define N 1000000000
+
+int main() {
+  philox4x64_ctr_t ctr = {{0, 0, 0, 0}};
+  philox4x64_key_t key = {{0, 0xDEADBEAF}};
+  philox4x64_ctr_t out;
+  uint64_t count = 0, sum = 0;
+  int i, j;
+  clock_t begin = clock();
+  for (i = 0; i < N / 4UL; i++) {
+    ctr.v[0]++;
+    out = philox4x64_R(philox4x64_rounds, ctr, key);
+    for (j = 0; j < 4; j++) {
+      sum += out.v[j];
+      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)(N / time_spent) / 1000000 * 1000000);
+}
diff --git a/numpy/random/src/philox/philox-test-data-gen.c b/numpy/random/src/philox/philox-test-data-gen.c
new file mode 100644
index 000000000..442e18b55
--- /dev/null
+++ b/numpy/random/src/philox/philox-test-data-gen.c
@@ -0,0 +1,82 @@
+/*
+ * Generate testing csv files
+ *
+ *  cl philox-test-data-gen.c /Ox
+ *  philox-test-data-gen.exe
+ *
+ *  gcc philox-test-data-gen.c -o philox-test-data-gen
+ *  ./philox-test-data-gen
+ *
+ * Requres the Random123 directory containing header files to be located in the
+ * same directory (not included).
+ *
+ */
+
+#include "../splitmix64/splitmix64.h"
+#include "Random123/philox.h"
+#include <inttypes.h>
+#include <stdio.h>
+
+#define N 1000
+
+int main() {
+  philox4x64_ctr_t ctr = {{0, 0, 0, 0}};
+  philox4x64_key_t key = {{0, 0}};
+  uint64_t state, seed = 0xDEADBEAF;
+  philox4x64_ctr_t out;
+  uint64_t store[N];
+  state = seed;
+  int i, j;
+  for (i = 0; i < 2; i++) {
+    key.v[i] = splitmix64_next(&state);
+  }
+  for (i = 0; i < N / 4UL; i++) {
+    ctr.v[0]++;
+    out = philox4x64_R(philox4x64_rounds, ctr, key);
+    for (j = 0; j < 4; j++) {
+      store[i * 4 + j] = out.v[j];
+    }
+  }
+
+  FILE *fp;
+  fp = fopen("philox-testset-1.csv", "w");
+  if (fp == NULL) {
+    printf("Couldn't open file\n");
+    return -1;
+  }
+  fprintf(fp, "seed, 0x%" PRIx64 "\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);
+
+  ctr.v[0] = 0;
+  state = seed = 0;
+  for (i = 0; i < 2; i++) {
+    key.v[i] = splitmix64_next(&state);
+  }
+  for (i = 0; i < N / 4UL; i++) {
+    ctr.v[0]++;
+    out = philox4x64_R(philox4x64_rounds, ctr, key);
+    for (j = 0; j < 4; j++) {
+      store[i * 4 + j] = out.v[j];
+    }
+  }
+
+  fp = fopen("philox-testset-2.csv", "w");
+  if (fp == NULL) {
+    printf("Couldn't open file\n");
+    return -1;
+  }
+  fprintf(fp, "seed, 0x%" PRIx64 "\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/philox/philox.c b/numpy/random/src/philox/philox.c
new file mode 100644
index 000000000..3382c60d6
--- /dev/null
+++ b/numpy/random/src/philox/philox.c
@@ -0,0 +1,29 @@
+#include "philox.h"
+
+extern INLINE uint64_t philox_next64(philox_state *state);
+
+extern INLINE uint32_t philox_next32(philox_state *state);
+
+extern void philox_jump(philox_state *state) {
+  /* Advances state as-if 2^128 draws were made */
+  state->ctr->v[2]++;
+  if (state->ctr->v[2] == 0) {
+    state->ctr->v[3]++;
+  }
+}
+
+extern void philox_advance(uint64_t *step, philox_state *state) {
+  int i, carry = 0;
+  uint64_t v_orig;
+  for (i = 0; i < 4; i++) {
+    if (carry == 1) {
+      state->ctr->v[i]++;
+      carry = state->ctr->v[i] == 0 ? 1 : 0;
+    }
+    v_orig = state->ctr->v[i];
+    state->ctr->v[i] += step[i];
+    if (state->ctr->v[i] < v_orig && carry == 0) {
+      carry = 1;
+    }
+  }
+}
diff --git a/numpy/random/src/philox/philox.h b/numpy/random/src/philox/philox.h
new file mode 100644
index 000000000..411404b55
--- /dev/null
+++ b/numpy/random/src/philox/philox.h
@@ -0,0 +1,253 @@
+#ifndef _RANDOMDGEN__PHILOX_H_
+#define _RANDOMDGEN__PHILOX_H_
+
+#include <inttypes.h>
+
+#ifdef _WIN32
+#define INLINE __inline __forceinline
+#else
+#define INLINE inline
+#endif
+
+#define PHILOX_BUFFER_SIZE 4L
+
+struct r123array2x64 {
+  uint64_t v[2];
+};
+struct r123array4x64 {
+  uint64_t v[4];
+};
+
+enum r123_enum_philox4x64 { philox4x64_rounds = 10 };
+typedef struct r123array4x64 philox4x64_ctr_t;
+typedef struct r123array2x64 philox4x64_key_t;
+typedef struct r123array2x64 philox4x64_ukey_t;
+
+static INLINE struct r123array2x64
+_philox4x64bumpkey(struct r123array2x64 key) {
+  key.v[0] += (0x9E3779B97F4A7C15ULL);
+  key.v[1] += (0xBB67AE8584CAA73BULL);
+  return key;
+}
+
+#ifdef _WIN32
+#include <intrin.h>
+/* TODO: This isn't correct for many platforms */
+#ifdef _WIN64
+#pragma intrinsic(_umul128)
+#else
+#pragma intrinsic(__emulu)
+static INLINE uint64_t _umul128(uint64_t a, uint64_t b, uint64_t *high) {
+
+  uint64_t a_lo, a_hi, b_lo, b_hi, a_x_b_hi, a_x_b_mid, a_x_b_lo, b_x_a_mid,
+      carry_bit;
+  a_lo = (uint32_t)a;
+  a_hi = a >> 32;
+  b_lo = (uint32_t)b;
+  b_hi = b >> 32;
+
+  a_x_b_hi = __emulu(a_hi, b_hi);
+  a_x_b_mid = __emulu(a_hi, b_lo);
+  b_x_a_mid = __emulu(b_hi, a_lo);
+  a_x_b_lo = __emulu(a_lo, b_lo);
+
+  carry_bit = ((uint64_t)(uint32_t)a_x_b_mid + (uint64_t)(uint32_t)b_x_a_mid +
+               (a_x_b_lo >> 32)) >>
+              32;
+
+  *high = a_x_b_hi + (a_x_b_mid >> 32) + (b_x_a_mid >> 32) + carry_bit;
+
+  return a_x_b_lo + ((a_x_b_mid + b_x_a_mid) << 32);
+}
+#endif
+static INLINE uint64_t mulhilo64(uint64_t a, uint64_t b, uint64_t *hip) {
+  return _umul128(a, b, hip);
+}
+#else
+#if __SIZEOF_INT128__
+static INLINE uint64_t mulhilo64(uint64_t a, uint64_t b, uint64_t *hip) {
+  __uint128_t product = ((__uint128_t)a) * ((__uint128_t)b);
+  *hip = product >> 64;
+  return (uint64_t)product;
+}
+#else
+static INLINE uint64_t _umul128(uint64_t a, uint64_t b, uint64_t *high) {
+
+  uint64_t a_lo, a_hi, b_lo, b_hi, a_x_b_hi, a_x_b_mid, a_x_b_lo, b_x_a_mid,
+      carry_bit;
+  a_lo = (uint32_t)a;
+  a_hi = a >> 32;
+  b_lo = (uint32_t)b;
+  b_hi = b >> 32;
+
+  a_x_b_hi = a_hi * b_hi;
+  a_x_b_mid = a_hi * b_lo;
+  b_x_a_mid = b_hi * a_lo;
+  a_x_b_lo = a_lo * b_lo;
+
+  carry_bit = ((uint64_t)(uint32_t)a_x_b_mid + (uint64_t)(uint32_t)b_x_a_mid +
+               (a_x_b_lo >> 32)) >>
+              32;
+
+  *high = a_x_b_hi + (a_x_b_mid >> 32) + (b_x_a_mid >> 32) + carry_bit;
+
+  return a_x_b_lo + ((a_x_b_mid + b_x_a_mid) << 32);
+}
+static INLINE uint64_t mulhilo64(uint64_t a, uint64_t b, uint64_t *hip) {
+  return _umul128(a, b, hip);
+}
+#endif
+#endif
+
+static INLINE struct r123array4x64 _philox4x64round(struct r123array4x64 ctr,
+                                                    struct r123array2x64 key);
+
+static INLINE struct r123array4x64 _philox4x64round(struct r123array4x64 ctr,
+                                                    struct r123array2x64 key) {
+  uint64_t hi0;
+  uint64_t hi1;
+  uint64_t lo0 = mulhilo64((0xD2E7470EE14C6C93ULL), ctr.v[0], &hi0);
+  uint64_t lo1 = mulhilo64((0xCA5A826395121157ULL), ctr.v[2], &hi1);
+  struct r123array4x64 out = {
+      {hi1 ^ ctr.v[1] ^ key.v[0], lo1, hi0 ^ ctr.v[3] ^ key.v[1], lo0}};
+  return out;
+}
+
+static INLINE philox4x64_key_t philox4x64keyinit(philox4x64_ukey_t uk) {
+  return uk;
+}
+static INLINE philox4x64_ctr_t philox4x64_R(unsigned int R,
+                                            philox4x64_ctr_t ctr,
+                                            philox4x64_key_t key);
+
+static INLINE philox4x64_ctr_t philox4x64_R(unsigned int R,
+                                            philox4x64_ctr_t ctr,
+                                            philox4x64_key_t key) {
+  if (R > 0) {
+    ctr = _philox4x64round(ctr, key);
+  }
+  if (R > 1) {
+    key = _philox4x64bumpkey(key);
+    ctr = _philox4x64round(ctr, key);
+  }
+  if (R > 2) {
+    key = _philox4x64bumpkey(key);
+    ctr = _philox4x64round(ctr, key);
+  }
+  if (R > 3) {
+    key = _philox4x64bumpkey(key);
+    ctr = _philox4x64round(ctr, key);
+  }
+  if (R > 4) {
+    key = _philox4x64bumpkey(key);
+    ctr = _philox4x64round(ctr, key);
+  }
+  if (R > 5) {
+    key = _philox4x64bumpkey(key);
+    ctr = _philox4x64round(ctr, key);
+  }
+  if (R > 6) {
+    key = _philox4x64bumpkey(key);
+    ctr = _philox4x64round(ctr, key);
+  }
+  if (R > 7) {
+    key = _philox4x64bumpkey(key);
+    ctr = _philox4x64round(ctr, key);
+  }
+  if (R > 8) {
+    key = _philox4x64bumpkey(key);
+    ctr = _philox4x64round(ctr, key);
+  }
+  if (R > 9) {
+    key = _philox4x64bumpkey(key);
+    ctr = _philox4x64round(ctr, key);
+  }
+  if (R > 10) {
+    key = _philox4x64bumpkey(key);
+    ctr = _philox4x64round(ctr, key);
+  }
+  if (R > 11) {
+    key = _philox4x64bumpkey(key);
+    ctr = _philox4x64round(ctr, key);
+  }
+  if (R > 12) {
+    key = _philox4x64bumpkey(key);
+    ctr = _philox4x64round(ctr, key);
+  }
+  if (R > 13) {
+    key = _philox4x64bumpkey(key);
+    ctr = _philox4x64round(ctr, key);
+  }
+  if (R > 14) {
+    key = _philox4x64bumpkey(key);
+    ctr = _philox4x64round(ctr, key);
+  }
+  if (R > 15) {
+    key = _philox4x64bumpkey(key);
+    ctr = _philox4x64round(ctr, key);
+  }
+  return ctr;
+}
+
+typedef struct s_philox_state {
+  philox4x64_ctr_t *ctr;
+  philox4x64_key_t *key;
+  int buffer_pos;
+  uint64_t buffer[PHILOX_BUFFER_SIZE];
+  int has_uint32;
+  uint32_t uinteger;
+} philox_state;
+
+static INLINE uint64_t philox_next(philox_state *state) {
+  uint64_t out;
+  int i;
+  philox4x64_ctr_t ct;
+
+  if (state->buffer_pos < PHILOX_BUFFER_SIZE) {
+    out = state->buffer[state->buffer_pos];
+    state->buffer_pos++;
+    return out;
+  }
+  /* generate 4 new uint64_t */
+  state->ctr->v[0]++;
+  /* Handle carry */
+  if (state->ctr->v[0] == 0) {
+    state->ctr->v[1]++;
+    if (state->ctr->v[1] == 0) {
+      state->ctr->v[2]++;
+      if (state->ctr->v[2] == 0) {
+        state->ctr->v[3]++;
+      }
+    }
+  }
+  ct = philox4x64_R(philox4x64_rounds, *state->ctr, *state->key);
+  for (i = 0; i < 4; i++) {
+    state->buffer[i] = ct.v[i];
+  }
+  state->buffer_pos = 1;
+  return state->buffer[0];
+}
+
+static INLINE uint64_t philox_next64(philox_state *state) {
+  return philox_next(state);
+}
+
+static INLINE uint32_t philox_next32(philox_state *state) {
+  uint64_t next;
+
+  if (state->has_uint32) {
+    state->has_uint32 = 0;
+    return state->uinteger;
+  }
+  next = philox_next(state);
+
+  state->has_uint32 = 1;
+  state->uinteger = (uint32_t)(next >> 32);
+  return (uint32_t)(next & 0xffffffff);
+}
+
+extern void philox_jump(philox_state *state);
+
+extern void philox_advance(uint64_t *step, philox_state *state);
+
+#endif