#include <cassert>
#include <iostream>
#include <iomanip>
#include <sstream>
#include <pthread.h>
#include <cmath>

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;
}

// parameters of the simpson() function as structure
struct simpson_parameters {
   double (*f)(double);
   double a, b;
   int n;
   double result;
};

// invoke simpson() with the given parameters from arg
// within the current thread
extern "C" void* run_simpson_thread(void* arg) {
   struct simpson_parameters* param = (struct simpson_parameters*) arg;
   param->result = simpson(param->f, param->a, param->b, param->n);
   return param;
}

// 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) {
   assert(n > 0 && a <= b && nofthreads > 0);
   int nofintervals = n / nofthreads;
   int remainder = n % nofthreads;
   int interval = 0;
   struct simpson_parameters* params = new simpson_parameters[nofthreads];
   pthread_t* thread = new pthread_t[nofthreads];
   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;
      params[i].f = f;
      params[i].a = xleft;
      params[i].b = x;
      params[i].n = intervals;
      if (pthread_create(&thread[i], 0, run_simpson_thread, &params[i])) {
	 cerr << "unable to create threads, aborting" << endl; exit(1);
      }
   }
   double sum = 0;
   for (int i = 0; i < nofthreads; ++i) {
      struct simpson_parameters* rp;
      if (pthread_join(thread[i], (void**) &rp)) {
	 cerr << "pthread_join failed, aborting" << endl; exit(1);
      }
      sum += rp->result;
   }
   delete[] params; delete[] thread;
   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 = simpson(f, a, b, N);
   double sum = mt_simpson(f, a, b, N, nofthreads);
   cout << setprecision(14) << sum << endl;
   cout << setprecision(14) << M_PI << endl;
}