/**
 * @file rng.c
 *
 * @brief Handle all the random number logic.
 *
 * Random numbers are currently generated using the mersenne twister.
 */

#include <stdint.h>
#include <math.h>
#include <unistd.h>
#include <time.h>
#include <errno.h>
#ifdef LINUX
#include <sys/time.h>
#include <fcntl.h>
#endif
#include <SDL.h>

#include "lephisto.h"
#include "rng.h"

static uint32_t MT[624];      /**< Mersenne twister state. */
static uint32_t mt_y;         /**< Internal mersenne twister variable. */
static int      mt_pos = 0;   /**< Current number being used. */

static uint32_t rng_timeEntropy(void);
/* Mersenne twister. */
static void mt_initArray(uint32_t seed);
static void mt_genArray(void);
static uint32_t mt_getInt(void);

/**
 * @fn void rng_init(void)
 *
 * @brief Initialize the random subsystem.
 */
void rng_init(void) {
  uint32_t i;
  int need_init;

  need_init = 1; /* Initialize by default. */
#ifdef LINUX
  int fd;
  fd = open("/dev/urandom", O_RDONLY); /* /dev/urandom is better then time seed. */
  if(fd != -1) {
    i = sizeof(uint32_t)*624;
    if(read(fd, &MT, i) == (ssize_t)i)
      need_init = 0;
    else
      i = rng_timeEntropy();
    close(fd);
  } else
    i = rng_timeEntropy();
#else
  i = rng_timeEntropy();
#endif
  if(need_init)
    mt_initArray(i);
  for(i = 0; i < 10; i++)
    /* Generate numbers to get away from poor initial values. */
    mt_genArray();
}

/**
 * @fn static uint32_t rng_timeEntropy(void)
 *
 * @brief Use time as a source of entropy.
 *    @return A 4 byte entopy seed.
 */
static uint32_t rng_timeEntropy(void) {
  int i;
  
#ifdef WIN32
  struct _timeb tb;
  _ftime(&tb);
  i = tb.time * 1000 + tb.millitm;
#else
  struct timeval tv;
  gettimeofday(&tv, NULL);
  i = tv.tv_sec * 1000000 + tv.tv_usec;
#endif
  return i;
}

/**
 * @fn static void mt_initArray(uint32_t seed)
 *
 * @brief Generate the initial mersenne twister based on seed.
 */
static void mt_initArray(uint32_t seed) {
  int i;
  MT[0] = seed;
  for(i = 1; i < 624; i++)
    MT[i] = 1812433253 * (MT[i-1] ^ (((MT[i-1])) + i) >> 30);
  mt_pos = 0;
}

/**
 * @fn static void mt_genArray(void)
 *
 * @brief Generate a new set of random numbers for the mersenne twister.
 */
static void mt_genArray(void) {
  int i;
  for(i = 0; i < 624; i++) {
    mt_y = (MT[i] & 0x80000000) + ((MT[i] % 624) & 0x7FFFFFFF);
    if(mt_y % 2)
      /* Odd. */
      MT[i] = (MT[(i+397) % 624] ^ (mt_y >> 1)) ^ 2567483615U;
    else
      /* Even. */
      MT[i] = MT[(i+397) % 624] ^ (mt_y >> 1);
  }
  mt_pos = 0;
}

/**
 * @fn static uint32_t mt_getInt(void)
 *    @return A random 4 byte number.
 */
static uint32_t mt_getInt(void) {
  if(mt_pos >= 624) mt_genArray();

  mt_y = MT[mt_pos++];
  mt_y ^= mt_y >> 11;
  mt_y ^= (mt_y << 7)   & 2636928640U;
  mt_y ^= (mt_y << 15)  & 4022730752U;
  mt_y ^= (mt_y >> 18);

  return mt_y;
}

/**
 * @fn unsigned int randint(void)
 * 
 * @brief Get a random integer.
 *    @return A random integer.
 */
unsigned int randint(void) {
  return mt_getInt();
}

/**
 * @fn double randfp(void)
 *
 * @brief Get a random float between 0 and 1 (inclusive).
 *    @return A random float between 0 and 1 (inclusive).
 */
static double m_div = (double)(0xFFFFFFFF); /**< Number to divide by. */
double randfp(void) {
  double m = (double)mt_getInt();
  return m / m_div;
}

