// Adaptive Simpson's rule integration
//   version of February 22nd, 2003 - 22:04

// ----------------------------------------------------
//     (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>


// weights for Simpson's Rule
    const double wa = 0.333333333333333  ;
    const double wb = 1.333333333333333  ;
    const int Nmax = 40 ;
    const int Nmax2 = 80 ;
// these are global constants

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

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( " %.10le %.10le %.10le\n", a, b, eps ) ;

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

  n = exp;            // point to exponential
  result = integral(a, b, eps, n );
  printf( " Integral with exp is %.10lf\n", result );


  return 0;
}


  #define simpson(A,B,C,D)  (D)*( wa * (A+C) + wb * (B) )

  double integral( double a, double b, double eps, double (*pn) (double) )
  {
    double x[Nmax2], E[Nmax], f[Nmax2], Int[Nmax] ;
    double oldI, finI, newI, deltaI ;
    int n, np ;

    /* initialize */

        n = 2;
        E[0] = eps;
        x[0] = a;
        x[2] = b;
        x[1] = 0.5 * (a + b);
        f[0] = (*pn)(a);
        f[2] = (*pn)(b);
        f[1] = (*pn)( x[1] );
        Int[0] = simpson( f[0], f[1], f[2], x[2]-x[1]  );
        finI = 0.0;

    while ( n>0 )       /* main loop */
    {
        /* subdivide */
            n = n + 2;
            if ( n > ( Nmax - 1 ) )     /* check n */
                {
                    printf( " Too many subdivisions" ) ;
                    return 0.0 ;
                }

            np = n/2 - 2 ;
            oldI = Int[ np ] ;
            E[ np + 1 ] = 0.5 * E[ np ] ;

            /* move down */
            x[ n ]   = x[ n-2 ] ;
            f[ n ]   = f[ n-2 ] ;
            x[ n-2 ] = x[ n-3 ] ;
            f[ n-2 ] = f[ n-3 ] ;

            /* new points */
            x[ n-1 ] = 0.5* ( x[ n ] + x[ n-2 ] ) ;
            x[ n-3 ] = 0.5* ( x[ n-2 ] + x[ n-4 ] ) ;
            f[ n-1 ] = (*pn) ( x[ n-1 ] ) ;
            f[ n-3 ] = (*pn) ( x[ n-3 ] ) ;

            Int[ np + 1 ] = simpson(f[n-2], f[n-1], f[n], x[n]-x[n-1] );
            Int[ np ] = simpson(f[n-4], f[n-3], f[n-2], x[n-2]-x[n-3] );

        /* converged? */
        newI = Int[np] + Int[np+1] ;
        deltaI = newI - oldI ;

        if ( fabs(deltaI) < E[np] )
            {
                finI = deltaI/15 + newI + finI ;
                n = n-4 ;
            }
    }

    return finI ;
  }