// Normally Distributed Random Number Generator

#include <stdio.h>

#include <math.h>

#include <conio.h>

//***************************RANDOMIZER************************

#include <time.h>

unsigned long urand0 (void);

unsigned long urand (void);

void init_generator(unsigned seed);

static unsigned long rand_x[56], rand_y[256], rand_z;

static long rand_j, rand_k;

void init_generator(unsigned seed)

{

int i;

rand_x[1] = 1;

if(seed)

rand_x[2] = seed;

else

rand_x[2] = time (NULL);

for (i=3; i<56; ++i) rand_x[i] = rand_x[i-1] + rand_x[i-2];

rand_j = 24;

rand_k = 55;

for (i=255; i>=0; --i)

urand0 (); //run loop for a while

for (i=255; i>=0; --i)

rand_y[i] = urand0 ();

rand_z = urand0 ();

}

unsigned long urand0 (void)

{

if (--rand_j == 0) rand_j = 55;

if (--rand_k == 0) rand_k = 55;

return rand_x[rand_k] += rand_x[rand_j];

}

unsigned long urand (void)

{

int i;

i = rand_z % 256;

rand_z = rand_y[i];

if (--rand_j == 0) rand_j = 55;

if (--rand_k == 0) rand_k = 55;

rand_y[i] = rand_x[rand_k] += rand_x[rand_j];

return rand_z;

}

//******************************END****************************

double normal_rand (void);

//*****************************NORMAL**************************

double normal_rand (void)

{

static int flag = 0;

static double z, a = 2147483648.0;

double v1, v2, s;

if (flag) {

flag = 0;

return z;

}

flag = 1;

do {

v1 = urand()/a - 1;

v2 = urand()/a - 1;

}

while ((s = v1*v1 + v2*v2) > 1.0);

s = sqrt (-2.0 * log(s) / s);

z = v1 * s;

return v2 * s;

}

//*****************************END*****************************

double dist[19] =

{

-1.645, -1.282, -1.037, -0.842, -0.674, -0.524, -0.385, -0.253, -0.126, 0.0,

0.126, 0.253, 0.385, 0.524, 0.674, 0.842, 1.037, 1.282, 1.645

};

int main()

{

const int M=20;

const int N=1000;

const int R=5;

int i, j, k;

int a[M];

double v, y;

clrscr();

init_generator(0);

for (i=0; i<R; ++i) {

for (j=0; j<M; ++j) a[j] = 0;

for (j=0; j<N; ++j) {

y = normal_rand();

for (k=0; k<M-1; ++k)

if (y < dist[k]) break;

++a[k];

}

for (j=0; j<M; ++j) printf("%4d", a[j]);

v = 0.0;

y = N, y /= M;

for (j=0; j<M; ++j) v += (a[j] - y) * (a[j] - y) / y;

printf("\n%d degrees of freedom\n", M-1);

printf("chi-square statistic = %6.2f\n\n", v);

}

return 0;

}