/* L-9 MCS 572 Mon 30 Jan 2012 : comptrap.c
 * The composite trapezoidal rule on sqrt(1-x^2) over [0,1] to compute pi. */

#include <stdio.h>
#include <math.h>

#define v 1      /* verbose flag */

double traprule 
 ( double (*f) ( double x ), double a, double b, int n );
/* Applies the composite Trapezoidal rule to approximate the integral
 * of f over [a,b] using n+1 function evaluations, use only for n > 0 */

double integrand ( double x );
/* Defines the function we integrate. */

int main ( int argc, char *argv[] )
{
   int n = 1000000;
   double my_pi = 0.0;
   double pi,error;

   my_pi = traprule(integrand,0.0,1.0,n);

   my_pi = 4.0*my_pi; pi = 2.0*asin(1.0); error = my_pi-pi;

   printf("Approximation for pi = %.15e with error = %.3e\n",my_pi,error);

   return 0;
}

double integrand ( double x )
{
   return sqrt(1.0 - x*x);
}

double traprule ( double (*f) ( double x ), double a, double b, int n )
{
   int i;
   double h = (b-a)/n; 
   double y = (f(a) + f(b))/2.0;
   double x;

   for(i=1,x=a+h; i < n; i++,x+=h) y += f(x);

   return h*y;
}