/************************************************
For AMS 595 assignment, Hongyan Chen, 106508158
************************************************/
#include <cstdlib>
#include <iostream>
#include <iomanip>
#include <ctime>
#include <cmath>
using namespace std;

// the code for my own power function
long power(int base, int powerN){
     long value = 1;
     for (int k =0; k < powerN; ++k){
          value *= base;
         }
     return value;
       }
       
// the code for random area algorithm
double RanArea (unsigned long trialN){
	 unsigned long hits = 0;
	 srand ((unsigned)time(0)); 
     for (unsigned int i = 0 ; i <= trialN; ++i ) {
       double x = double(rand())/RAND_MAX;         // define x ,y as uniform random variables in [0,1]
       double y = double(rand())/RAND_MAX;
       if ( x*x + y*y < 1 ) 
		   ++hits; 
	 }
  return double(hits)*4/trialN;                    // pi/4 = p{x^2+y^2<1} 
}

// the code for geometrical approximation
double Geom (int n){
	double x = 1., y, z;
	for (int j=0; j <= n; ++j){
		x = sqrt (x+2);
	    y = sqrt (2-x);
		z = double (power (2, j+1));
	}
	return 3*y*z;                                  // pi = 3* 2^n * sqrt(2-sqrt(2+sqrt(2+sqrt(2+...sqrt(2+1)))))
}

// the code for series method. pi/2 = 1 + 1/3 + 1/3 * 2/5 + 1/3 * 2/5 * 3/7 +... 
void series (int n){
	const int ArraySize = 6000;
	int a[ArraySize], b[ArraySize];               //use arrays to store numbers with large digits, each element contains 1 digit
	memset (a, 0, sizeof(a));
	memset (b, 0, sizeof(b));
	a[0] = 2;
	b[0] = 2;
	int x = 1, y = 3, z = 0, w = 0;
	int Count = 0;
	while ( Count < (n+2)){                        //significant digits No. determines iteration cycles
		for (int i=ArraySize-1; i >= 0; --i){      //perform the multiplication from the last digit/element
			z = b[i]*x + w;
			b[i] = z % 10;
			w = z / 10;
		}
		w = 0;
		for (int j = 0; j < ArraySize; ++j){       //perform the division from the first digit/element
			z = b[j] + (w*10);
			b[j] = z / y;
			w = z % y;
		}
		w = 0;
		for (int k = ArraySize-1; k >= 0; --k){     //perform the addition from the last digit/element
			z = a[k] + b[k] +w;
			a[k] = z % 10;
			w = z / 10; 	
		}
		Count = 0;
		for (int l = 0; l < ArraySize; ++l){        //check how many significant digit already obtained
			if (b[l] == 0)
			    ++Count;
			else
				break;
		}
		++x;
		y+=2;
	}
	cout << "PI = "<< a[0] << '.' << endl;          //output pi with 30 digits each line
	for (int i = 1; i < n; ++i){
		if ((i%30) == 0)
			cout << a[i] <<endl;
		else
			cout << a[i];
	}
	
}


int main (){ 
  int MethodNumber, SigNumber;
  double PiValue;
  char flag = 'y';
  cout << "PI CALCULATION PROGRAM\nCreated by Hongyan Chen, Nov 2008.\n"<< endl;
  cout << "1. Random area algorithm (significant digits<10)" << endl;
  cout << "2. Geometrical approximation method (significant digits<18)" << endl;
  cout << "3. Series Method (100<=significant digits<=5000)" << endl;
  while (flag == 'y'){
     cout << "Please select one of the above methods (1, 2, 3):" << endl;
     do {                                                                                 //user input checking
         cin >> MethodNumber;
         if (!(MethodNumber==1||MethodNumber==2||MethodNumber==3))
            cout << "Sorry, please input 1, 2 or 3 for the method selection" << endl;
	 } while (!(MethodNumber==1||MethodNumber==2||MethodNumber==3));  
     switch (MethodNumber){
	        case 1:{
				do {                                                                      //input checking
					cout << "How many significant digits do you want?" <<endl;
                    cin >> SigNumber;
					if (SigNumber > 9 || SigNumber < 1)
						cout << "Invalid value, try again(significant digits<10)"<<endl;
				} while (SigNumber > 9 || SigNumber < 1);					      
	             unsigned long trial = 4* power(10, SigNumber-1);               
                 PiValue = RanArea(trial);
				 cout << "PI = ";
			     cout << setprecision(SigNumber) << PiValue << endl;
                 break;
			}
	        case 2:{
				do {                                                                    //input checking
					cout << "How many significant digits do you want?" <<endl;
                    cin >> SigNumber;
					if (SigNumber > 17 || SigNumber < 1)
						cout << "Invalid value, try again(significant digits<18)"<<endl;
				} while (SigNumber > 17 || SigNumber < 1);
		         double threshold = 1./power(10, SigNumber-2);
		         double pi0 = Geom(0);
		         double pi1 = Geom(1);
				 for (int i = 0; pi1-pi0 >= threshold; ++i){                             //significant digits No. determines iteration cycles
			         pi0 = Geom(i+1);
			         pi1 = Geom(i+2);
				 }
		    PiValue = pi1;
			cout << "PI = ";
			cout << setprecision(SigNumber) << PiValue << endl;
			break;
			}
	        case 3:{
				do {                                                                       //input checking
					cout << "How many significant digits do you want?" <<endl;
                    cin >> SigNumber;
					if (SigNumber > 5000 || SigNumber < 100)
						cout << "I expect a value between 100 and 5000"<<endl;
				} while (SigNumber > 5000 || SigNumber < 100);
			     series(SigNumber);
			     break;
			}
	 }
     
	 cout << "\nPress 'y' to recalculate pi, any other key to exit: ";
	 cin >> flag;                                                                   // A "getchar()? = statement" would be better,
  }                                                                                 // But I met some problems here.
  cout << "Goodbye!" << endl;
  system("pause");
  return 0;
  
}