/******************************************/
/* random.c 0.0.2 (1999-Dec-07-Tue)       */
/* Adam M. Costello <amc@cs.berkeley.edu> */
/******************************************/

/* Implementation of the interface defined by random.h 0.0.x. */

/* This is ANSI C code. */


#include "random.h"

#include <assert.h>
#include <float.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>


static int initialized = 0;  /* boolean: init() has been called      */
static int rand_calls;       /* number of times frand() calls rand() */
static double frand_max;     /* max return value of frand()          */
static double pi;            /* circumference/diameter of a circle   */


/* Initialization function must be called before  */
/* using any static variables or calling frand(): */

static void init(void)
{
  double max, m = RAND_MAX + 1.0, p;
  int i;

  srand((unsigned int) time(0));

  for (i = 0, max = 0.0;  i < DBL_DIG;  max = 10 * max + 9, ++i);
  /* All integral doubles 0.0 to max can be distinguished. */

  for (i = 0, p = 1.0;  p * m + 1.0 <= max;  ++i, p *= m);
  /* invariant:  p == pow(m,i) */

  rand_calls = i;
  assert(i > 0);
  frand_max = p - 1.0;
  /* All integral doubles 0.0 to frand_max + 2.0 can be distinguished. */
  pi = acos(-1.0);
  initialized = 1;
}


/* Return an integral double from 0 to frand_max: */

static double frand(void)
{
  int i;
  double r = 0.0;

  for (i = 0;  i < rand_calls;  ++i) {
    r = r * (RAND_MAX + 1.0) + rand();
  }

  return r;
}


double random_uniform(double x0, int include0, double x1, int include1)
{
  double n, d;

  if (!initialized) init();
  n = !include0 + frand();
  d = !include0 + frand_max + !include1;
  return x0 + (x1 - x0) * (n / d);
}


double random_uniform01(void)
{
  if (!initialized) init();
  return frand() / frand_max;
}


long random_discrete_uniform(long n0, long n1)
{
  long min, max;

  if (n0 <= n1) {
    min = n0;
    max = n1;
  }
  else {
    min = n1;
    max = n0;
  }

  return min + (long) random_uniform(0.0, 1, max - min + 1, 0);
}


int random_bernoulli(double p)
{
  return random_uniform(0.0, 1, 1.0, 0) < p;
}


int random_bit(void)
{
  if (!initialized) init();
  return (rand() & 0x100) != 0;
}


double random_exponential(double mean)
{
  return -mean * log(random_uniform(0.0, 0, 1.0, 1));
}


double random_rayleigh(double mean)
{
  double u;

  u = log(random_uniform(0.0, 0, 1.0, 1));
  /* Now pi must be set. */
  return sqrt(-4 * mean * mean / pi * u);
}