/* rng/tt.c
 * 
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 James Theiler, Brian Gough
 * 
 * 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.
 */

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

/* This is the TT800 twisted GSFR generator for 32 bit integers. It
   has been superceded by MT19937 (mt.c). The period is 2^800.

   This implementation is based on tt800.c, July 8th 1996 version by
   M. Matsumoto, email: matumoto@math.keio.ac.jp

   From: Makoto Matsumoto and Yoshiharu Kurita, "Twisted GFSR
   Generators II", ACM Transactions on Modelling and Computer
   Simulation, Vol. 4, No. 3, 1994, pages 254-266. */

static inline unsigned long int tt_get (void *vstate);
static double tt_get_double (void *vstate);
static void tt_set (void *state, unsigned long int s);

#define N 25
#define M 7

typedef struct
  {
    int n;
    unsigned long int x[N];
  }
tt_state_t;

static inline unsigned long int
tt_get (void *vstate)
{
  tt_state_t *state = (tt_state_t *) vstate;

  /* this is the magic vector, a */

  const unsigned long mag01[2] =
  {0x00000000, 0x8ebfd028UL};
  unsigned long int y;
  unsigned long int *const x = state->x;
  int n = state->n;

  if (n >= N)
    {
      int i;
      for (i = 0; i < N - M; i++)
        {
          x[i] = x[i + M] ^ (x[i] >> 1) ^ mag01[x[i] % 2];
        }
      for (; i < N; i++)
        {
          x[i] = x[i + (M - N)] ^ (x[i] >> 1) ^ mag01[x[i] % 2];
        };
      n = 0;
    }

  y = x[n];
  y ^= (y << 7) & 0x2b5b2500UL;         /* s and b, magic vectors */
  y ^= (y << 15) & 0xdb8b0000UL;        /* t and c, magic vectors */
  y &= 0xffffffffUL;    /* you may delete this line if word size = 32 */

  /* The following line was added by Makoto Matsumoto in the 1996
     version to improve lower bit's correlation.  Delete this line
     to use the code published in 1994.  */

  y ^= (y >> 16);       /* added to the 1994 version */

  state->n = n + 1;

  return y;
}

static double
tt_get_double (void * vstate)
{
  return tt_get (vstate) / 4294967296.0 ;
}

static void
tt_set (void *vstate, unsigned long int s)
{
  tt_state_t *state = (tt_state_t *) vstate;

  const tt_state_t init_state =
  {0,
   {0x95f24dabUL, 0x0b685215UL, 0xe76ccae7UL,
    0xaf3ec239UL, 0x715fad23UL, 0x24a590adUL,
    0x69e4b5efUL, 0xbf456141UL, 0x96bc1b7bUL,
    0xa7bdf825UL, 0xc1de75b7UL, 0x8858a9c9UL,
    0x2da87693UL, 0xb657f9ddUL, 0xffdc8a9fUL,
    0x8121da71UL, 0x8b823ecbUL, 0x885d05f5UL,
    0x4e20cd47UL, 0x5a9ad5d9UL, 0x512c0c03UL,
    0xea857ccdUL, 0x4cc1d30fUL, 0x8891a8a1UL,
    0xa6b7aadbUL}};


  if (s == 0)   /* default seed is given explicitly in the original code */
    {
      *state = init_state;
    }
  else
    {
      int i;

      state->n = 0;

      state->x[0] = s & 0xffffffffUL;

      for (i = 1; i < N; i++)
        state->x[i] = (69069 * state->x[i - 1]) & 0xffffffffUL;
    }

  return;
}

static const gsl_rng_type tt_type =
{"tt800",                       /* name */
 0xffffffffUL,                  /* RAND_MAX */
 0,                             /* RAND_MIN */
 sizeof (tt_state_t),
 &tt_set,
 &tt_get,
 &tt_get_double};

const gsl_rng_type *gsl_rng_tt800 = &tt_type;