"""
trapezoid(f, a, b, n, s = 0)
Computes the `n`th refinement of the extended trapezoid rule.
For `n=1` it returns the approximation by the trapezoid rule.
For `n>1` it improves the accuracy by adding `2^(n-2)` additional points.
`s` is the previous approximation.
"""
function trapezoid(f, a, b, n, s = 0)
if n == 1
return (f(a) + f(b)) / 2 * (b - a)
end
k = 2 ^ (n - 2) # Number of new points
h = (b - a) / k # Distance between new points (2x the distance between all points)
x = a + h / 2 # New points start here
(s + h * sum(f, x:h:b)) / 2
end
"""
simpson(f, a, b)
Computes the discrete integral of `f` between `a` and `b` using Simpson's rule.
The algorithm uses dynamic refinement with the `trapezoid` function.
"""
function simpson(f, a, b, iterations = 100, ɛ = 1.0e-7)
local s
s_prev = st_prev = -Inf
for i = 1:iterations
st = trapezoid(f, a, b, i, st_prev)
s = (4 * st - st_prev) / 3
if i > 5 && # Avoid exiting too early
(abs(s - s_prev) < ɛ * abs(s_prev) ||
(s == 0 && s_prev == 0))
break
end
s_prev, st_prev = s, st
end
s
end