// Adaptive Simpson's rule integration
//   version of February 27th, 2003 - 10:20

// ----------------------------------------------------
//     (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  ;
// these are global constants

  #define Nmax 40
  #define Nmax2 Nmax*2

  struct enchilada
  {
    double x[Nmax2], E[Nmax], f[Nmax2], Int[Nmax] ;
    double oldI, finI, newI ;
    int n, np ;
  } i ;


// prototypes
  double integral( double, double, double, double (*pn) (double) );
  void initialize( struct enchilada *, double, double, double, double (*pn) (double));
  void subdivide( struct enchilada * );
  void MoveDown( struct enchilada * );
  void NewPoints( struct enchilada *, double (*pn) (double) );
  void converged( struct enchilada *);



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


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

  void initialize( struct enchilada *i, double a, double b, double eps, double (*pn)(double) )
  {
    i->n = 2;
    i->E[0] = eps;
    i->x[0] = a;
    i->x[2] = b;
    i->x[1] = 0.5 * (a + b);
    i->f[0] = (*pn)(a);
    i->f[2] = (*pn)(b);
    i->f[1] = (*pn)( i->x[1] );
    i->Int[0] = simpson( i->f[0], i->f[1], i->f[2], i->x[2] - i->x[1] );
    i->finI = 0.0;
  }

  void subdivide( struct enchilada *i )
  {
    i->n = i->n + 2;
    if ( (i->n) > Nmax - 1 )        /* check n */
    {
        printf( " Too many subdivisions" ) ;
        abort() ;
    }
        i->np = i->n / 2 - 2 ;
        i->oldI = i->Int[ i->np ] ;
        i->E[ i->np + 1 ] = 0.5 * i->E[ i->np ] ;
  }


  void MoveDown( struct enchilada *i )
  {
    int n;
    n = i->n;

    i->x[ n ]   = i->x[ n - 2 ] ;
    i->f[ n ]   = i->f[ n - 2 ] ;
    i->x[ n - 2 ] = i->x[ n - 3 ] ;
    i->f[ n - 2 ] = i->f[ n - 3 ] ;
  }

  void NewPoints( struct enchilada *i, double (*pn) (double) )
  {
    int n, np;
    double dx;

    n  = i->n  ;
    np = i->np ;

    i->x[ n-1 ] = 0.5* ( i->x[ n ] + i->x[ n-2 ] ) ;
    i->x[ n-3 ] = 0.5* ( i->x[ n-2 ] + i->x[ n-4 ] ) ;

    i->f[ n-1 ] = (*pn) ( i->x[ n-1 ] ) ;
    i->f[ n-3 ] = (*pn) ( i->x[ n-3 ] ) ;

    dx = i->x[n]-i->x[n-1] ;
    i->Int[ np + 1 ] = simpson( i->f[n-2], i->f[n-1], i->f[n], dx );

    dx = i->x[n-2] - i->x[n-3];
    i->Int[ np ] = simpson( i->f[n-4], i->f[n-3], i->f[n-2], dx );
  }

    void converged( struct enchilada *i)
    {
        double deltaI;
        i->newI = i->Int[i->np] + i->Int[i->np + 1] ;
        deltaI = i->newI - i->oldI ;
        if ( fabs(deltaI) < i->E[i->np] )
            {
                i->finI = deltaI/15 + i->newI + i->finI ;
                i->n = i->n - 4 ;
            }
    }


  double integral( double a, double b, double eps, double (*pn) (double) )
  {
    initialize( &i, a, b, eps, pn); 
    while ( i.n > 0 )
    {
        subdivide( &i);
        MoveDown( &i);
        NewPoints( &i, pn);
        converged( &i);
    }
    return i.finI ;
  }