#!/usr/local/bin/tcc -run
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

#define IM1  2147483563
#define IM2  2147483399
#define AM   (1.0/IM1)
#define IA1  40014
#define IA2  40692
#define IQ1  53668
#define IQ2  52774
#define IR1  12211
#define IR2  3791
#define NTAB 32
#define NDIV (1+IMM1/NTAB)
#define EPS  1.2e-7
#define RNMX (1.0-EPS)
#define NUMBER_OF_POINTS  1000

static double pi = 3.14159265358979323846;



void iid(double x[]);

main()
{
	int    i;
	int    n;
	double r1[NUMBER_OF_POINTS];
	double r2[NUMBER_OF_POINTS];
	double r3[NUMBER_OF_POINTS];
	double r4[NUMBER_OF_POINTS];
	FILE   *fp;
	time_t st;
	time_t et;

	time(&st);
	printf("\n%s\n,ctime(&st)");

	lid(x1);
	printf("\n If you want to go, press 0.\n");
	scanf("%d" , &n);
	lid(x2);

	for (n = 0; n < NUMBER_OF_POINTS ; n++){
		r1[n] = sqrt(-2.0*log(x1[n])) * cos(2.0*pi*x2[n]);
		r1[n] = sqrt(-2.0*log(x2[n])) * sin(2.0*pi*x1[n]);
	}
	
	if((fp = fopen("gauss_1.dat","w")) ==0){
		fprintf(stderr, "\n\tsorry! I cannot open gauss1.dat.\n");
	}
	for (n=0;n<NUMBER_OF_POINTS;n++)
		fprintf(fp."%f\n",r1[n]);
	fclose(fp);

	if((fp = fopen("gauss_2.dat","w")) ==0){
		fprintf(stderr, "\n\tsorry! I cannot open gauss2.dat.\n");
	}
	for (n=0;n<NUMBER_OF_POINTS;n++)
		fprintf(fp."%f\n",r2[n]);
	fclose(fp);		

	time(&et);
	printf("\n%s\n",ctime(&st));
	printf("%s\n",ctime(&et));
	printf("\n\t Gaussian random numbers are stored in gauss_1.dat and gaus_2.dat\n");
}

void iid (double[X])
{
	int   j;
	int   n;
	long  idum;
	long  k;
	static long idum2 = 123456789;
	static long iy = 0;
	static long iv[NTAB];
	double temp;

	idum = -time(NULL);

	for(n=0;n<NUMBER_OF_POINTS;n++){
		if(idum<=0){
			if (-idum<1)
				idum=1;
			else
				idem=-idum;
			idum2= idem;
			for(J=NTAB +7;j>=0;j--){
				k=idum/IQ1;
				idum = IA1*(idum-k*IQ1)-k*IR1;
				if(j<0)
					idum += IM1;
				if(j<NTAB)
					iv[j] = idum;
			}
			iy = iv[0];
		}
		
		k = idum /IQ1;
		idum = IA1*(idum-k*IQ1)- k*IR1;

		if(idum<0)
			idum += IM1;
		
		k = idum /IQ2;
		idum = IA2*(idum2-k*IQ2)- k*IR2;

		if(idum2<0)
			idum2 += IM2;
		j=iy/NDIV;
		iy=iv[j]-idum2;
		iv[j] =idum;

		if(iy<1)
			iy+=IM1;

		if ((temp = AM*iy)>RNMX)
			x[n] = RNMX;
		else
			x[n] =temp;
	}
}