/****************************************************************************/
/*                                                                          */
/*    Copyright (c) 2000 by Donald E. Knuth                                 */
/*    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 pages 171 and following of Volume 2).               */
/*                                                                          */
/*  The functions of this file implement a random number generator.         */
/*  There are two functions: ranf_array() and ranf_start().                 */
/*  The logical-and operation '&' is used  for efficiency, so this          */
/*  random generator is not strictly portable unless the computer uses      */
/*  two's complement representation for integer. It does not limit          */
/*  application of the package because almost all modern computers are      */
/*  based on two's complement arithmetic.                                   */
/*                                                                          */
/*    The function ranf_start (long} seed)                                  */
/*  initializes the generator when given a seed number between              */
/*  0 and 2^30 = 1,073,741,821.                                             */
/*    The function ranf_array (double aa[], int n) generates n new random   */
/*  numbers and places them into a given array aa, using the recurrence     */
/*    X(j) = [X(j-100) - X(j-37)] mod 2^30,                                 */
/*  that is particularly well studied to modern computers. The value n must */
/*  be at least 100; larger values like 1000 are recommended, by default,   */
/*    n = NUM_RND = 1009.                                                   */
/*  (see also the header file rnd_gen.h)                                    */
/*                                                                          */
/*********** See the book of D.Knuth for explanations and caveats! **********/

#include "rnd_gen.h"

double rnd_num[NUM_RND]; /* array of random numbers */

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

void ranf_array(double aa[], int n) /* put n new random fractions in aa */
  /* double *aa  - destination */
  /* int n       - array length (must be at least KK) */
{
  int i, j;
  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]);
/* All information about numbers that will be generated by future calls   */
/* to ranf_array() appears in ran_u, so you can make a copy of that array */
/* in the midst of a computation if you want to restart at the same point */
/* later without going all the way back to the beginning of the sequence. */
  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]);
} /* ranf_array() */


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

  for (j=0;j<KK;j++) {
    u[j]=ss; ul[j]=0.0;                  /* bootstrap the buffer */
    ss+=ss; if (ss>=1.0) ss-=1.0-2*ulp;  /* cyclic shift of 51 bits */
  }
  for (;j<KK+KK-1;j++) u[j]=ul[j]=0.0;
  u[1]+=ulp;ul[1]=ulp;         /* make u[1] (and only u[1]) "odd" */
  s=seed&0x3fffffff;
  t=TT-1; while (t) {
    for (j=KK-1;j>0;j--) ul[j+j]=ul[j],u[j+j]=u[j];   /* "square" */
    for (j=KK+KK-2;j>KK-LL;j-=2)
        ul[KK+KK-1-j]=0.0,u[KK+KK-1-j]=u[j]-ul[j];
    for (j=KK+KK-2;j>=KK;j--) if(ul[j]) {
      ul[j-(KK-LL)]=ulp-ul[j-(KK-LL)],
        u[j-(KK-LL)]=mod_sum(u[j-(KK-LL)],u[j]);
      ul[j-KK]=ulp-ul[j-KK],u[j-KK]=mod_sum(u[j-KK],u[j]);
    }
    if (is_odd(s)) {                          /* "multiply by z" */
      for (j=KK;j>0;j--)  ul[j]=ul[j-1],u[j]=u[j-1];
      ul[0]=ul[KK],u[0]=u[KK];       /* shift the buffer cyclically */
      if (ul[KK]) ul[LL]=ulp-ul[LL],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];
} /* ranf_start() */