quadrature <- function(ftn, a, b, tol = 1e-8, trace = FALSE) {
  #numerical integration using adaptive quadrature
  c = (a + b)/2
  fa <- ftn(a)
  fb <- ftn(b)
  fc <- ftn(c)
  h <- (b - a)/2
  I.start <- h/3*(fa + 4*fc +fb)
  q.out <- q.recursion(ftn, a, b, c, fa, fb, fc, I.start, tol, 1, trace)
  q.out[2] <- q.out[2] + 3
  if (trace) {
    cat("final value is", q.out[1], "in", q.out[2], "function evaluations\n")
  }
  return(q.out)
}

q.recursion <- function(ftn, a, b, c, fa, fb, fc, I.old, tol, level, trace) {
#refinement of the numerical integration of ftn from a to b
  
  level.max <- 64
  if (level > level.max) {
    cat("recurion limit reached\n")
  } else {
    h <- (b-a)/4
    f1 <- ftn(a+h)
    f2 <- ftn(b-h)
    I.left <- h/3*(fa + 4*f1 + fc)
    I.right <- h/3*(fc + 4*f2 + fb)
    I.new <- I.left + I.right
    f.count <- 2
    
    if (abs(I.new - I.old) > tol) {
      q.left <- q.recursion(ftn, a, c, a+h, fa, fc, f1, I.left,
                            tol/2, level + 1, trace)
      q.right <- q.recursion(ftn, c, b, b-h, fc, fb, f2, I.right,
                            tol/2, level + 1, trace)
      I.new <- q.left[1] + q.right[1]
      f.count <- f.count + q.left[2] + q.right[2]
    } else {
      if (trace) {
        cat("integral over [", a, ",", b, "] is ", I.new,
            "(at level ", level, ")\n", sep = "" )
      }
    }
    return(c(I.new, f.count))
  }
}