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

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

void traprule ( double (*f) ( double x ), double a, double b, int n,
                double *integral, double *error_estimate );
/* applies the composite Trapezoidal rule to approximate the integral
 * of f over [a,b] using n+1 function evaluations, use only for n > 0;
 * returns the integral with an estimate for the error */

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,err,my_pi;

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

   printf("Integral %.15e with error estimate %.2e and error %.2e\n",
          my_pi,err,my_pi-pi);

   return 0;
}

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

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

   if(n==1)
   {
      *integral = h*y;
      *error_estimate = *integral;
   }
   else
   {
      dy = 0.0;
      for(i=1,x=a+h; i < n; i++,x+=h)
         if (i % 2 == 0) 
            y += f(x);     /* approximate with step size 2*h */
         else
            dy += f(x);    /* add for finer approximation, with step h */
      
      *integral = h*(y+dy);
      *error_estimate = *integral - 2.0*h*y;
   }
}