/* This file is part of the alrand library.
   alrand is Copyright by Vincent Penquerc'h <lyrian -at- kezako -dot- net>.

   The alrand library is free software; you can redistribute it and/or
   modify it under the terms of the GNU Lesser General Public
   License as published by the Free Software Foundation; either
   version 2.1 of the License, or (at your option) any later version.

   The alrand library 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
   Lesser General Public License for more details.

   You should have received a copy of the GNU Lesser General Public
   License along with the alrand Library; if not, write to the Free
   Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
   02111-1307 USA.

   This file is an adaptation of the original integer version of the
   Mersenne Twister algorithm. Its original header follows. */


/* A C-program for MT19937: Integer version (1999/10/28)          */
/*  genrand() generates one pseudorandom unsigned integer (32bit) */
/* which is uniformly distributed among 0 to 2^32-1  for each     */
/* call. sgenrand(seed) sets initial values to the working area   */
/* of 624 words. Before genrand(), sgenrand(seed) must be         */
/* called once. (seed is any 32-bit integer.)                     */
/*   Coded by Takuji Nishimura, considering the suggestions by    */
/* Topher Cooper and Marc Rieffel in July-Aug. 1997.              */

/* This library is free software; you can redistribute it and/or   */
/* modify it under the terms of the GNU Library General Public     */
/* License as published by the Free Software Foundation; either    */
/* version 2 of the License, or (at your option) any later         */
/* version.                                                        */
/* This library 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 Library General Public License for more details.    */
/* You should have received a copy of the GNU Library General      */
/* Public License along with this library; if not, write to the    */
/* Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA   */ 
/* 02111-1307  USA                                                 */

/* Copyright (C) 1997, 1999 Makoto Matsumoto and Takuji Nishimura. */
/* Any feedback is very welcome. For any question, comments,       */
/* see http://www.math.keio.ac.jp/matumoto/emt.html or email       */
/* matumoto@math.keio.ac.jp                                        */

/* REFERENCE                                                       */
/* M. Matsumoto and T. Nishimura,                                  */
/* "Mersenne Twister: A 623-Dimensionally Equidistributed Uniform  */
/* Pseudo-Random Number Generator",                                */
/* ACM Transactions on Modeling and Computer Simulation,           */
/* Vol. 8, No. 1, January 1998, pp 3--30.                          */

#include <malloc.h>
#include <assert.h>
#include "mersenne.h"

/* Period parameters */  
#define N 624
#define M 397
#define MATRIX_A 0x9908b0df   /* constant vector a */
#define UPPER_MASK 0x80000000 /* most significant w-r bits */
#define LOWER_MASK 0x7fffffff /* least significant r bits */

/* Tempering parameters */   
#define TEMPERING_MASK_B 0x9d2c5680
#define TEMPERING_MASK_C 0xefc60000
#define TEMPERING_SHIFT_U(y)  (y >> 11)
#define TEMPERING_SHIFT_S(y)  (y << 7)
#define TEMPERING_SHIFT_T(y)  (y << 15)
#define TEMPERING_SHIFT_L(y)  (y >> 18)

typedef struct {
  alrand_uint32_t mt[N];
  int mti;
} alrandint_mersenne_twister_state;

void alrand_mersenne_twister_destroy(alrand_state *state)
{
  if (state) {
    assert(state->type==alrand_type_mersenne_twister);
    free(state);
  }
}

/* Create a state for the Mersenne Twister algorithm */
void alrand_mersenne_twister_seed(alrand_state *state,alrand_uint32_t seed)
{
  size_t i;
  alrandint_mersenne_twister_state *priv;
  assert(state);
  assert(state->type==alrand_type_mersenne_twister);
  priv=(alrandint_mersenne_twister_state*)&state->priv;

  priv->mti=N;
  for (i=0;i<N;i++) {
    priv->mt[i]=seed&0xffff0000;
    seed=69069*seed+1;
    priv->mt[i]|=(seed&0xffff0000)>>16;
    seed=69069*seed+1;
  }
}

