/*
 *  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 <inttypes.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);
}

#if HAVE_STRTOULL
#  define STRTOULL strtoull
#elif HAVE___STRTOULL
/* HP-UX 11 has __strtoull in <inttypes.h> */
#  define STRTOULL __strtoull
#else
/*
 * Just fall back to strtoul -- in the worst case we just lose the ability
 * to set all 64 bits of the seed.
 */
#  define STRTOULL strtoul
#endif

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;
}