/*

	Author: Cyle Riggs 
	
	
	3mail:
	cyle !!at!!
	riggs family web %d0t% net

*/

#include <iostream>
#include <iomanip>
#include <cstdlib>
#include <fstream>
#include <cmath>
#include <map>
#include <vector>
#include "mpi.h"

#define TAG_POINTS 11
#define TAG_INCIRCLE 12
#define REALPI 3.14159265358979323846264338327950288419716939937510

using namespace std;

struct CPoint {
	double x;
	double y;
};

double drand( double low, double high );

int main( int argc, char* argv[] ) {
	MPI_Status status;

	MPI_Init( &argc, &argv );
	
	int size = -1;
	int rank = -1;
	MPI_Comm_size( MPI_COMM_WORLD, &size );
	MPI_Comm_rank( MPI_COMM_WORLD, &rank );
	
	srand( (rank+1) * time(NULL) + clock() );
	
	unsigned long int done = 0;
	
	double maxCount = 0;
	if( rank == 0 ) {
		cout << "max # of points? " << flush;
		cin >> maxCount;
	}
	
	MPI_Bcast( &maxCount, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD );
	
	//If the user chose 0 quit here
	if( maxCount == 0 ) {
		MPI_Finalize();
		return 0;
	}
	
	ofstream data_file;
	if( rank == 0 ) {
		data_file.open( "pi.dat", ios::in | ios::trunc );
		if( !data_file.good() ) {
			cerr << "Could not open pi.dat" << endl;
		} else {
			data_file << "#count\t#wtime\t#pi\t#error" << endl;
		}
	}
	
	for( double count = 10; count <= maxCount; count++ ) {
		//Bcast 'count' to stay in sync
		MPI_Bcast( &count, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD );
		
		double start_time = MPI_Wtime();
	
		map<unsigned long int,unsigned long int> job_map;
		for( unsigned long int i = 0; i < count; i++ ) {
			job_map[i%size]++;
		}
		
		vector<CPoint> mypoints;
		for( unsigned long int i = 0; i < job_map[rank]; i++ ) {
			CPoint p;
			p.x = drand( -1.0, 1.0 ); //Generate X for the point
			p.y = drand( -1.0, 1.0 ); //Generate Y for the point
			
			mypoints.push_back( p );
		}
		
		double my_in_circle = 0.0;
		
		for( vector<CPoint>::iterator point = mypoints.begin(); point != mypoints.end(); point++ ) {
			if( (pow(point->x,2)+pow(point->y,2)) < 1 ) {
				my_in_circle++;
			}
		}
		
		double final_in_circle = 0.0;
		
		//Add together everybody's count of points in the circle
		MPI_Reduce( &my_in_circle, &final_in_circle, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD );
		
		if( rank == 0 ) {
			//PI is the 4 * (# of points that fell in circle / the total number of points)
			double est_pi = (4.0 * (final_in_circle/count));
			
			double error = est_pi - REALPI;
			if( error < 0 ) {
				error *= -1;
			}
			
			double end_time = MPI_Wtime();
			
			data_file << setprecision(50);
			data_file << count << "\t";
			data_file << end_time - start_time << "\t";
			data_file << est_pi << "\t";
			data_file << error << "\t";
			data_file << endl << flush;
		}
	}
	
	if( rank == 0 ) {
		data_file.close();
	}
	
	MPI_Finalize();
	
	return 0;
}

double drand( double low, double high ) {
	double retval = static_cast<double>( rand() ) / RAND_MAX;
	return low + ((high - low) * retval);
}