// f(x) = 4/(1+x^2)
// Domain of integration: [0,1]
// Integral over the domain = pi =(approx) 3.14159

#include <iostream>
#include <omp.h>
#include <vector>
#include <algorithm>
#include <functional>
#include <numeric>


int main (void)
{
  //Information common to serial and parallel computation.
  int    num_steps = 2e7;
  double dx        = 1.0/num_steps;



  //Serial Computation
   double start = omp_get_wtime();

   double serial_sum = 0;
   double x          = 0;
   for (int i=0;i< num_steps; ++i)
      {
         serial_sum += 4.0*dx/(1.0+x*x);
              x += dx;
     }

    double end = omp_get_wtime();
    std::cout << "Time taken for the serial computation: "      << end-start         << " seconds";
    std::cout << "\t\tPi serial: "                              << serial_sum        <<   std::endl;







   //OpenMP computation. Method of integration, just a plain Riemann sum
    std::cout << "How many OpenMP threads do you need for parallel computation? ";
    int t;//number of openmp threads
    std::cin >> t;

    start  = omp_get_wtime();
    double  parallel_sum = 0; //will be modified atomically
    #pragma omp parallel num_threads(t)
    {
      int threadIdx = omp_get_thread_num();
      int begin = threadIdx * num_steps/t; //integer index of left end point of subinterval
      int end   = begin + num_steps/t;   // interger index of right-endpoint of sub-interval
      double dx_local = dx;
      double temp = 0;
      double x    = begin*dx;

      for (int i = begin; i < end; ++i)
    {
         temp += 4.0*dx_local/(1.0+x*x);
         x    += dx_local;
    }
     #pragma omp atomic
      parallel_sum += temp;
     }
    end   = omp_get_wtime();
    std::cout << "Time taken for the parallel computation: "    << end-start << " seconds";
    std::cout << "\tPi parallel: "                                << parallel_sum        <<   std::endl;

    return 0;
}