#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#define INITSEED 19660809UL

/*************************************************
 C code : Chi square test in 1-dimensional
  Random numbers are generated by Mersenne Twister.
**************************************************/

/* Period parameters */
#define P 624
#define Q 397
#define MATRIX_A   0x9908b0dfUL /* constant vector a */
#define UPPER_MASK 0x80000000UL /* most significant w-r bits */
#define LOWER_MASK 0x7fffffffUL /* least significant r bits */
#define nofseq 1000             /* No. of random numbers to generate */
#define k0 10                   /* No of equally spaced intervals between 0 and 1 */

static unsigned long mt[P]; /* the array for the state vector */
static int mti = P + 1;     /* mti==P+1 means mt[P] is not initialized */

/* initializes mt[P] with a seed */
void init_genrand(unsigned long s)
{
    mt[0] = s & 0xffffffffUL;
    for (mti = 1; mti < P; mti++) {
        mt[mti] = (1664525UL * mt[mti-1] + 1UL);
        mt[mti] &= 0xffffffffUL;
    }
}


/* generates a random number on [0, 0xffffffff]-interval */
unsigned long genrand(void)
{
    unsigned long y;
    static unsigned long mag01[2] = {0x0UL, MATRIX_A};
    /* mag01[x] = x * MATRIX_A for x=0, 1 */

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

        if (mti == P + 1) /* if init_genrand() has not been called, */
            init_genrand(5489UL); /* a default initial seed is used */

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

        mti = 0;
    }

    y = mt[mti++];

    /* 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_31(void)
{
    return (long)(genrand() >> 1);
}

void chi2_1dim(int n,int k,float rseq[],float *chi2,int *df,int freq[]){
	float interval[k0+1];
	int i;
	int j;
	interval[0] = 0.0;
	for (i=1; i<k+1; i++){
		freq[i] = 0;
		interval[i] = (1.0/k)*i;
	}
	*df = k -1;
	for (i=0; i<n; i++){
		for (j=1; j<k+1; j++){
			if (interval[j-1]<=rseq[i]&&rseq[i]<interval[j]){
				freq[j]++;
			}
		}
	}
	*chi2=0.0;
	for (j=1; j<k+1; j++){
		*chi2=*chi2+pow((freq[j]-(float)n/k),2)/((float)n/k);
	}
}

int main(int argc, char* argv[])
{
	long r;
	int i;

	int rmax;
	int rmin;
	int outn;
	int irseq[nofseq];
	float rseq[nofseq];
	int freq[k0+1];         /* No of equally spaced intervals between 0 and 1 */
	float chi2;
	int df;
	int k;

	k=k0;
	outn=nofseq;

	init_genrand(INITSEED);
	rmax=0;
	rmin=1000000;
	/* Convert to uniform random numbers on [0,1) */
	for (i = 0; i < outn; i++){
		r = genrand_31();
		if(rmax<r){
			rmax=r;
		}
		if(rmin>r){
			rmin=r;
		}
		irseq[i]=r;
		rseq[i]=r;
	}
	rmax=rmax+10000;
	for (i= 0; i<outn; i++){
		rseq[i]=rseq[i]/rmax;
		/* printf(fp,"%10.7f\n",rseq[i]); */
	}
	chi2_1dim(outn,k,rseq,&chi2,&df,freq);
	printf("Chi square value:%10.7f\n",chi2);
	printf("Degrees of freedom:%d\n",df);
	printf("No. of random numbers in each interval\n");
	for (i=1; i<k+1; i++){
		printf("%d\n",freq[i]);
	}

	return 0;
}