/* aquad3t.c from book page 301, with modifications */

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

#undef  abs
#define abs(a)  ((a)<0.0?(-(a)):(a))

static double stack[400][7];   /* a place to store/retrieve */
static int top, hitop, maxtop; /* top points to where stored */
static int funeval=0;          /* count function evaluations */

void store(double s0, double s1, double s2, double s3, double s4,
           double s5, double s6)
{
  stack[top][0] = s0;
  stack[top][1] = s1;
  stack[top][2] = s2;
  stack[top][3] = s3;
  stack[top][4] = s4;
  stack[top][5] = s5;
  stack[top][6] = s6;
}

void retrieve(double* s0, double* s1, double* s2, double* s3, double* s4,
           double* s5, double* s6)
{
  *s0 = stack[top][0];
  *s1 = stack[top][1];
  *s2 = stack[top][2];
  *s3 = stack[top][3];
  *s4 = stack[top][4];
  *s5 = stack[top][5];
  *s6 = stack[top][6];
} /* end retrieve */

double Sn(double F0, double F1, double F2, double h)
{
  return h*(F0 + 4.0*F1 + F2)/3.0; /* 2h/6 */

} /* end Sn */

double RS(double F0, double F1, double F2, double F3,
          double F4, double h) /* 4h/16 */
{
  return h*(14.0*F0 +64.0*F1 + 24.0*F2 + 64.0*F3 + 14.0*F4)/45.0; 
  /* error term 8/945 h^7 f^(8)(c) */
} /* end RS */

double aquad3(double f(double x), double xmin, double xmax, double eps)
{
  double a, b, c, d, e;
  double Fa, Fb, Fc, Fd, Fe;
  double h1, h2, hs, hmin;
  double Sab, Sac, Scb, S2ab;
  double tol;
  double val, value;
  int i, ns; /* starting intervals */

  maxtop = 399;
  hitop = 0;
  top = 0;
  value = 0.0;
  tol = eps;
  funeval = 0;

  ns = 9;
  a = xmin;
  hs = (xmax-xmin)/(double)(ns);
  b = a + hs;
  for(i=0; i<ns; i++) /* hueristic starter set */
  {
    h1 = (b-a)/2.0;
    hmin = h1;
    c = a + h1;
    Fa = f(a); funeval++;
    Fc = f(c); funeval++;
    Fb = f(b); funeval++;
    Sab = Sn(Fa, Fc, Fb, h1);
    store(a, Fa, Fc, Fb, h1, tol, Sab);
    top++;
    a = b;
    b = a + hs;
  }
  while(top > 0)
  {
    top--;
    retrieve(&a, &Fa, &Fc, &Fb, &h1, &tol, &Sab);
    c = a + h1;
    b = a + 2.0*h1;
    h2 = h1/2;
    d = a + h2;
    e = a + 3.0*h2;
    Fd = f(d); funeval++;
    Fe = f(e); funeval++;
    Sac = Sn(Fa, Fd, Fc, h2);
    Scb = Sn(Fc, Fe, Fb, h2);
    S2ab = Sac + Scb;
    /* printf("S2ab=%15.6f, Sab=%15.6f, err=%g \n", S2ab, Sab, S2ab-Sab); */
    if(abs(S2ab-Sab) < tol || h2 < 1.0e-13)
    {
      val = RS(Fa, Fd, Fc, Fe, Fb, h2);
      value = value + val;
    }
    else
    {
      h1 = h2;
      /* printf("splitting %15.6f to %15.6f to %15.6f \n", a, c, b); */
      tol = tol/2.0;
      store(a, Fa, Fd, Fc, h1, tol, Sac);
      top++;
      store(c, Fc, Fe, Fb, h1, tol, Scb);
      top++;
      if(h1 < hmin) hmin = h1;
    }
    if(top>hitop) hitop = top;
    if(top>=maxtop) break;
  } /* end while */
  printf("aquad3t hitop=%d, funeval=%d, hmin=%g \n", hitop, funeval, hmin);
  return value;
} /* end main of aquad3t.c */