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

[Albert Chan]

qthsh version 2, using expm1() instead of exp(), for possibly more accuracy.
We can use below identity, to scale from x = expm1(a), y = expm1(b), to expm1(a+b)

(x+1)*(y+1) - 1 = x*y + x + y

Code:

expm1 = require'mathx'.expm1

function qthsh2(f, a, b, c, n, eps) -- tanh-sinh quadrature
    c = c or 1                      -- decay factor, 0 to 1
    n = n or 6                      -- 6 is optimal, 7 take longer
    eps = eps or 1e-9               -- relative error
    local s, d = f((b+a)/2), (b-a)/2
    local h, k, fp, fm, err = 2, 0, 0, 0
    repeat
        h = h/2
        local p, t = 0, expm1(h)
        local t0 = k==0 and t or t*t + 2*t
        repeat
            local u = exp(-t - t/(t+1))
            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=c*fp end
            end
            if b-x ~= b then
                y = f(b-x)
                if y-y == 0 then fm=y else fm=c*fm end
            end
            y = r/(1+u) * (fp+fm)
            y = t*t/(t+1)*y + 2*y   -- apply weight
            p, t = p+y, t0*t + t0 + t
        until abs(y) <= eps*abs(p)
        s, err, k = s+p, abs(s-p), k+1
    until err <= 10*eps*abs(s) or k>n
    return d*s*h, err/abs(s)
end

Redo example from https://www.hpmuseum.org/forum/thread-79...#pid145564

lua> function f(x) n=n+1; local y=expm1(x); return (x+x*y)/(1+y*y) end
lua> function f01(t) local y=t*t-t; return f((2*t-1)/y) * (2*y+1)/(y*y) end

lua> n = 0
lua> qthsh2(f01, 0, 1, 0)
0.8165947838638208       3.915835799526105e-009
lua> n
93

∫(f(x), x=-inf .. inf) = ∫(f01(t), t=0 .. 1) = 3/8*log(2)*pi ≈ 0.8165947838638508

f(x) will decay to 0.0 for both ends, even if intermediate calculations overflowed.

Thus, default 0.0 is appropriate here. Had we use previous point, we get this:

lua> n = 0
lua> qthsh2(f01, 0, 1, 1)
0.8165983039584289       4.3106914389618785e-006
lua> n
178

Comment: too bad expm1 version make no difference.

lua> qthsh(f01, 0, 1, 0)
0.816594783863821         3.915834575907313e-009
lua> qthsh(f01, 0, 1, 1)
0.8165983039584306       4.310691436990493e-006