

#define IM1 2147483563
#define IM2 2147483399
#define IA1 40014
#define IA2 40692
#define IQ1 53668
#define IQ2 52774
#define IR1 12211
#define IR2 3791
#define PRECISION 1.0e-300
#define WUP 12

#include "random.h"
#include <cmath>
#include <ctime>
#include <limits>

using namespace std;
//Random number generator of Numerical recipes in C, "ran2," period of about 2e18

randomnumber::randomnumber() {
	IM1INV = 1.0e0 / IM1;
	IMM1 = IM1-1;
	DIVISOR = (1+IMM1/TABSIZE);
	//RNMX=1.0e0 - PRECISION;

	// TODO: (RMW) Perhaps we should seed with time() or simply 1, instead of an arbitrary number like 1234
	//			 -- This might require updating regression tests if an existing program does not perform a seed.
	seed(1234);  //Seed with a basic number in case the user does not seed

}

void randomnumber::seed(long seeddouble) {
	int i;
	long j;
	
	idum = seeddouble;

	// TODO: (RMW) Perhaps seeding with 0 should in fact seed with time()
	//			and seeding with a negative number should seed with seeddouble+IM1 (compare to roll() -- "if (idum<0) idum = idum + IM1;")
	//			 -- This might require updating regression tests if an existing program does not perform a seed.

	if (idum<1) idum = 1; //prevent seeding with zero
	idum2 = idum;
	for (i=TABSIZE+WUP;i>=0;i--) {
		j = idum/IQ1;
		idum = IA1*(idum-j*IQ1)-j*IR1;
		if (idum < 0) idum = idum+IM1;
		if (i<TABSIZE) iv[i]=idum;

	}
	iy = iv[0];
}

double randomnumber::roll() {
	long j;
	int i;

	j = idum/IQ1;
	idum=IA1*(idum-j*IQ1)-j*IR1;
	if (idum<0) idum = idum + IM1;
	j = idum2/IQ2;
	idum2=IA2*(idum2-j*IQ2)-j*IR2;
	if (idum2<0) idum2 = idum2 + IM2;
	i = (int) (iy/DIVISOR);
	iy=iv[i]-idum2;
	iv[i] = idum;
	if (iy<1) iy=iy+IMM1;
	return IM1INV*iy;
}

int randomnumber::roll_int(int min, int max){
	if (min>max) return min;
	double range = max - min+1;
	return floor(roll()*range)+min;
}


// A Pseudo-Random Number Generator that is extremely fast on modern computer architectures
// and has very good random-number qualities (distribution, period etc) 
// This is based on xorshift* generator as suggested by Marsaglia.
// It is a 64-bit generator with 64 bits of state has a maximal period of 2^64−1 and passes all the "Diehard tests" and fails only the MatrixRank test of BigCrush.
// Author: Richard M. Watson 2017
UINT64 rand64::calc_next(UINT64 &x) {
	// generate the next number in their sequence by repeatedly taking 
	// the exclusive or of a number with a bit-shifted version of itself
	x ^= x >> 12; // a
	x ^= x << 25; // b
	x ^= x >> 27; // c
	// apply an invertible multiplication as a non-linear transformation
	return x * 0x2545F4914F6CDD1D; 
}

rand64::rand64(const UINT64 initSeed) {
	seed(initSeed);
}
void rand64::seed(const UINT64 seed) {
	if (seed==0)
		state=time(0);  // TODO: add more entropy by using a 64-bit timestamp or use milliseconds
	else
		state=seed;
}

// return the next random 64-bit unsigned integer
UINT64 rand64::next() {
	return calc_next(state);
}
//return the current random 64-bit unsigned integer (i.e. the last one generated by a call to next*)
UINT64 rand64::current() const { return state; }

//	generate a random number in the range [0, 1.0)
double rand64::nextDouble() {
	return next() / (double)numeric_limits<UINT64>::max();
}

//	generate a random int in the range [-INT_MAX, INT_MAX)
int rand64::nextInt() { return nextInt(numeric_limits<int>::min(), numeric_limits<int>::max()); }
//	generate a random int in the range [0, max)
int rand64::nextInt(const int max) { return nextInt(0, max); }
//	generate a random int in the range [min, max)  i.e.  min <= value < max
int rand64::nextInt(const int min, int max) {
	UINT64 diff = (UINT64)max - min;
	return (int)(floor(nextDouble()*diff) + min);
}

// function operator returns the same as next()
UINT64 rand64::operator()() { return next(); }
// function operator returns the same as nextInt(int max)
int rand64::operator()(const int max) { return nextInt(max); }