/**
 * @fn double Normal(double x)
 *
 * @brief Calculates the normal distribution.
 *
 * Calculate N(x) where N is the normal distribution.
 *
 * Approximates to a power series:
 *
 *  N(x) = 1 - n(x)*(b1*t + b2*t^2 + b3*t^3 + b4*t^4 + b5*t^5) + Err
 *  where t = 1 / (1 + 0.2316419*x)
 *
 * Maximum absolute error is 7.5e^-8.
 *
 *    @param x Value to calculate the normal of.
 *    @return The value of the Normal.
 */
double Normal(double x) {
  double t;
  double series;
  const double b1 =  0.319381530;
  const double b2 = -0.356563782;
  const double b3 =  1.781477937;
  const double b4 = -1.821255978;
  const double b5 =  1.330274429;
  const double p  =  0.2316419;
  const double c  =  0.39894228;

  t = 1. / (1. + p * FABS(x));
  series = (1. - c * exp( -x * x / 2.) * t *
      (t*(t*(t*(t*b5 + b4) + b3) + b2) + b1));

  return (x > 0.) ? 1. - series : series;
}

/**
 * @fn double NormalInverse(double p)
 *
 * @brief Calculate the inverse of the normal.
 *
 * Lower tail quantile for standard disribution function.
 *
 * This function returns an approximation of the inverse cumulative
 * standard normal distribution system. I.e., given P, it returns
 * an approximation to the X satisfying P = Pr{Z <= X} Where Z is a
 * random variable from the standard normal distribution.
 *
 * The algorithm uses a minimax approximation by rational functions and
 * the result has a relative error whose absolute value is less than
 * 1.15e-9.
 */
/* Coefficients in rational approximations. */
static const double a[] = {
  -3.969683028665376e+01,
  2.209460984245205e+02,
  -2.759285104469687e+02,
  1.383577518672690e+02,
  -3.066479806614716e+01,
  2.506628277459239e+00
}; /**< Inverse normal coefficients. */

static const double b[] = {
  -5.447609879822406e+01,
  1.615858368580409e+02,
  -1.556989798598866e+02,
  6.680131188771972e+01,
  -1.328068155288572e+01
}; /**< Inverse normal coefficients. */

static const double c[] = {
  -7.784894002430293e-03,
  -3.223964580411365e-01,
  -2.400758277161838e+00,
  -2.549732539343734e+00,
  4.374664141464968e+00,
  2.938163982698783e+00
}; /**< Inverse normal coefficients. */

static const double d[] = {
  7.784695709041462e-03,
  3.224671290700398e-01,
  2.445134137142996e+00,
  3.754408661907416e+00
}; /**< Inverse normal coefficients. */


#define LOW   0.02425 /**< Low area threshold. */
#define HIGH  0.97575 /**< High area threshold. */

double NormalInverse( double p ) {
  double x, e, u, q, r;

  /* Check for errors. */
  errno = 0;
  if((p < 0) || (p > 1)) {
    errno = EDOM;
    return 0.;
  }
  else if(p == 0.) {
    errno = ERANGE;
    return -HUGE_VAL /* Minus "infinity". */;
  }
  else if(p == 1.) {
    errno = ERANGE;
    return HUGE_VAL /* "Infinity". */;
  }
  /* Use different aproximations for different parts. */
  else if(p < LOW) {
    /* Rational approximation for lower region. */
    q = sqrt(-2*log(p));
    x = (((((c[0]*q+c[1])*q+c[2])*q+c[3])*q+c[4])*q+c[5]) /
       ((((d[0]*q+d[1])*q+d[2])*q+d[3])*q+1);
  }
  else if(p > HIGH) {
    /* Rational approximation for upper region. */
    q  = sqrt(-2*log(1-p));
    x = -(((((c[0]*q+c[1])*q+c[2])*q+c[3])*q+c[4])*q+c[5]) /
       ((((d[0]*q+d[1])*q+d[2])*q+d[3])*q+1);
  }
  else {
    /* Rational approximation for central region. */
    q = p - 0.5;
    r = q*q;
    x = (((((a[0]*r+a[1])*r+a[2])*r+a[3])*r+a[4])*r+a[5])*q /
       (((((b[0]*r+b[1])*r+b[2])*r+b[3])*r+b[4])*r+1);
  }

  /* Full machine precision. */
  e = 0.5 * erfc(-x / M_SQRT2) - p;
  u = e * 2.5066282746310002 /* sqrt(2*pi) */ * exp((x*x)/2);
  x =  x - u/(1 + x*u/2);

  return x;
}