/* randpi.c
 * Richard Kay, April 2006
 * Computes PI based on proportion of scattered random X,Y points 
 * being within the area of a quarter circle of radius 1 within
 * a square of side 1. The area of the quarter circle is PI/4
 * and the area of the square is 1. If there is a random distribution
 * of points within the square, then PI/4 of all points should be
 * within the quarter circle. If the proportion within the quarter
 * circle diverges from PI/4 this indicates that the random data is
 * not evenly scattered. The opposite result does not follow, i.e.
 * it cannot be deduced that because the value of PI computed is
 * close to the true value of PI that the data is genuinely random,
 * but this program does give an indication. */

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

#define RFILE "random.bin" /* binary file containing random data */
#define PI 3.1415926 /* actual value to 7DP */


int getxycoord(FILE *rfh,float *x,float *y);
		/* reads next 2 random 32 bit ints from random file,
		 * converts to an x,y coordinate within a square of 
		 * side 1.0 and returns true if successful or false if EOF */

int main(void){
	FILE *rfh;
	unsigned long int inqcirc=0, counted=0, max;
	float x,y,z,myPI;
	
	printf("enter max sample size\n");
	scanf("%lu",&max);	
	if((rfh=fopen(RFILE,"r"))==NULL){
               printf("error opening random data file\n");
               return 1;
        }
	while(getxycoord(rfh,&x,&y) && counted < max){
		z=sqrt(x*x+y*y); /* pythagarus theorem */
		counted++;
		if(z<=1.0)
			inqcirc++;
	}
	myPI= 4.0*((float)inqcirc/counted);
	printf("Sample: %lu Estimate of PI: %f Actual PI: %f\n",
		counted,myPI,PI);
	printf("Estimated PI/Actual PI: %f\n",myPI/PI);
	return 0;
}

int getxycoord(FILE *rfh,float *x,float *y){

        /* reads next 2 random 32 bit int from random file. 
	 * Converts pair to a floating point X,Y coordinate
	 * with X and Y values in the range 0.0 to 1.0 . 
	 *
	 * Returns 1 if successful, 0 if file exhausted. */
	static double *twoxx32=NULL; 
	unsigned int random;
	float randf;
	if(!twoxx32){ /* only compute 2**32 on first call */
		twoxx32=(double*)malloc(sizeof(double));
		*twoxx32=pow(2,32);
	}
	/* read and convert x coordinate */
        if(feof(rfh))
        	return 0; /* end of file */ 
        fread(&random,sizeof(unsigned int),(size_t)1,rfh);
	randf=(float)random;
	*x=(float) (randf/ *twoxx32);
	/* read and convert y coordinate */
        if(feof(rfh))
        	return 0; /* end of file */ 
        fread(&random,sizeof(unsigned int),(size_t)1,rfh);
	randf=(float)random;
	*y=(float) (randf/ *twoxx32);
        return 1;
}