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

[Albert Chan]

Just to show expt make no difference, this version eliminated t = exp(h), replaced by cheaper sqrt().

Code:

function qthsh3(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 k, fp, fm, t0, err = 0, 0, 0, 7.38905609893065
    repeat
        local p, t = 0, sqrt(t0)
        local t1 = k==0 and t or t0
        t0 = t                      -- e, e^.5, e^.25, ...
        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=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 = (t+1/t)*r/(1+u) * (fp+fm)   -- apply weight
            p, t = p+y, t*t1
        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*ldexp(s,1-k), err/abs(s)
end

lua> f = function(x) return 1/sqrt(-log(x)) end

lua> qthsh (f, 0, 1) -- exp version
1.7724538240988943       6.586211614189097e-009
lua> qthsh2(f, 0, 1) -- expm1 version
1.7724538240988927       6.58621067462469e-009
lua> qthsh3(f, 0, 1) -- sqrt version
1.7724538240988927       6.586210486711807e-009