# aquad3.py
import math
stack = [[0,1,2,3,4,5,6]]   # stack to store/retrieve 

def store(s0, s1, s2, s3, s4, s5, s6):
  a = [s0, s1, s2, s3, s4, s5, s6]
  stack.append(a)

def retrieve():
  a = stack.pop()
  return a[0], a[1], a[2], a[3], a[4], a[5], a[6]

def Sn(F0, F1, F2, h):
  return h*(F0 + 4.0*F1 + F2)/3.0

def RS(F0, F1, F2, F3, F4, h):
  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)

def aquad3(f, xmin, xmax, eps):
  top = 0
  value = 0.0
  tol = eps
  ns = 32
  a = xmin
  hs = (xmax-xmin)/ns
  b = a + hs
  stack.pop() # get rid of initial junk
  for i in range(ns): # hueristic starter set
    h1 = (b-a)/2.0
    c = a + h1
    Fa = f(a)
    Fc = f(c)
    Fb = f(b)
    Sab = Sn(Fa, Fc, Fb, h1)
    store(a, Fa, Fc, Fb, h1, tol, Sab)
    top += 1
    a = b
    b = a + hs

  while(top > 0):
    top -= 1
    a, Fa, Fc, Fb, h1, tol, Sab = retrieve()
    c = a + h1
    b = a + 2.0*h1
    h2 = h1/2
    d = a + h2
    e = a + 3.0*h2
    Fd = f(d)
    Fe = f(e)
    Sac = Sn(Fa, Fd, Fc, h2)
    Scb = Sn(Fc, Fe, Fb, h2)
    S2ab = Sac + Scb
    if abs(S2ab-Sab) < tol or h2 < 1.0e-13:
      val = RS(Fa, Fd, Fc, Fe, Fb, h2)
      value += val
    else:
      h1 = h2
      tol = tol/2.0
      store(a, Fa, Fd, Fc, h1, tol, Sac)
      top += 1
      store(c, Fc, Fe, Fb, h1, tol, Scb)
      top += 1


  return value
# end aquad3.py