/* aquad2.c from book page 301 */

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

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

static double stack[100][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*(F0 +4.0*F1 + 6.0*F2 + 4.0*F3 + F4)/4.0; 

} /* end RS */

double aquad2(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, hmin;
  double Sab, Sac, Scb, S2ab;
  double tol, reqtol;
  double val, value, exact;

  maxtop = 99;
  hitop = 0;
  top = 0;
  value = 0.0;
  reqtol = eps;
  tol = reqtol/100.0;
  a = xmin;
  b = xmax;
  funeval = 0;
  exact = log(2.0)-log(0.1);
  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 = 1;
  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)
    {
      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("aquad2 hitop=%d, funeval=%d, hmin=%g \n", hitop, funeval, hmin);
  return value;
} /* end main of aquad2.c */