#include <cmath>
#include <iostream>
#include <mpi.h>
using namespace std;

static const double exact = 1.3932039296856768591842462603255;

double g(double* x)
{
  return 1./(1.0 - cos(x[0])*cos(x[1])*cos(x[2]))/(M_PI*M_PI*M_PI);  
}


int main(int argc, char *argv[])
{
  MPI::Init ( argc, argv );
  int mpi_size = MPI::COMM_WORLD.Get_size ( );  // How many cores?
  int my_rank = MPI::COMM_WORLD.Get_rank ( );   // Which processor is this?

  int iseed = time(0)+10*my_rank;              // Each processor should start from different random sequence.
  srand48(iseed);                              // Set iseed
  
  int N = 10000000;
  double ave=0.0;  // average
  double av2=0.0;  // average^2
  for (int i=0; i<N; i++){
    double kv[3] = {drand48()*M_PI, drand48()*M_PI, drand48()*M_PI};
    double f = g(kv);
    ave += f;
    av2 += f*f;
  }
  ave=ave/N;  // now we have average
  av2=av2/N;  //     and average^2
  double Vol = M_PI*M_PI*M_PI;  // Volume of the region
  double Int = Vol*ave;         // Integrals is Vol*<f>
  double Int2 = Vol*Vol*av2;    // For error we also need Vol^2 * <f^2>
  double err = sqrt((Int2-Int*Int)/N); // This is standard deviation  sigma^2 =  ( (<f>*Vol)^2 - <f^2>*Vol^2 )/N
  cout<<"Integral="<<Int<<" Error="<<err<<" approximation-exact="<<Int-exact<<endl;

  double res[2]={Int,Int2};
  double res_sum[2];
  //void MPI::Comm::Reduce(const void* sendbuf, void* recvbuf, int count, const MPI::Datatype& datatype, const MPI::Op& op, int root);
  MPI::COMM_WORLD.Reduce(res, res_sum, 2, MPI::DOUBLE, MPI::SUM, 0);
  
  if (my_rank==0){
    cout<<"*** Average over all processors ***"<<endl;
    double Int = res_sum[0]/mpi_size;          // Average just sums up
    double Int2 = res_sum[1]/mpi_size;         // Also average^2 just sums up
    double err = sqrt((Int2-Int*Int)/(N*mpi_size)); // New error, which is now divided by (N*mpi_size), because we have more measurements.
    cout<<"Total Int="<<Int<<" Total err="<<err<<" approximation-exact="<<Int-exact<<endl;
  }
  
  MPI::Finalize ( );
  return 0;
}