// Cauchy principal value integration using
// 6-point Gauss-Legendre quadrature

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


// Description of algorithm:
//   integrate f(x)/(x - x0) from x0-delta to x0+delta
//   using even-order Gauss-Legendre quadrature formula
//


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

// prototypes
  double pval( double a , double b, double (*pn) (double) );
  double myfunc( double x );

// points and weights for 3 point Gauss-Legendre integration
    const double xi1 = 0.238619186083197;
    const double xi2 = 0.661209386466265;
    const double xi3 = 0.932469514203152;

    const double wt1 = 0.467913934572691;
    const double wt2 = 0.360761573048139;
    const double wt3 = 0.171324492379170;
// these are global constants


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

  double delta, x0, result;

  printf("What are delta and x0?");
  scanf(" %lf%lf", &delta, &x0 );

  printf( "%le\n%le\n", delta, x0 ) ;

  n = exp;            // point to exponential
  result = pval(delta, x0, n );
  printf( " Integral with exp  is %.15le\n", result );

  printf("What are delta and x0?");
  scanf(" %lf%lf", &delta, &x0 );

  printf( "%le\n%le\n", delta, x0 ) ;

  n = myfunc;         // point to myfunc
  result = pval(delta, x0, n );
  printf( " Integral with myfunc is %.15le\n", result );

  return 0;
}


double myfunc( double x )
{ return -1/( 1+x*(1+x) ) ;
}

double pval( double delta, double x0, double (*pn) (double) )
{
    double  y1, y1m, y2, y2m, y3, y3m;
    
    y1 = x0 + delta*xi1;                 // scale interval
    y2 = x0 + delta*xi2;
    y3 = x0 + delta*xi3;

    y1m = x0 - delta*xi1;
    y2m = x0 - delta*xi2;
    y3m = x0 - delta*xi3;

    return wt1 * ( (*pn) (y1) - (*pn) (y1m) ) / xi1 +
           wt2 * ( (*pn) (y2) - (*pn) (y2m) ) / xi2 +
           wt3 * ( (*pn) (y3) - (*pn) (y3m) ) / xi3 ;
}