#include <cassert>
#include <mpi.h>
#include <printf.hpp>
#include <hpc/aux/slices.hpp>

using namespace std;
using namespace hpc::aux;

constexpr double a = 0;
constexpr double b = 1;
constexpr unsigned int N = 100; /* number of intervals */

// numerical integration according to the Simpson rule
// for f over the interval [a,b] using n subintervals
template <typename T, typename F>
T simpson(F f, T a, T b, int n) {
   assert(n > 0 && a <= b);
   T value = f(a)/2 + f(b)/2;
   T x = a;
   for (int i = 1; i < n; ++i) {
      T xleft = x; x = a + i * (b - a) / n;
      value += f(x) + 2 * f((xleft + x)/2);
   }
   value += 2 * f((x + b)/2);
   value *= (b - a) / n / 3;
   return value;
}

int main(int argc, char** argv) {
   MPI_Init(&argc, &argv);

   int nof_processes; MPI_Comm_size(MPI_COMM_WORLD, &nof_processes);
   int rank; MPI_Comm_rank(MPI_COMM_WORLD, &rank);

   UniformSlices<unsigned int> slices(nof_processes, N);
   auto first_interval = slices.offset(rank);
   auto nofintervals = slices.size(rank);
   auto next_interval = first_interval + nofintervals;

   double x0 = a + first_interval * (b - a) / N;
   double x1 = a + next_interval * (b - a) / N;
   double value = simpson([](double x) -> double {
	 return 4 / (1 + x*x);
   }, x0, x1, nofintervals);

   if (rank) {
      MPI_Send(&value, 1, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
   } else {
      double sum = value;
      for (int i = 0; i + 1 < nof_processes; ++i) {
	 MPI_Status status;
	 double value;
	 MPI_Recv(&value, 1, MPI_DOUBLE, MPI_ANY_SOURCE,
	    0, MPI_COMM_WORLD, &status);
	 sum += value;
      }
      fmt::printf("%.15lf\n", sum);
   }
   MPI_Finalize();
}