/*
  OpenMP example program that computes the value of Pi using the trapeziod rule.
  Compile with gcc -fopenmp -O3 omp_pi.c -o omp_pi
*/
// Online source: http://users.abo.fi/mats/PP2012/examples/OpenMP/omp_pi.c
// Modified pi.c to add an orphan
#include <omp.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>

#ifdef _CIVL
$input int N=100; 
#else
#define N 1000
#endif

double starttime, d, x, sum = 0.0;
int nthreads;
int n=N;

void print_usage(char *s) {
  printf("Usage: %s -i <nr of intervals>\n", s);
  exit(0);
}

/* This is the function to integrate */
double f(double x) {
  return (4.0 / (1.0 + x*x));
}

float piOrphan(){
  /* The master thread checks how many there are */
  int i;
#pragma omp master
  {
    nthreads = omp_get_num_threads();
    printf("Number of threads = %d\n", nthreads);
    starttime = omp_get_wtime();  /* Measure the execution time */
  }

  /* This loop is executed in parallel by the threads */
#pragma omp for reduction(+:sum) 
  for (i=0; i<n; i++) {
    x = d*(double)i;
    sum += f(x) + f(x+d);
  }
}

int main (int argc, char *argv[]) {

  const double PI24 = 3.141592653589793238462643;  
  int tid;
  double pi;
  double stoptime;
  char c;

#ifndef _CIVL
  /* Check if we have at least one argument */
  if (argc <=1 ) {
    print_usage(argv[0]);
  }
  else {
    /* Parse the arguments for a -h or -i flag */
    while ((c=getopt(argc, argv, "hi:")) != EOF) {
      switch (c) {
      case 'h':
  	  print_usage(argv[0]);
      case 'i':
	  n = atoi(optarg);    /* Get number of intervals  */
	  break;
      default:
	  print_usage(argv[0]);
      }
    }
  }
#endif

  /* Compute the size of intervals */
  d = 1.0/(double) n;

  /* Start the threads */
#pragma omp parallel default(shared) private(nthreads, tid, x)
  {
    /* Get the thread number */
    tid = omp_get_thread_num();
	
	piOrphan();

  }  /* The parallel section ends here */

  /* The multiplication by d and division by 2 is moved out of the loop */  
  pi = d*sum*0.5;

  stoptime = omp_get_wtime();
  printf("The computed value of Pi is %2.24f\n", pi);
  printf("The  \"exact\" value of Pi is %2.24f\n", PI24);
  printf("The difference is %e\n", fabs(PI24-pi));
#ifdef _CIVL
  assert(fabs(PI24-pi) < 0.1);
#endif
  printf("Time: %2.4f seconds \n", stoptime-starttime);

  exit(0);
}