#!/usr/bin/luajit
-- Adaptive quadrature for numerical integration
-- <https://en.wikipedia.org/wiki/Adaptive_quadrature> because
-- AnimalMuppet made it sound so easy in
-- <https://news.ycombinator.com/item?id=45162618>.

-- It was indeed easy but required looking up the Newton–Cotes
-- quadrature formulas for efficient high precision.

-- I’m not satisfied with the termination test.

-- CC0 <https://creativecommons.org/publicdomain/zero/1.0/>

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, math.abs(Q5 - 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