/* ranuni.c: Knuth's rng_dbl.c made into a library, no other changes */
/*****************************************************************/
/*
BQ, December, 2003
* code made into a library
* generator self-initializes if user forgets to do so
*/
/*****************************************************************/

#include <float.h>
#include "ran.h"

/*    This program by D E Knuth is in the public domain and freely copyable
 *    AS LONG AS YOU MAKE ABSOLUTELY NO CHANGES!
 *    It is explained in Seminumerical Algorithms, 3rd edition, Section 3.6
 *    (or in the errata to the 2nd edition --- see
 *        http://www-cs-faculty.stanford.edu/~knuth/taocp.html
 *    in the changes to Volume 2 on pages 171 and following).              */

/*    N.B. The MODIFICATIONS introduced in the 9th printing (2002) are
      included here; there's no backwards compatibility with the original. */

/*    This version also adopts Brendan McKay's suggestion to
      accommodate naive users who forget to call ranf_start(seed).         */

/*    If you find any bugs, please report them immediately to
 *                 taocp@cs.stanford.edu
 *    (and you will be rewarded if the bug is genuine). Thanks!            */

/************ see the book for explanations and caveats! *******************/
/************ in particular, you need two's complement arithmetic **********/

#include <float.h>
#include "ran.h"

#define KK 100                     /* the long lag */
#define LL  37                     /* the short lag */
#define mod_sum(x,y) (((x)+(y))-(int)((x)+(y)))   /* (x+y) mod 1.0 */

#define QUALITY 1009 /* recommended quality level for high-res use */
static double ranf_arr_buf[QUALITY];
static double ranf_arr_dummy=-1.0, ranf_arr_started=-1.0;
static double *ranf_arr_ptr=&ranf_arr_dummy; /* the next random fraction, or -1 */
static int    not_ready = 1;

#define TT  70   /* guaranteed separation between streams */
#define is_odd(s) ((s)&1)

#define SEED0 314159L           /* selector for different streams */

static double ran_u[KK];           /* the generator state */

/*****************************************************************/

void RanUniArray(double aa[], int n) /* ranf_array(double aa[], int n)*/
/*
aa[0:n-1] <- new random fractions
n = Array length, must be at least KK
*/
{
  register int i,j;
  const double one=1, zero=0;

  if (n<KK) {for (j=0;j<n;j++) aa[j]= -one/zero; return;} /* ERROR */

  if (not_ready) RanStart(SEED0); /* the user forgot to initialize */

  for (j=0;j<KK;j++) aa[j]=ran_u[j];
  for (;j<n;j++) aa[j]=mod_sum(aa[j-KK],aa[j-LL]);
  for (i=0;i<LL;i++,j++) ran_u[i]=mod_sum(aa[j-KK],aa[j-LL]);
  for (;i<KK;i++,j++) ran_u[i]=mod_sum(aa[j-KK],ran_u[i-LL]);
}
/*****************************************************************/
/* the following routines are adapted from exercise 3.6--15 */
/* after calling ranf_start, get new randoms by, e.g., "x=ranf_arr_next()" */
/*****************************************************************/

void RanStart(long seed)
/*
call this function before using RanUniArray
seed = selector for different streams
*/
{
  register int t,s,j;
  double u[KK+KK-1];
  double ulp=(1.0/(1L<<30))/(1L<<22);        /* 2 to the -52 */
  double ss=2.0*ulp*((seed&0x3fffffff)+2);

  not_ready=0;  /* marks that user has initialized */
  for (j=0;j<KK;j++) {
    u[j]=ss;                                /* bootstrap the buffer */
    ss+=ss; if (ss>=1.0) ss-=1.0-2*ulp;  /* cyclic shift of 51 bits */
  }
  u[1]+=ulp;                     /* make u[1] (and only u[1]) "odd" */
  for (s=seed&0x3fffffff,t=TT-1; t; ) {
    for (j=KK-1;j>0;j--)
      u[j+j]=u[j],u[j+j-1]=0.0;                         /* "square" */
    for (j=KK+KK-2;j>=KK;j--) {
      u[j-(KK-LL)]=mod_sum(u[j-(KK-LL)],u[j]);
      u[j-KK]=mod_sum(u[j-KK],u[j]);
    }
    if (is_odd(s)) {                             /* "multiply by z" */
      for (j=KK;j>0;j--) u[j]=u[j-1];
      u[0]=u[KK];                    /* shift the buffer cyclically */
      u[LL]=mod_sum(u[LL],u[KK]);
    }
    if (s) s>>=1; else t--;
  }
  for (j=0;j<LL;j++) ran_u[j+KK-LL]=u[j];
  for (;j<KK;j++) ran_u[j-LL]=u[j];
  for (j=0;j<10;j++) RanUniArray(u,KK+KK-1);  /* warm things up */
  ranf_arr_ptr = &ranf_arr_started; /* marks that user has initialized */
}

/*****************************************************************/

double RanUni(void)     /* ranf_arr_next() */
{
  if (*ranf_arr_ptr >= 0)  return  *ranf_arr_ptr++;

  if (not_ready) RanStart(SEED0); /* the user forgot to initialize */

  RanUniArray(ranf_arr_buf, QUALITY);
  ranf_arr_buf[100] = -1;
  ranf_arr_ptr = ranf_arr_buf+1;
  return ranf_arr_buf[0];
}
/*****************************************************************/