// This program calculates pi with a Monte Carlo method i.e.
// integrate f(x) = 4/(1+x**2) on the range [0,1] by
// approximating with a sum of f(x), with N random x values
// on the interval, divided by N.

#include "myroutines.hpp"

#include <cmath>
#include <iostream>
#include <ios>
#include <iomanip>

using namespace std;


int main()
{
  // careful: x[0]....x[n-1], (stack size too...)
  size_t const n = 65536;
  double x[n];

  // use user-supplied subroutines "random_number", "sum", and "func"
  random_number(x, n);   // x vector of psuedorandom numbers in [0,1]

  double fx[n];
  func(fx, x, n);
  double pi = 4.0 * sum(fx, n) / n;





  // "atan" is an intrinsic function
  streamsize prec = cout.precision();
  cout << setprecision(8)
       << pi << '\t' << pi - 4*atan(1.0) << '\t' << M_PI
       << setprecision(prec) << endl;
}