#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <pthread.h>
#include <unistd.h>

#define MAX_WORK 500000

int Nthreads;
int myid[513];
int N;
double sum[513];
pthread_t tid[513];

double h;

double result[MAX_WORK];
int gLower[MAX_WORK];
int gUpper[MAX_WORK];
int total_work;
int cur_index = 0;

pthread_mutex_t work=PTHREAD_MUTEX_INITIALIZER;

void getWork(int *lower, int *upper, int *index)
{
	pthread_mutex_lock(&work);
	if (cur_index < total_work) {
		*lower = gLower[cur_index];
		*upper = gUpper[cur_index];
		*index = cur_index;
		cur_index = cur_index + 1;   // remove the work
	} else {
		*lower = *upper = *index = -1;
	}
	pthread_mutex_unlock(&work);
}
	

void *mycpi(void *arg) {
  int myrank;
  // int i, j;
  int i;
  double x;
  double local_sum;

  myrank = *(int *)arg;
  
  h   = 1.0 / (double) (N-1);

  result[10] = 100.0;
  
  for (;;) {
	  int lower, upper, index;

	  getWork(&lower, &upper, &index);
	  //printf("thread %d get work: %d %d %d\n", myrank, lower ,upper, index); 

	  if (lower < 0) pthread_exit(NULL);
	  local_sum = 0.0;
	 
  
	  for (i = lower; i < upper; i ++) {
		  x = h * ((double)i - 0.5);
		  local_sum += 4.0 / (1.0 + x*x);
	  }
	  result[index] = local_sum;
	  //printf("in thread result [%d] = %f\n", index, result[index]);
  }
}


int main(int argc, char *argv[])
{
  // int done = 0,  numprocs, i, rc;
  int i;
  double PI25DT = 3.141592653589793238462643;
  // double mypi, pi, h, x, a;
  double mypi, h;
 
  /*
  if (argc != 2) {
    printf("Usage: a.out N\n");
    exit(0);
  }
  */

  if (argc > 2) {
    N = atoi(argv[2]);
    Nthreads = atoi(argv[1]);
  } else if (argc > 1) {
    N = 1000;
    Nthreads = atoi(argv[1]);
  } else {N=1000; Nthreads = 1;}


  // compute work
  
  if (N < 20000) {
	  gLower[0] = 1;
	  gUpper[0] = N;
	  total_work = 1;
	  cur_index = 0;
  } else {
	  int cur;
	  gLower[0] = 1;
	  gUpper[0] = 20000;
          total_work = 1;
	  
	  cur = 20000;
	  while (cur < N) {
		  if (cur + 10000 < N) {
			  gLower[total_work] = cur;
			  gUpper[total_work] = cur+ 10000;
			  cur = cur + 10000;
			  total_work = total_work + 1;
		  } else {
			  gLower[total_work] = cur;
			  gUpper[total_work] = N;
			  cur = cur + 10000;
			  total_work = total_work + 1;
		  }			  
	  }
	  cur_index = 0;
  }
  
  for (i=0; i<Nthreads; i++) {
    myid[i] = i;
    pthread_create(&tid[i], NULL, &mycpi, &myid[i]);
  } 

  
  for (i=0; i<Nthreads; i++) {
    pthread_join(tid[i], NULL);
  }
  
  h   = 1.0 / (double) (N-1);

  mypi = 0.0;
  for (i=0; i<total_work; i++) {
    mypi += h * result[i];
    //    printf("%f result[%d] = %.16f\n", h, i, result[0]);
  }

  printf("pi is approximately %.16f, Error is %.16f\n", mypi, fabs(mypi - PI25DT));
  return 0;
}