/* AQ.C
xref:
input:
output: stdout
does    - Globally adaptive quadrature
        - Kronrod 7-point rule (Gauss-Legendre n=3).
        - univariate integration over finite intervals
*/
/************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "pq.h"
#include "aq.h"
/***********************************************************/

static double Kronrod (double f(double), double a, double b, double* ehat);

static double split (double f(double), INTERVAL x, double* err, PQUEUE pq);

static int  conv (double x, double e, double abserr, double relerr);

/*************************************************/

double AdaptQuad (double f(double), double a, double b,
              double abserr,
              double relerr,
              double maxneval,
              int*   rc,
              int*   neval,
              double* ehat)

/*

Globally adaptive quadrature over a finite interval
Uses a 7-point Gauss-Legendre-Kronrod as the basic building block.
Input:
f       = function to be integrated
a, b    = limits of integration, a < b works properly
abserr  = absolute error tolerance
relerr  = relative error tolerance
maxneval= max acceptable # of function evaluations f(x)
Output:
returns value of integral,
*rc     <- 1 if converged, 0 otherwise,
*neval  <- # function evaluations
*ehat   <- error estimate (error := true - estimate)

*/

{

  double i,             /* i = estimate of integral*/
         e,             /* e = estimate of error   */
         v,             /* v = |e|  */
         sn=1;          /* +1 if a < b, -1 if a > b */
  INTERVAL x;
  PQUEUE pq;            /* priority queue */
  int done;
  const int k = 7; /* # evaluations per call to the core routine */

  /* handle special cases first */
  if (a == b) {*ehat = *neval = 0; *rc=1; return(0);}
  if (a >  b) {sn = -1; i=a; a=b; b=i;}
  /* now a < b, so we proceed */

  if (maxneval < k) maxneval = k;

  /* setup the queue with one interval */
  *neval = k;
  i = x.i = Kronrod (f, a, b, ehat);
  v = x.v = fabs(x.e = *ehat);
  x.a = a;
  x.b = b;

  pq = pq_construct(maxneval/k);
  pq_insert(pq, x);

  maxneval -= (k+k);    /* 2k evals each time through the loop */
  done = conv(i, v, abserr, relerr);
  while ((!done) && *neval <= maxneval) {
    x = pq_remove(pq);  /* get the interval with max abs error */
    i -= x.i;           /* subtract it out */
    *ehat -= x.e;       /* subtract the error */
    i += split(f, x, &e, pq);  /* split in two and add back */
    v = fabs(*ehat += e);  /* update error estimate    */
    done = conv(i, v, abserr, relerr);  /* converged ? */
    *neval += (k+k);            /* 2k per loop */
  }

  *rc = done;   /* success indicator */
  pq_destruct(pq);      /* release memory */

  *ehat *= sn;       /* correct the sign of the error */
  return i*sn;       /* and the integral */
}

/*************************************************/
static
double split (double f(double), INTERVAL x, double* err, PQUEUE pq)

/*
split the interval into two, evaluate each piece, then place in
the queue. Returns the sum of two pieces, *err <- error estimate.
*/

{
  double mid, b, i, le, re;

/*printf ("split: [%g, %g]\n", x.a, x.b);*/

  b = x.b;

  /* the left half */
  x.b = mid = 0.5*(x.a + x.b);
  i = x.i = Kronrod (f, x.a, mid, &le);
  x.v = fabs(x.e = le);
  pq_insert(pq, x);

  /* the right half */
  i += x.i = Kronrod (f, mid, b, &re);
  x.a = mid;
  x.b = b;
  x.v = fabs(x.e = re);
  pq_insert(pq, x);

  *err = le + re;
  return i;
}
/*************************************************/

static
int  conv (double x, double e, double abserr, double relerr)

/* returns 1 if converged, 0 otherwise */
{
  if (x != 0) return ((e <= relerr*x)  && (e <= abserr));
  else return (e <= abserr);
}
/*************************************************/

static
double Kronrod (double f(double), double a, double b, double* ehat)

/*  Kronrod's added points for n=3 Gauss-Legendre  */
/*  7 function evaluations inside (a, b), d=11 */
/*  error = 7.14e-17  (b-a)^13  f^(12)(z)         */

/*  returns estimate of integral of f() from a to b */
/*  and estimate of error */
/*  error := exact - approx  */
/*  b < a works properly */

{
  const double x1 = sqrt(0.6);
  const double x2 = 0.9604912687080202;
  const double x3 = 0.4342437493468026;
  const double w1 = 0.2684880898683334;
  const double w2 = 0.1046562260264672;
  const double w3 = 0.4013974147759622;
  const double w4 = 0.4509165386584744;

  double h, mid, integ, f1m, f1p, f0, q1, q2, sn=1;

  if (a == b) {*ehat = 0; return(0);}
  if (a >  b) {sn = -1; h=a; a=b; b=h;}

  h   = (b-a) * 0.5;            /* half width */
  mid = a + h;
  f1m = f(mid-h*x1);
  f1p = f(mid+h*x1);
  f0  = f(mid);

  q1 = (5*(f1m+f1p) + 8*f0) / 9;     /* GL, n=3 */
  q1 *= (sn * h);

  q2     =  w1*(f1m+f1p)
          + w2*(f(mid-h*x2) + f(mid+h*x2))
          + w3*(f(mid-h*x3) + f(mid+h*x3))
          + w4*f0;      /* Kronrod's 7-point formula */

  q2 *= (sn * h);

  *ehat = q2 - q1;     /* the signed error estimate */

  return q2;

}

/*************************************************/