/* Initialization by "sgenrand()" is an example. Theoretically,      */
/* there are 2^19937-1 possible states as an intial state.           */
/* This function allows to choose any of 2^19937-1 ones.             */
/* Essential bits in "seed_array[]" is following 19937 bits:         */
/*  (seed_array[0]&UPPER_MASK), seed_array[1], ..., seed_array[N-1]. */
/* (seed_array[0]&LOWER_MASK) is discarded.                          */ 
/* Theoretically,                                                    */
/*  (seed_array[0]&UPPER_MASK), seed_array[1], ..., seed_array[N-1]  */
/* can take any values except all zeros.                             */
void alrand_mersenne_twister_seed_v(alrand_state *state,alrand_uint32_t *seed)
    /* the length of seed_array[] must be at least N */
{
  size_t i;
  alrandint_mersenne_twister_state *priv;
  assert(state);
  assert(state->type==alrand_type_mersenne_twister);
  priv=(alrandint_mersenne_twister_state*)&state->priv;

  assert(seed);
  for (i=0;i<N;i++) {
    priv->mt[i]=seed[i];
  }
  priv->mti=N;
}

alrand_uint32_t alrand_mersenne_twister_generate(alrand_state *state)
{
  alrand_uint32_t y;
  alrandint_mersenne_twister_state *priv;
  static alrand_uint32_t mag01[2]={0x0, MATRIX_A};
  /* mag01[x] = x * MATRIX_A  for x=0,1 */

  assert(state);
  assert(state->type==alrand_type_mersenne_twister);
  priv=(alrandint_mersenne_twister_state*)&state->priv;

  if (priv->mti >= N) { /* generate N words at one time */
      int kk;

      for (kk=0;kk<N-M;kk++) {
          y = (priv->mt[kk]&UPPER_MASK)|(priv->mt[kk+1]&LOWER_MASK);
          priv->mt[kk] = priv->mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1];
      }
      for (;kk<N-1;kk++) {
          y = (priv->mt[kk]&UPPER_MASK)|(priv->mt[kk+1]&LOWER_MASK);
          priv->mt[kk] = priv->mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1];
      }
      y = (priv->mt[N-1]&UPPER_MASK)|(priv->mt[0]&LOWER_MASK);
      priv->mt[N-1] = priv->mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1];

      priv->mti = 0;
  }

  y = priv->mt[priv->mti++];
  y ^= TEMPERING_SHIFT_U(y);
  y ^= TEMPERING_SHIFT_S(y) & TEMPERING_MASK_B;
  y ^= TEMPERING_SHIFT_T(y) & TEMPERING_MASK_C;
  y ^= TEMPERING_SHIFT_L(y);

  return y; 
}

alrand_state *alrand_mersenne_twister_load(
  alrand_uint32_t (*read)(void*),void *ptr
)
{
  int n;
  alrand_state *state;
  alrandint_mersenne_twister_state *priv;

  assert(read);

  state=alrand_mersenne_twister_create(0);
  priv=(alrandint_mersenne_twister_state*)&state->priv;
  for (n=0;n<N;++n) priv->mt[n]=read(ptr);
  priv->mti=read(ptr);

  return state;
}

int alrand_mersenne_twister_save(
  alrand_state *state,void (*write)(alrand_uint32_t,void*),void *ptr
)
{
  int n;
  alrandint_mersenne_twister_state *priv;

  assert(write);
  assert(state);
  assert(state->type==alrand_type_mersenne_twister);

  priv=(alrandint_mersenne_twister_state*)&state->priv;
  for (n=0;n<N;++n) write(priv->mt[n],ptr);
  write(priv->mti,ptr);
  return 0;
}

static alrand_vtable alrandint_mersenne_twister_vtable={
  &alrand_mersenne_twister_seed,
  &alrand_mersenne_twister_generate,
  &alrand_mersenne_twister_destroy,
  &alrand_mersenne_twister_save,
};

alrand_state *alrand_mersenne_twister_create(alrand_uint32_t seed)
{
  size_t size=sizeof(alrand_state)+sizeof(alrandint_mersenne_twister_state);
  alrand_state *state=malloc(size);
  if (!state) return NULL;
  state->type=alrand_type_mersenne_twister;
  state->vtable=&alrandint_mersenne_twister_vtable;
  alrand_mersenne_twister_seed(state,seed);
  return state;
}

