Hey, this is awesome! I'd never heard of adaptive quadrature before. Thanks for the tip! Here's what I hacked together from your description and looking up the Newton–Cotes formulas to get a reasonable order of convergence:

    local function q(f, start, fstart, mid, fmid, stop, fstop, ε)
       -- Scaled approximation of integral on 3 points using Simpson’s
       -- rule.
       -- <https://en.wikipedia.org/wiki/Newton%E2%80%93Cotes_formulas#Closed_Newton%E2%80%93Cotes_formulas>
       local Q3 = (1/6)*(fstart + 4*fmid + fstop)

       -- Scaled approximation of integral on 5 points.
       local mid1, mid2 = (1/2)*(start + mid), (1/2)*(mid + stop)
       local fmid1, fmid2 = f(mid1), f(mid2)
       -- Boole’s rule, giving exact results for up to cubic polynomials
       -- and an O(h⁷) error.  This means that every additional level of
       -- recursion, adding 1 bit of precision to the independent
       -- variable, gives us 7 more bits of precision for smooth
       -- functions!  So we need like 256 evaluations to get results to
       -- 53-bit machine precision?
       local Q5 = (7/90)*(fstart+fstop) + (32/90)*(fmid1+fmid2) + (12/90)*fmid

       --print(mid, Q3, Q5)

       -- Approximate the error by comparing the two.
       if (stop - start)*math.abs(Q5 - Q3) <= ε then return (stop - start) * Q5 end

       -- Recursively subdivide the interval.
       return
          q(f, start, fstart, mid1, fmid1, mid, fmid, ε) +
          q(f, mid, fmid, mid2, fmid2, stop, fstop, ε)
    end

    function adaquad(f, start, stop, ε)
       local mid = (1/2)*(start + stop)
       -- In some quick testing, ε/2 seems to work okay.
       return q(f, start, f(start), mid, f(mid), stop, f(stop), (1/2)*ε)
    end

    return adaquad

http://canonical.org/~kragen/sw/dev3/adaquad.lua

It's kind of insane that 16 lines of code can do this!