/* L-13 MCS 572 Wednesday 8 Feb 2006: the composite trapezoidal rule 
 * compile as "gcc -o /tmp/traprule0 traprule0.c -lm" */

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

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 ); /* the function we integrate */

int main ( int argc, char *argv[] )
{
   const double pi = 2.0*asin(1.0);  /* the program will approximate pi */
   int n = 100000000;                /* 100 million function evaluations */

   double y = traprule(integrand,0.0,1.0,n);   /* approximation for pi/4 */
   double my_pi = 4.0*y;

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

   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;
}