/*
* mppi2.c
* Neil Gershenfeld 6/21/19
* OpenMP pi calculation benchmark, with multiple iterations
* pi = 3.14159265358979323846
*/

#include <stdio.h>
#include <time.h>
#include <omp.h>
#include <stdint.h>

#define NPTS 1000000000
#define NLOOP 10

void main() {
   uint64_t i;
   int j;
   double a,b,c,pi,dt,mflops,max;
   struct timespec tstart,tend;
   a = 0.5;
   b = 0.75;
   c = 0.25;
   max = 0;
   for (j = 0; j < NLOOP; ++j) {
      pi = 0;
      clock_gettime(CLOCK_REALTIME,&tstart);
      #pragma omp parallel for reduction(+:pi)
      for (i = 1; i <= NPTS; ++i)
         pi += a/((i-b)*(i-c));
      clock_gettime(CLOCK_REALTIME,&tend);
      dt = (tend.tv_sec+tend.tv_nsec/1e9)-(tstart.tv_sec+tstart.tv_nsec/1e9);
      mflops = NPTS*5.0/(dt*1e6);
      printf("NPTS = %ld, pi = %f, threads = %d\n",NPTS,pi,omp_get_max_threads());
      if (mflops > max) max = mflops;
      printf("time = %f, estimated MFlops = %f, max = %f\n",dt,mflops,max);
      }
   }