/* MODIFICATION TO the procedure written by Makoto Matsumoto, Takujhi Nishimura */
/* Shawn Cokus, Matthe Bellew, and Isaku Wada */

#include <stdio.h>

/* Period parameters */
#define N 624
#define M 397
#define MATRIX_A 0x9908b0dfUL		/* constant vector a */
#define UMASK 0x80000000UL		/* most significant w-r bits */
#define LMASK 0x7fffffffUL		/* least significant r bits */
#define MIXBITS(u,v) ( ((u) & UMASK) | ((v) & LMASK) )
#define TWIST(u,v) ((MIXBITS(u,v) >> 1) ^ ((v)&1UL ? MATRIX_A : 0UL))

static unsigned long state[N];		/* the array for the state vector  */
static int left = 1;
static int initf = 0;
static unsigned long *next;

/* initializes state[N] with a seed */
void    init_genrand(unsigned long s)
{
    int     j;

    state[0] = s & 0xffffffffUL;
    for (j = 1; j < N; j++) {
	state[j] = (1812433253UL * (state[j - 1] ^ (state[j - 1] >> 30)) + j);
	/* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
	/* In the previous versions, MSBs of the seed affect   */
	/* only MSBs of the array state[].                        */
	/* 2002/01/09 modified by Makoto Matsumoto             */
	state[j] &= 0xffffffffUL;		/* for >32 bit machines */
    }
    left = 1;
    initf = 1;
}

/* initialize by an array with array-length */
/* init_key is the array for initializing keys */
/* key_length is its length */
void    init_by_array(init_key, key_length)
unsigned long init_key[],
        key_length;
{
    int     i,
            j,
            k;

    init_genrand(19650218UL);
    i = 1;
    j = 0;
    k = (N > key_length ? N : key_length);
    for (; k; k--) {
	state[i] = (state[i] ^ ((state[i - 1] ^ (state[i - 1] >> 30)) * 1664525UL))
	    + init_key[j] + j;			/* non linear */
	state[i] &= 0xffffffffUL;		/* for WORDSIZE > 32 machines */
	i++;
	j++;
	if (i >= N) {
	    state[0] = state[N - 1];
	    i = 1;
	}
	if (j >= key_length)
	    j = 0;
    }
    for (k = N - 1; k; k--) {
	state[i] = (state[i] ^ ((state[i - 1] ^ (state[i - 1] >> 30)) * 1566083941UL))
	    - i;				/* non linear */
	state[i] &= 0xffffffffUL;		/* for WORDSIZE > 32 machines */
	i++;
	if (i >= N) {
	    state[0] = state[N - 1];
	    i = 1;
	}
    }

    state[0] = 0x80000000UL;			/* MSB is 1; assuring
						 * non-zero initial array */
    left = 1;
    initf = 1;
}

static void next_state(void)
{
    unsigned long *p = state;
    int     j;

    /* if init_genrand() has not been called, */
    /* a default initial seed is used         */
    if (initf == 0)
	init_genrand(5489UL);

    left = N;
    next = state;

    for (j = N - M + 1; --j; p++)
	*p = p[M] ^ TWIST(p[0], p[1]);

    for (j = M; --j; p++)
	*p = p[M - N] ^ TWIST(p[0], p[1]);

    *p = p[M - N] ^ TWIST(p[0], state[0]);
}

/* generates a random number on [0,0xffffffff]-interval */
unsigned long genrand_int32(void)
{
    unsigned long y;

    if (--left == 0)
	next_state();
    y = *next++;

    /* Tempering */
    y ^= (y >> 11);
    y ^= (y << 7) & 0x9d2c5680UL;
    y ^= (y << 15) & 0xefc60000UL;
    y ^= (y >> 18);

    return y;
}

/* generates a random number on [0,0x7fffffff]-interval */
long    genrand_int31(void)
{
    unsigned long y;

    if (--left == 0)
	next_state();
    y = *next++;

    /* Tempering */
    y ^= (y >> 11);
    y ^= (y << 7) & 0x9d2c5680UL;
    y ^= (y << 15) & 0xefc60000UL;
    y ^= (y >> 18);

    return (long) (y >> 1);
}

