// adaptive.cpp : Defines the entry point for the console application.
// This is an example of adaptive numerical quadrature 
//

#include "stdafx.h"
#include <math.h>
double f(double x) {
	/* This user-provided function defines the integrand     */
	return cos(x);
}
double integrate(double a, double b, double errTol) {
	/* This function determines an approximate integral of the  */
	/* function defined in f(x) over the range [a, b] to within */
	/* an accuracy of errTol using a total of n points or nodes */
	/* w[] defines the weights and x[] defines the nodes for    */
	/* the formula used to provide an approximate integral of f */
	int n = 5;
	double w[] = { 0.19797979797979798, 0.49090909090909090, 0.6222222222222222,
		0.49090909090909090, 0.19797979797979798 };
	double x[] = { -0.9258200997725514, -0.5773502691896257, 0.0000000000000000,
		0.5773502691896257, 0.9258200997725514 };
	double fval[5];
	double ans, errEst;
	int i;
	double sum = 0.0;
	double midpt, halfWidth;

	printf("integrate called with a= %lf, b=%lf, and errTol=%e\n", a, b, errTol);
	midpt = 0.5 * (a + b);
	halfWidth = 0.5 * (b - a);
	for (i = 0; i<n; i++) {
		fval[i] = f(midpt + halfWidth * x[i]);
		sum += w[i] * fval[i];
	}
	ans = halfWidth * sum;
	errEst = halfWidth * (fval[1] + fval[3]) - ans;

	if (fabs(errEst) > errTol) {  /* if the estimated error, errEst, is too large then */
		/* call integrate recursively over each half interval and add the results      */
		ans = integrate(a, midpt, 0.5 * errTol) + integrate(midpt, b, 0.5 * errTol);
	}
	else {  /* otherwise, the estimated error was not greater than the error tolerance */
		/* so printout the estimated error and interval, since they may be of interest */
		printf(" Estimated error is %6.2e  for interval [%6.3lf, %6.3lf]\n\n", errEst, a, b);
	}
	return ans;
}

int _tmain(int argc, _TCHAR* argv[])
{
	/* This is the main program.  User needs to provide values for:  */
	/* a                -  the left edge of the interval             */
	/* b                -  the right edge of the interval            */
	/* errTolerance     -  the desired accuracy or error tolerance   */
	/* The function to be integrated is defined in f(x)              */
	/* The function, integrate(a, b, errTolerance), returns an       */
	/* approximate value of the integral, approxAnswer               */
	double a, b, errTolerance, approxAnswer, trueAnswer;
	a = 0.0;
	b = 1.0;
	errTolerance = 1.0e-05;
	printf("Integrating from %lf to %lf with error tolerance of %6.2e\n", a, b, errTolerance);
	trueAnswer = 0.8414709848078965;
	approxAnswer = integrate(a, b, errTolerance);
	printf("Approximate answer for the integral is %20.16lf\n", approxAnswer);
	printf("Actual error is  %6.2e\n", approxAnswer - trueAnswer);
	return 0;
}