// 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 <cmath>
#include <cstdlib>
#include <iostream>
#include <ios>
#include <iomanip>
#include <vector>
#include <algorithm>
#include <numeric>



template <typename T> struct F_of_x
{ T operator() (T x) { return 1.0 / (1.0 + x*x); } };

template <typename T> struct RandomNumber
{ T operator() () { return static_cast<T>(std::rand())/RAND_MAX; } };

using namespace std;
int main() {
  const size_t n = 65536;
  vector<double> x(n);
  // Fill with random numbers in range [0,1]
  generate(x.begin(), x.end(), RandomNumber<double>());

  // Calculate F_of_x for each element and store back into x
  transform(x.begin(), x.end(), x.begin(), F_of_x<double>());
  double pi = 4.0 * accumulate(x.begin(), x.end(), 0.0)/n;

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