/* -*- Mode: C; indent-tabs-mode: t; c-basic-offset: 4; tab-width: 4 -*-  */
/*
 * mersenne_twister.c
 * Copyright (C) 2016 drdev <gualtieri@ieee.org>

 * Mersenne Twister is Copyright (C) 1997 Makoto Matsumoto and Takuji Nishimura,
 * who released it under the terms of the GNU Library General Public License,
 * version 2.                                                        
 * Ref: 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.
 *                          
 * The vast majority of this code was written by Takuji Nishimura
 * 
 * mersenne_twister 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.
 * 
 * mersenne_twister 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, see <http://www.gnu.org/licenses/>.
 */

#include<stdio.h>
#include <time.h>

/* Period parameters */  
#define N_param 624
#define M_param 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)

int r; /* random number from C library */
unsigned long r_number;
unsigned char r_byte;
int stat_bin[256];
char fn1[64] = {"MT_stats.txt"};
FILE *outdata;
static unsigned long mt[N_param]; /* the array for the state vector  */
static int mti=N_param+1; /* mti==N+1 means mt[N] is not initialized */

/* initializing the array with a NONZERO seed */
void sgenrand(seed)
    unsigned long seed;	
{
    /* setting initial seeds to mt[N] using         */
    /* the generator Line 25 of Table 1 in          */
    /* [KNUTH 1981, The Art of Computer Programming */
    /*    Vol. 2 (2nd Ed.), pp102]                  */
    mt[0]= seed & 0xffffffff;
    for (mti=1; mti<N_param; mti++)
        mt[mti] = (69069 * mt[mti-1]) & 0xffffffff;
}

unsigned long genrand()
{
    unsigned long y;
    static unsigned long mag01[2]={0x0, MATRIX_A};
    /* mag01[x] = x * MATRIX_A  for x=0,1 */

    if (mti >= N_param) { /* generate N words at one time */
	int kk;
    /* if sgenrand() has not been called, */
    /* a default initial seed is used   */
	if (mti == N_param+1) 
	{
	srand(time(0));
    	r = rand();
	printf("\n Random Number = %d\n",r);
	sgenrand(r);
        }

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

        mti = 0;
    }
  
    y = mt[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; 
}

/* main() outputs a file as an example of the uniformity of the distribution  */
main()
{ 
    long j;

    printf("\nOutput file selected = %s\n",fn1);
    if ((outdata = fopen(fn1,"wb"))==NULL)
	{printf ("\nOutput file cannot be opened.\n");
	exit (1);}

    srand(time(0));
    r = rand();
    /*printf("\n Random Number = %d\n",r);*/
    sgenrand(r); /* any nonzero integer can be used as a seed */

    /* Initialize histogram bins */

    for (j=0;j<256;j++)
    stat_bin[j]=0;


    for (j=0; j<10000000; j++)
    {
    r_number = genrand();
    r_byte = r_number%256;
	stat_bin[r_byte]++;
    r_number = r_number>>8;
    r_byte = r_number%256;
	stat_bin[r_byte]++;
    r_number = r_number>>8;
    r_byte = r_number%256;
	stat_bin[r_byte]++;
    r_number = r_number>>8;
    r_byte = r_number;
	stat_bin[r_byte]++;
    }

/* Write bins to file */
    for (j=0;j<256;j++)
    	fprintf(outdata,"%d\n", stat_bin[j]);

fclose(outdata);
printf("\nDone.\n");

    return(0);
}