/* rng/knuthran2002.c
 * 
 * Copyright (C) 2007 Brian Gough
 * Copyright (C) 2001, 2007 Brian Gough, Carlo Perassi
 * 
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

/*
 * This generator is taken from
 *
 * Donald E. Knuth, The Art of Computer Programming, Volume 2, Section 3.6
 * Third Edition, Addison-Wesley, 
 * 
 * The modifications introduced in the 9th printing (2002) are
 * included here; there's no backwards compatibility with the
 * original.  [ see http://www-cs-faculty.stanford.edu/~knuth/taocp.html ] 
 * 
 */

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_rng.h>

#define BUFLEN 1009             /* length of the buffer aa[] */
#define KK 100                  /* the long lag */
#define LL 37                   /* the short lag */
#define MM (1L << 30)           /* the modulus */
#define TT 70                   /* guaranteed separation between streams */

#define is_odd(x) ((x) & 1)                       /* the units bit of x */
#define mod_diff(x, y) (((x) - (y)) & (MM - 1))   /* (x - y) mod MM */

static inline void ran_array (long int aa[], unsigned int n,
                              long int ran_x[]);
static inline unsigned long int ran_get (void *vstate);
static double ran_get_double (void *vstate);
static void ran_set (void *state, unsigned long int s);

typedef struct
{
  unsigned int i;
  long int aa[BUFLEN]; 
  long int ran_x[KK];  /* the generator state */
}
ran_state_t;

static inline void
ran_array (long int aa[], unsigned int n, long int ran_x[])
{
  unsigned int i;
  unsigned int j;

  for (j = 0; j < KK; j++)
    aa[j] = ran_x[j];

  for (; j < n; j++)
    aa[j] = mod_diff (aa[j - KK], aa[j - LL]);

  for (i = 0; i < LL; i++, j++)
    ran_x[i] = mod_diff (aa[j - KK], aa[j - LL]);

  for (; i < KK; i++, j++)
    ran_x[i] = mod_diff (aa[j - KK], ran_x[i - LL]);
}

static inline unsigned long int
ran_get (void *vstate)
{
  ran_state_t *state = (ran_state_t *) vstate;

  unsigned int i = state->i;
  unsigned long int v;

  if (i  == 0)
    {
      /* fill buffer with new random numbers */
      ran_array (state->aa, BUFLEN, state->ran_x);
    }

  v = state->aa[i];

  state->i = (i + 1) % KK;

  return v;
}

static double
ran_get_double (void *vstate)
{
  ran_state_t *state = (ran_state_t *) vstate;

  return ran_get (state) / 1073741824.0;        /* RAND_MAX + 1 */
}

static void
ran_set (void *vstate, unsigned long int s)
{
  ran_state_t *state = (ran_state_t *) vstate;

  long x[KK + KK - 1];          /* the preparation buffer */

  register int j;
  register int t;
  register long ss;

  if (s == 0 ) 
    s = 314159;                 /* default seed used by Knuth */

  ss = (s + 2)&(MM-2);

  for (j = 0; j < KK; j++)
    {
      x[j] = ss;                /* bootstrap the buffer */
      ss <<= 1;
      if (ss >= MM)             /* cyclic shift 29 bits */
        ss -= MM - 2;
    }
  x[1]++;                       /* make x[1] (and only x[1]) odd */

  ss = s & (MM - 1);
  t = TT - 1;
  while (t)
    {
      for (j = KK - 1; j > 0; j--)      /* square */
        {
          x[j + j] = x[j];
          x[j + j - 1] = 0;
        }

      for (j = KK + KK - 2; j >= KK; j--)
        {
          x[j - (KK - LL)] = mod_diff (x[j - (KK - LL)], x[j]);
          x[j - KK] = mod_diff (x[j - KK], x[j]);
        }

      if (is_odd (ss))
        {                       /* multiply by "z" */
          for (j = KK; j > 0; j--)
            {
              x[j] = x[j - 1];
            }
          x[0] = x[KK];         /* shift the buffer cyclically */
          x[LL] = mod_diff (x[LL], x[KK]);
        }

      if (ss)
        ss >>= 1;
      else
        t--;
    }

  for (j = 0; j < LL; j++)
    state->ran_x[j + KK - LL] = x[j];
  for (; j < KK; j++)
    state->ran_x[j - LL] = x[j];


  for (j = 0; j< 10; j++) 
    ran_array(x, KK+KK-1, state->ran_x);  /* warm things up */

  state->i = 0;

  return;
}

static const gsl_rng_type ran_type = {
  "knuthran2002",                /* name */
  0x3fffffffUL,                 /* RAND_MAX = (2 ^ 30) - 1 */
  0,                            /* RAND_MIN */
  sizeof (ran_state_t),
  &ran_set,
  &ran_get,
  &ran_get_double
};

const gsl_rng_type *gsl_rng_knuthran2002 = &ran_type;