/* 
% 
% Written by:
% -- 
% John L. Weatherwax                2005-05-26
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

Recursive trapezoidal rule (serial).

*/ 

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

/* Global storage for the integral: */
double sum = 0.0; 
/* Count the number of times we recurse: */ 
int rec_indx = 0; 

/* The function to integrate: */ 
double f( double x ){
  return 4/(1+x*x);
  /* return sin(x)*exp(x); */ 
  /* return 1.0; */
}

double 
rec_quad( double tol, double xL, double fL, double xR, double fR, double area0 ){

  double xM,fM,aL,aR;
  
  /* 
     1) Split the interval into two
     2) compute two areas with the trapezoidal rule
     3) Stop if the relative tolerance is met
  */
  rec_indx = rec_indx+1; 

  xM = 0.5*(xL+xR); fM = f(xM); 
  aL = 0.5*(xM-xL)*(fM+fL); aR = 0.5*(xR-xM)*(fM+fR); 
  printf( "rec_quad #(%d): aL=%f, aR=%f, aL+aR=%f, area0=%f\n", rec_indx, aL, aR, aL+aR, area0 );
  printf( "              : fabs(area0-(aL+aR))=%f; fabs(area0-(aL+aR))/area0=%f\n", 
	  fabs(area0-(aL+aR)), fabs(area0-(aL+aR))/area0 ); 
  if( fabs(area0-(aL+aR))/area0 <= tol ) 
    return aL+aR;
  else
    return rec_quad( tol, xL, fL, xM, fM, aL ) + rec_quad( tol, xM, fM, xR, fR, aR ); 

}

int 
main( int argc, char* argv[] ){

  double tol,xL,fL,xR,fR,area0;
  
  /* Argument checking */
  if( argc != 4 ){
    fprintf( stderr, "Expecting three input arguments\n" ); 
    exit(1); 
  }
  /* Input the left/right end points of the inteval */ 
  xL = atof(argv[1]); 
  xR = atof(argv[2]);
  /* Input the tolerance required: */ 
  tol = atof(argv[3]);

  fL = f(xL); fR = f(xR); 
  area0 = 0.5*(xR-xL)*(fL+fR); 
  printf( "main: fL=%f, fR=%f, area0=%f\n", fL, fR, area0 );
  sum = rec_quad( tol, xL, fL, xR, fR, area0 ); 

  /* Display the resulting integral: */ 
  printf("Integral = %18.10f\n", sum); 

}