/* generates a random number on [0,1]-real-interval */
double  genrand_real1(void)
{
    unsigned long y;

    if (--left == 0)
	next_state();
    y = *next++;

    /* Tempering */
    y ^= (y >> 11);
    y ^= (y << 7) & 0x9d2c5680UL;
    y ^= (y << 15) & 0xefc60000UL;
    y ^= (y >> 18);

    return (double) y *(1.0 / 4294967295.0);

    /* divided by 2^32-1 */
}

/* generates a random number on [0,1)-real-interval */
double  genrand_real2(void)
{
    unsigned long y;

    if (--left == 0)
	next_state();
    y = *next++;

    /* Tempering */
    y ^= (y >> 11);
    y ^= (y << 7) & 0x9d2c5680UL;
    y ^= (y << 15) & 0xefc60000UL;
    y ^= (y >> 18);

    return (double) y *(1.0 / 4294967296.0);

    /* divided by 2^32 */
}

/* generates a random number on (0,1)-real-interval */
double  genrand_real3(void)
{
    unsigned long y;

    if (--left == 0)
	next_state();
    y = *next++;

    /* Tempering */
    y ^= (y >> 11);
    y ^= (y << 7) & 0x9d2c5680UL;
    y ^= (y << 15) & 0xefc60000UL;
    y ^= (y >> 18);

    return ((double) y + 0.5) * (1.0 / 4294967296.0);
    /* divided by 2^32 */
}

/* generates a random number on [0,1) with 53-bit resolution*/
double  genrand_res53(void)
{
    unsigned long a = genrand_int32() >> 5,
            b = genrand_int32() >> 6;

    return (a * 67108864.0 + b) * (1.0 / 9007199254740992.0);
}

/* These real versions are due to Isaku Wada, 2002/01/09 added */
/*
**************************************************************
**************************************************************
rng_afm.c
The following  subprograms, written by I. Adan, G. Fishman, and
M. Matsumoto, allow a user to perform successive Monte Carlo
experiments using nonoverlapping pseudorandom number subsequences
from the 2002/2/10 improved version of MT19937.   It also allows
a user to keep track of the initial seeds in an experiment so that
he/she may use them to rerun the experiment either after code
modification or with a variance-reducing procedure.

Executing fload_seeds(); prior to an MC experiment:

    1. Loads 624 seeds into state(1),...,state(624).

    2. Loads the value of left.

Executing fsave_seeds(); after MC execution using mt19937: 

    3. Saves the final values of state(1),...,state(624) and
       the value of left.

    4. fload_seeds() and fsave_seeds() use seed_file as input
       and output files.

    5. left is a variable in MT19937.

The first time the generator is used, 624 seeds need to be provided
as input. One source of this initalization is the statement

   init_genrand(5489UL);

Note that any integer from 1 to 2^31-1 can be used as the seed.

********************************************************************
********************************************************************
*/
void    fload_seeds(void)
{					/* loading from file seed_file */
    int     i;
    FILE   *fp;

    fp = fopen("seed_file", "r");

    for (i = 0; i < N; i++)
	fscanf(fp, "%lu\n", &state[i]);		/* loading the states */
    fscanf(fp, "%d\n", &left);			/* load the number of left
						 * words */
    initf = 1;
    /*    next = state + N - left + 1;*/		/* load the pointer to the
						 * state array */
    next = state + N - left + 1;

    fclose(fp);
}

/* *************************************************************** */
void    fsave_seeds(void)
{					/* saving to file seed_file */
    int     i;
    FILE   *fp;

    fp = fopen("seed_file", "w");

    for (i = 0; i < N; i++)
	fprintf(fp, "%lu\n", state[i]);		/* save the states */
    fprintf(fp, "%d\n", left);			/* save the number of left
						 * words */
    fclose(fp);
}

/* *************************************************************** */


int main()
{
  unsigned long int seed = 6;
  unsigned long int temp;
  int i;

  FILE *fp;

  fp = fopen("seed_file", "r");

  if (NULL == fp) {
    init_genrand(seed);
    puts("MERSENNE TWISTER INITIALIZED.");
  }
  else {
    fclose(fp);
    puts("MERSENNE TWISTER CONTINUATION.");
    fload_seeds();
  }

  for (i = 0; i < 4; i++)
    printf("idat= %i, xr= %10.9f \n", i, genrand_real2());


  fsave_seeds();

  printf("extra    xr= %10.9f\n", genrand_real2());

  return 0;
}