#include <iostream>
#include <random>

#include "mpi.h"

const int nb_desired_trials{100000000};

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

  // Get the number of processes
  int world_size;
  MPI_Comm_size(MPI_COMM_WORLD, &world_size);

  // Get the rank of the process
  int world_rank;
  MPI_Comm_rank(MPI_COMM_WORLD, &world_rank);

  // Print hello from each process
  std::cout << "Hello from process " << world_rank << " of "
            << world_size << std::endl;

  // Number of trials on current process
  const int nb_local_trials{nb_desired_trials / world_size};

  // Initialize random number generator
  std::mt19937_64 rng;
  rng.seed(world_rank);
  std::uniform_real_distribution<double> uniform(0, 1);

  // Count number of successful trials
  int local_successful_trials{0};
  for (int i{0}; i < nb_local_trials; ++i) {
    const double x{uniform(rng)};
    const double y{uniform(rng)};
    if (x*x + y*y < 1) {
      local_successful_trials++;
    }
  }

  // Sum over all MPI processes
  int successful_trials{0};
  int nb_trials{1};
  MPI_Reduce(static_cast<void *>(&local_successful_trials),
             static_cast<void *>(&successful_trials),
             1,            // Number of elements to send
             MPI_INTEGER,  // Data type
             MPI_SUM,      // Reduction operation
             0,            // Receive rank
             MPI_COMM_WORLD);
  // We need to compute the number of trials again here as `nb_desired_trials`
  // may not be divisible by the number of processes
  MPI_Reduce(static_cast<const void *>(&nb_local_trials),
             static_cast<void *>(&nb_trials),
             1,            // Number of elements to send
             MPI_INTEGER,  // Data type
             MPI_SUM,      // Reduction operation
             0,            // Receive rank
             MPI_COMM_WORLD);

  // Print result
  if (world_rank == 0) {
    std::cout << "pi = " << 4*static_cast<double>(successful_trials)/nb_trials
              << std::endl;
  }

  // Finalize MPI
  MPI_Finalize();
  
  return 0;
}