(PC-12xx~14xx) qthsh Tanh-Sinh quadrature

[Albert Chan]

Lua implementation of quad(), using qthsh code to handle finite intervals.
Note: Q.quad is just a shorthand for Q["quad"], a function stored in table Q

Code:

function Q.quad(f, a, b, d, n, eps)
    d = d or 1
    n = n or 6
    eps = eps or 1e-9
    local k, c, t0, s, e = 0, 0, 7.38905609893065
    local mode = (a-a==0 and 0 or 1) + (b-b==0 and 0 or 2)

    if mode == 0 then               -- tanh-sinh
        c = d                       -- decay factor, 0 to 1
        d = (b-a)/2
        s = f((a+b)/2)
    elseif mode == 3 then           -- sinh-sinh
        s = f(c)
    else                            -- exp-sinh
        c = (mode==2) and a or b    -- c = finite edge
        if (a>b) == (mode==2) then d = -d end
        s = f(c+d)
    end

    repeat
        local p, t = 0, sqrt(t0)
        local t1 = k==0 and t or t0
        t0 = t                      -- t0 = e, e^.5, e^.25, ...
        if mode == 0 then           -- tanh-sinh
            local fp, fm = 0, 0
            repeat
                local u = exp(1/t-t)
                local r = 2*u/(1+u)
                local x, y = d*r
                if a+x ~= a then
                    y = f(a+x)
                    if y-y == 0 then fp=y else fp=fp*c end
                end
                if b-x ~= b then
                    y = f(b-x)
                    if y-y == 0 then fm=y else fm=fm*c end
                end
                y = (t+1/t)*r/(1+u) * (fp+fm)   -- apply weight
                p, t = p+y, t*t1
            until abs(y) <= eps*abs(p)
        else
            t = t * 0.5
            repeat
                local r = exp(t-0.25/t) -- = exp(sinh(j*h))
                local w, q, y = r, 0
                if mode == 3 then       -- sinh-sinh
                    r = 0.5*r - 0.5/r   -- = sinh(sinh(j*h))
                    w = 0.5*w + 0.5/w   -- = cosh(sinh(j*h))
                    y = f(c - d*r)
                    if y-y==0 then q = y*w end
                else                    -- exp-sinh
                    y = c + d/r         -- y approaching c
                    if y == c then break end
                    y = f(y)
                    if y-y==0 then q = y/w end
                end
                y = f(c + d*r)
                if y-y==0 then q = q + y*w end
                q = q * (t+0.25/t)      -- q *= cosh(j*h)
                p, t = p+q, t*t1
            until abs(q) <= eps*abs(p)
        end
        s, e, k = s+p, abs(s-p), k+1
    until e <= 10*eps*abs(s) or k>n
    if mode==1 or mode==3 and a>b then d=-d end
    return d*ldexp(s,1-k), e/abs(s)
end

---

I noticed exp-sinh is actually doing 2 integrals, simultaneously.

\(\int_a^∞ = \int_a^{a+d} + \int_{a+d}^∞ \)

Both integrals are using same number of points, but the last one is much harder.
If possible, we should pick d so that first one do most of the work.

This version allowed to input d, instead of hard-coded d=1.
If we have a clue about the shape of curve, it helps. Example, Standard Normal Distribution.

lua> Q = require 'quad'
lua> k = 1/sqrt(2*pi)
lua> Inf = math.huge
lua> pdf = function(x) return k*exp(-x*x/2) end

lua> Q.quad(Q.count(pdf), 1, Inf), Q.n
0.15865525392847651      129
lua> Q.quad(Q.count(pdf), 1, Inf, 4), Q.n
0.15865525393145705      69

Above is now as efficient as integrating closed interval of [0, 1]

lua> 1/2 - Q.quad(Q.count(pdf), 0, 1), Q.n
0.15865525393145719      58

Same idea for \(\int_{-∞}^∞\), we wanted \(\int_{-d}^d\) to do the most work.

lua> Q.quad(Q.count(pdf), -Inf, Inf), Q.n
0.999999999999999          119
lua> Q.quad(Q.count(pdf), -Inf, Inf, 4), Q.n
0.9999999999999996        43

Had we optimized specifically for even functions, above only need (n+1)/2 points !