/* rng/knuthran.c

 * 

 * 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

 * Third Edition

 * Addison-Wesley

 * Section 3.6

 *

 * The comments are taken from the book

 * Our comments are signed

 */

#include <config.h>

#include <stdlib.h>

#include <gsl/gsl_rng.h>

#define BUFLEN 2009             /* [Brian]: 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 evenize(x) ((x) & (MM - 2))     /* make x even */

#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 (unsigned long int aa[], unsigned int n,

                              unsigned 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;

  unsigned long int aa[BUFLEN]; /* [Carlo]: I can't pass n to ran_array like

                                   Knuth does */

  unsigned long int ran_x[KK];  /* the generator state */

}

ran_state_t;

static inline void

ran_array (unsigned long int aa[], unsigned int n, unsigned 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;

  if (i == 0)

    {

      /* fill buffer with new random numbers */

      ran_array (state->aa, BUFLEN, state->ran_x);

    }

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

  return state->aa[i];

}

static double

ran_get_double (void *vstate)

{

  ran_state_t *state = (ran_state_t *) vstate;

  return ran_get (state) / 1073741824.0;        /* [Carlo]: 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 = evenize (s + 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;

    }

  for (; j < KK + KK - 1; j++)

    x[j] = 0;

  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];

      for (j = KK + KK - 2; j > KK - LL; j -= 2)

        x[KK + KK - 1 - j] = evenize (x[j]);

      for (j = KK + KK - 2; j >= KK; j--)

        if (is_odd (x[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 */

          if (is_odd (x[KK]))

            x[LL] = mod_diff (x[LL], x[KK]);

        }

      if (ss)

        ss >>= 1;

      else

        t--;

    }

  state->i = 0;

  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];

  return;

}

static const gsl_rng_type ran_type = {

  "knuthran",                   /* name */

  0x3fffffffUL,                 /* RAND_MAX *//* [Carlo]: (2 ^ 30) - 1 */

  0,                            /* RAND_MIN */

  sizeof (ran_state_t),

  &ran_set,

  &ran_get,

  &ran_get_double

};

const gsl_rng_type *gsl_rng_knuthran = &ran_type;