#include <iostream>
#include <iomanip>
#include <cstdlib>
#include <mpi.h>
#include <unistd.h>

/*
Estimate pi by Monte Carlo integration of
 pi/4 = \int_{0}^{1} = \frac{1}{1 + x^2}.

Example from:
 Pitt-Francis and Whiteley, Guide to Scientific Computing in C++, Springer 2012, ch.11

$ mpic++ -o integrate.out integrate.cpp


Result on my laptop.

Speedup    = T_1 / T_p   = [1, 1.751, 1.880, 2.171]
Efficiency = Speedup / p = [1, 0.876, 0.627, 0.543]


$ time mpirun -np 1 ./integrate.out
pi is approximated as 3.141582652!

real    0m4.107s
user    0m3.925s
sys     0m0.062s

$ time mpirun -np 2 ./integrate.out
pi is approximated as 3.141439028!

real    0m2.345s
user    0m4.302s
sys     0m0.067s

$ time mpirun -np 3 ./integrate.out
pi is approximated as 3.141557082!

real    0m2.185s
user    0m5.446s
sys     0m0.086s

$ time mpirun -np 4 ./integrate.out
pi is approximated as 3.1416283!

real    0m1.892s
user    0m6.557s
sys     0m0.136s

*/


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

  MPI::Init(argc, argv);
  int num_procs = MPI::COMM_WORLD.Get_size();
  int rank = MPI::COMM_WORLD.Get_rank();

  // set seed, make sure each process uses different one
  srand(getpid());
  long long total_n = 100000000;
  double sum = 0;

  for (long long i=0; i<total_n/num_procs; i++) {
    // generate from uniform on [0, 1]
    double x = rand() / (double)RAND_MAX;
    double f = 1.0 / (1 + x * x);
    sum += (f / total_n);
  }

  double result;
  MPI::COMM_WORLD.Reduce(&sum, &result, 1, MPI::DOUBLE, MPI::SUM, 0);
  MPI::Finalize();

  if (rank == 0) {
    std::cout << std::setprecision(10) <<
      "pi is approximated as " << 4.0*result << "!\n";
  }


  return 0;
}