#include <iostream>
#include <omp.h>
#include <random>

using namespace std;

int main() {
  unsigned int seed = 0;
  default_random_engine re{seed};
  uniform_real_distribution<double> zero_to_one{0.0, 1.0};

  int n = 100000000; // number of points to generate
  int counter = 0; // counter for points lying in the first quadrant of a unit circle
  auto start_time = omp_get_wtime(); // omp_get_wtime() is an OpenMP library routine

#pragma omp parallel 
  {
    int num_threads = omp_get_num_threads();
    int thread_id = omp_get_thread_num();
    
    int thread_shift = n / num_threads; // the remainder will just get handled by the last thread
    int thread_start = thread_shift*thread_id;
    int thread_end = thread_id == num_threads-1 ? n : thread_shift*(thread_id+1);
#pragma omp critical
    { // Just a pretty lil debug print out to make sure my division logic is logicing
    std::cout << "Thread " << thread_id << " of " << num_threads << ": Computing " << thread_start << " to " << thread_end << std::endl;
    }
    
    int thread_counter = 0;
    
    // compute n points and test if they lie within the first quadrant of a unit circle
    for (int i = thread_start; i < thread_end; i++) {
      auto x = zero_to_one(re); // generate random number between 0.0 and 1.0
      auto y = zero_to_one(re); // generate random number between 0.0 and 1.0
      if (x * x + y * y <= 1.0) { // if the point lies in the first quadrant of a unit circle
        thread_counter++;
      }
    }
    
#pragma omp critical
    {
      counter = counter+thread_counter;
    }
  }

  auto run_time = omp_get_wtime() - start_time;
  auto pi = 4 * (double(counter) / n);

  cout << "pi: " << pi << endl;
  cout << "run_time: " << run_time << " s" << endl;
  cout << "n: " << n << endl; 
}