#include <cassert>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <sstream>
#include <thread>

using namespace std;

// to be integrated function
double f(double x) {
   return 4 / (1 + x*x);
}

// numerical integration according to the Simpson rule
// for f over the interval [a,b] using n subintervals
double simpson(double (*f)(double), double a, double b, int n) {
   assert(n > 0 && a <= b);
   double value = f(a)/2 + f(b)/2;
   double xleft;
   double x = a;
   for (int i = 1; i < n; ++i) {
      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;
}

// function object class representing a thread that invokes simpson
class SimpsonThread {
   public:
      SimpsonThread(double (*f)(double), double a, double b, int n,
	       double& rp) :
	    f(f), a(a), b(b), n(n), rp(rp) {
      }
      void operator()() {
	 rp = simpson(f, a, b, n);
      }
   private:
      double (*f)(double);
      double a, b;
      int n;
      double& rp;
};

// multithreaded version of simpson() that distributes the
// computation over nofthreads threads
double mt_simpson(double (*f)(double), double a, double b, int n,
      int nofthreads) {
   // divide the given interval into nofthreads partitions
   assert(n > 0 && a <= b && nofthreads > 0);
   int nofintervals = n / nofthreads;
   int remainder = n % nofthreads;
   int interval = 0;

   std::thread threads[nofthreads];
   double results[nofthreads];

   // fork off the individual threads
   double x = a;
   for (int i = 0; i < nofthreads; ++i) {
      int intervals = nofintervals;
      if (i < remainder) ++intervals;
      interval += intervals;
      double xleft = x; x = a + interval * (b - a) / n;
      threads[i] = std::thread(SimpsonThread(f,
	 xleft, x, intervals, results[i]));
   }

   // join threads and sum up their results
   double sum = 0;
   for (int i = 0; i < nofthreads; ++i) {
      threads[i].join();
      sum += results[i];
   }
   return sum;
}

char* cmdname;
void usage() {
   cerr << "Usage: " << cmdname << " [# intervals [# threads]]" << endl;
   exit(1);
}

int main(int argc, char** argv) {
   const double a = 0;
   const double b = 1;
   int N = 10000;
   int nofthreads = 10;

   cmdname = *argv++; --argc;
   if (argc > 0) {
      istringstream arg(*argv++); --argc;
      if (!(arg >> N) || N <= 0) usage();
   }
   if (argc > 0) {
      istringstream arg(*argv++); --argc;
      if (!(arg >> nofthreads) || nofthreads <= 0) usage();
   }
   if (argc > 0) usage();

   double sum = mt_simpson(f, a, b, N, nofthreads);
   cout << setprecision(14) << sum << endl;
   cout << setprecision(14) << M_PI << endl;
}