/* circlepi.c */

/* © Copyright 2005 Dauger Research, Inc. 
    
Our lawyers made us say this: 
DISCLAIMER: We provide the following on an "AS IS" basis. Use it at your own risk.  */

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

/* Produces random double-precision numbers in [0,1) */ 
inline double myrandom(); 
inline double myrandom()
	{/* You're welcome to substitute a better random number generator */
	return ((double)rand()+((double)rand())/RAND_MAX)/RAND_MAX; 
	}

inline unsigned long throwdart(double bottom, double top); 
inline unsigned long throwdart(double bottom, double top)
	{
	double y=myrandom(); 
	double x=myrandom(); 
	y=bottom+y*(top-bottom); 
	return x*x+y*y<1?1:0; 
	}
	
void performiterations(unsigned long long *numeratorOut, unsigned long long *denominatorOut, unsigned long long iterationcount, double bottom, double top); 
void performiterations(unsigned long long *numeratorOut, unsigned long long *denominatorOut, unsigned long long iterationcount, double bottom, double top)
	{
	unsigned long long numerator=0; 
	unsigned long long denominator=0; 
	unsigned long long index; 
	
	for(index=0; index<iterationcount; index++) {
		numerator+=throwdart(bottom, top); 
		numerator+=throwdart(bottom, top); 
		numerator+=throwdart(bottom, top); 
		numerator+=throwdart(bottom, top); 
		denominator++; 
		if (!(0xfffffL&denominator)) printf("%llu...\n", iterationcount-denominator); 
		}
	
	numeratorOut[0]+=numerator; 
	denominatorOut[0]+=denominator; 
	}


#define MaxCount		0x1000000LL

void performcalculation(unsigned long long *numeratorOut, unsigned long long *denominatorOut); 
void performcalculation(unsigned long long *numeratorOut, unsigned long long *denominatorOut)
	{
	unsigned long long numerator=0; 
	unsigned long long denominator=0; 
	
	performiterations(&numerator, &denominator, MaxCount, 0, 1); 
	
	numeratorOut[0]+=numerator; 
	denominatorOut[0]+=denominator; 
	}

void fprintresult(FILE *stream, unsigned long long numerator, unsigned long long denominator); 
void fprintresult(FILE *stream, unsigned long long numerator, unsigned long long denominator)
	{
	double picalculation=numerator/(double)denominator; 
	double errorestimate=sqrt((double)numerator)/(double)denominator; 
	char dashS[]="--------------------"; 
	dashS[(long)(log((double)numerator+1)/log(10))]=0; 
	fprintf(stream, "                     %llu\n", numerator); 
	fprintf(stream, "pi is approximately -%s-- = %.9g +/- %.5g\n", dashS, 
		picalculation, errorestimate); 
	fprintf(stream, "                     %llu\n", denominator); 
	}

void saveresult(char *filename, unsigned long long numerator, unsigned long long denominator); 
void saveresult(char *filename, unsigned long long numerator, unsigned long long denominator)
	{
	if (filename) {
		FILE *fp=fopen(filename, "a+"); 
		if (fp) {
			fprintresult(fp, numerator, denominator); 
			fclose(fp); 
			}
		}
	}

#include <time.h>

int main(void)
	{
	unsigned long long numerator=0; 
	unsigned long long denominator=0; 
	{
	long throwoutcount=time(NULL)&0x3ffL; 
	printf("Initializing random number generator by throwing out %d random numbers...\n", throwoutcount);
	while (throwoutcount--) myrandom(); 
	}
	printf("Beginning calculation... (first random is %g)\n", myrandom()); 
	
	performcalculation(&numerator, &denominator); 
	
	fprintresult(stdout, numerator, denominator); 
	saveresult("output", numerator, denominator); 
	
	return 0;
	}