Remove randomdecl test dependency on GSL.

[cdecl99.git] / test / rng.c

diff --git a/test/rng.c b/test/rng.c
new file mode 100644 (file)
index 0000000..884d981
--- /dev/null
+++ b/test/rng.c
@@ -0,0 +1,143 @@
+/*
+ *  Simple random number generator for testing.
+ *  Copyright © 2022 Nick Bowler
+ *
+ *  This program is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ *  This implementation is adapted from xoshiro256+ and splitmix64
+ *  by David Blackman and Sebastiano Vigna, originally distributed
+ *  under the Creative Commons Zero public domain dedication.
+ */
+
+#include <config.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <errno.h>
+#include <limits.h>
+#include <float.h>
+#include <math.h>
+
+#include "test.h"
+
+#define B64(x) ((x) & 0xffffffffffffffff)
+
+struct test_rng {
+       unsigned long long state[4];
+};
+
+/* Rotate val left by n bits.  The behaviour is undefined if n is zero. */
+static unsigned long long rot_left64(unsigned long long val, int n)
+{
+       return B64( (val << n) | (val >> (64 - n)) );
+}
+
+static unsigned long long xoshiro256p(unsigned long long *s)
+{
+       unsigned long long ret = s[0];
+       unsigned long long t;
+
+       t = B64(s[1] << 17);
+
+       s[2] ^= s[0];
+       s[3] ^= s[1];
+       s[1] ^= s[2];
+       s[0] ^= s[3];
+       s[2] ^= t;
+       s[3] = rot_left64(s[3], 45);
+
+       return ret;
+}
+
+static unsigned long long splitmix64(unsigned long long *state)
+{
+       unsigned long long z;
+
+       z = B64(*state += 0x9e3779b97f4a7c15);
+       z = B64((z ^ (z >> 30)) * 0xbf58476d1ce4e5b9);
+       z = B64((z ^ (z >> 27)) * 0x94d049bb133111eb);
+
+       return z ^ (z >> 31);
+}
+
+struct test_rng *test_rng_alloc(const char *seed_str)
+{
+       unsigned long long seed;
+       struct test_rng *rng;
+       char *end;
+
+       errno = 0;
+       seed = strtoull(seed_str, &end, 0);
+       if (*end != 0) {
+               fprintf(stderr, "%s: invalid seed\n", seed_str);
+               return NULL;
+       } else if (errno != 0) {
+               fprintf(stderr, "%s: invalid seed: %s\n",
+                               seed_str, strerror(errno));
+               return NULL;
+       }
+
+       rng = malloc_nofail(sizeof *rng);
+       rng->state[0] = splitmix64(&seed);
+       rng->state[1] = splitmix64(&seed);
+       rng->state[2] = splitmix64(&seed);
+       rng->state[3] = splitmix64(&seed);
+
+       return rng;
+}
+
+void test_rng_free(struct test_rng *rng)
+{
+       free(rng);
+}
+
+#define BITS_PER_FP_DIGIT ( FLT_RADIX <  4 ? 1 \
+                         : FLT_RADIX <  8 ? 2 \
+                         : FLT_RADIX < 16 ? 3 \
+                         : FLT_RADIX < 32 ? 4 \
+                         : 5 )
+
+double test_rng_uniform(struct test_rng *rng)
+{
+       unsigned long long val = xoshiro256p(rng->state);
+       int prec = MIN(64, DBL_MANT_DIG*BITS_PER_FP_DIGIT);
+
+       return ldexp(val >> (64-prec), -prec);
+}
+
+/* Calculate the least power of two greater than val, minus 1. */
+static unsigned rng_mask(unsigned val)
+{
+       val |= val >> 1;
+       val |= val >> 2;
+       val |= val >> 4;
+       val |= val >> 8;
+
+       if (UINT_MAX >= 65536)
+               val |= val >> 16;
+
+       return val;
+}
+
+unsigned test_rng_uniform_int(struct test_rng *rng, unsigned max)
+{
+       unsigned mask = rng_mask(max-1);
+       unsigned long long val;
+
+       do {
+               val = ( xoshiro256p(rng->state) >> 32 ) & mask;
+       } while (val >= max);
+
+       return val;
+}