/****************************************************************
Copyright (C) :    University of Michigan
Project:           ENG101/200 - Class example
File:              randomPi.cpp
Purpose:           Calculates Pi by random sampling
Compiler:          g++
Programmer:        Alex Bielajew
Revision Date:     10/01/2103
****************************************************************/

//File: randomPi.cpp
#include <iostream>
#include <iomanip>
#include <cstdlib> // Need for rand() and srand()
#include <cmath>   // Need for pow() and atan()
#include <ctime>   // Need for time()

using namespace std;

/* Define a conversion function that converts the int returned from
   rand() to a double between -1.0 and 1.0 inclusive
*/
#define MYRAND (2.*rand()/RAND_MAX - 1.)


int main(void)
{  unsigned int nDarts  = 0; // Number of darts thrown
   unsigned int nInside = 0; // Number of darts that land inside
   int power = 1;            // Used to control output
   srand(static_cast<int>(time(NULL)));// Seed the random number generator

   double Pi = 4.0*atan(1.0); // This is how you can get a very
                              // accurate estimate of the real Pi!
   cout << "cmath gives " << Pi << " as the value of Pi\n\n";

   double xbar; // fraction that land inside
   double pi; // Estimate for Pi
   
   while(nDarts < 2147483648)
   {  double xDart = MYRAND; 
      double yDart = MYRAND;
      if (pow(xDart,2) + pow(yDart,2) <= 1.0) nInside = nInside + 1;
      nDarts = nDarts + 1;
      if (nDarts%power == 0)
      {  if (power*2 != 0) power = power*2;
         xbar = nInside/static_cast<double>(nDarts);
         pi = 4.0*xbar; // Estimate for Pi
         cout << "After " << setw(10) << nDarts
              << scientific << setprecision(6) << " P~=" << pi 
              << " d="  << setw(13) << Pi - pi 
              << "(+/-)" << sqrt(xbar*(1-xbar)/(nDarts-1)) // statisical error
              << "\n";
      }
      
   }
   nDarts = nDarts - 1; // Actual number of darts thrown
   cout << "Summary:\n";
   cout << "########\n";
   cout << "After " << setw(10) << nDarts
        << scientific << setprecision(6) << " P~=" << pi 
        << " d="  << setw(13) << Pi - pi 
        << "(+/-)" << sqrt(xbar*(1-xbar)/(nDarts-1))
        << "\n";
        
   return 0;
}