// --------------------------------------------------------------------------

// RECURSIVE ADAPTIVE INTEGRATION  v. of January 20th, 2003
// (with thanks to R. Heathfield for elucidation of pointers)
// (and to Bryan Wright for elucidation of gcc's interface)

// ----------------------------------------------------
//     (c) Copyright 2003  Julian V. Noble.          //
//       Permission is granted by the author to      //
//       use this software for any application pro-  //
//       vided this copyright notice is preserved.   //
// ----------------------------------------------------


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

// prototypes
  double integral( double a , double b, double eps, double (*pn) (double) );

  double gauss3( double a, double b, double (*pn) (double) );

// Note: one can use any quadrature rule -- 3-point Gauss is not de rigeur.

// points and weights for 3 point Gauss-Legendre integration
    const double wa = 0.5555555555556 ;
    const double wb = 0.8888888888889 ;
    const double xi = 0.7745966692415 ;
// these are global constants

int main(void)
{
  typedef double (*PFI) ( double );
  PFI n;               // define a pointer to simple functions

  double a, b, eps, result;

  printf("What are a, b, and epsilon? ");
  scanf(" %lf%lf%lf", &a, &b, &eps);

  printf( "%le\n%le\n%le\n", a, b, eps ) ;

  n = sqrt;            // point to square root
  result = integral(a,b,eps, n );
  printf( " Integral with sqrt is %.12le\n", result );

  n = exp;            // point to exponential
  result = integral(a,b,eps, n );
  printf( " Integral with exp  is %.12le\n", result );
  return 0;
}


  double integral( double a, double b, double eps, double (*pn) (double) )
  // recursive
  {  
     double c, I1, I2, eps2, deltaI;

     c  = 0.5*( (a) + (b) );

     I1 = gauss3( a, b, pn );
     I2 = gauss3( a , c, pn ) + gauss3( c, b, pn );

     deltaI = fabs( I2 - I1 ) ;

     if (  deltaI < eps )
        { return I2 + ( I2 - I1)/63.0 ; }
     else
        { eps2 = 0.5*(eps);
          return integral( a, c, eps2, pn )
               + integral( c, b, eps2, pn ) ;
        }
  }


  double gauss3( double a, double b, double (*pn) (double) )
  { 
    double c, dx, xa, xb;
    c  = 0.5*( a + b ) ;
    dx = b - c ;
    xa = c - dx * xi;
    xb = c + dx * xi;
    return dx * ( wa*( (*pn) (xa) + (*pn) (xb) ) + wb * (*pn) (c) ) ; 
  }
// --------------------------------------------------------------------------