/**********************************************************
** trapezoid_OpenMP.cpp     parallelized wiht OpenMP     **
***********************************************************
** Integrate x^2 - 2x + 7 [1,10] using trapezoidal rule  **
***********************************************************
** Drive the number of subintervals into range where the **
** accumulated round-off error swamps the approximation  **
** Time the whole process and report elapsed wall clock  **
** time in anticipation of parallelizing the code.       **
**********************************************************/

#include <stdio.h>
#include <math.h>
#include <limits.h>
#include <omp.h>
#include <time.h>

const int MAX = (int) 1e9;	//maximum number of subintervals
const float X_MIN = 1.0;
const float X_MAX = 10.0;	//limits of our integration

float f(float x) {	//function to integrate
	return static_cast<float>((x-2.0)*x+7.0);	//we'll use x^2 -2x + 7
}

float integralF(float x) {//the definite integral of f(x)
	return static_cast<float>(((((1.0/3.0)*x)-1.0)*x+7.0)*x);
}

//The workhorse function: trapezoidal rule integration
float trap(float xLo, float xHi, int numIntervals) {
	float area;		//area under the curve (the integral)
	float width;	//width of each trapezoid
	float x;		//our points of evaluation
	float sum;		//summing up the f(x)'s

	sum = 0.0;		//init our summing var
	width = (xHi-xLo)/numIntervals;	//compute width of each trap

#pragma omp parallel for default(none) \
					shared(numIntervals,xLo,width) \
					private(x) \
					reduction(+:sum)

	for(int i=1; i<numIntervals; i++) {//get the interior points
		x = xLo + i*width;	//makes each iteration independent of others
		sum += f(x);		//add the interior value
	}//for

	//add the endpoints of trapezoidal rule
	sum += (f(xLo) + f(xHi))/2.0f;
	area = sum*width;

 	return area;	//return the approximation
}


/**********************************************************
*******************************  M A I N   ****************
**********************************************************/
int main () {
	int n=64;			//start with a power of 2
	float trueAnswer;	//we have the luxury of knowing the right answer
	clock_t startTime, stopTime;	//for our stopwatch function
	double elapsedTime;			//runtime for all our approximations

	trueAnswer = integralF(X_MAX)-integralF(X_MIN); //297.0;

	printf("Trapezoidal Approximation\n");
	printf("\tIntegrating x^2 - 2x + 7 from 1.0 to 10.0\n\tThe trueAnswer should be %f\n\n", trueAnswer);
	printf("numIntervals    Approx\t\tabs(Error)\n");

	startTime = clock();	//start our stopwatch
	while (n < MAX) {
		float approx = trap (X_MIN, X_MAX, n);
		printf ("%12d\t%f\t%.3E\n", n, approx, abs(trueAnswer-approx));
		n=(3*n)/2;	//increase intervals for next pass
	}
	stopTime = clock();		//stop the clock
	elapsedTime = (double)(stopTime-startTime)/CLOCKS_PER_SEC;

	printf("\n\n\nTotal running time = %.3E seconds.\n\n",elapsedTime);

	return 0;
}//main