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

[Albert Chan]

My translated qthsh C code to lua

Code:

function qthsh(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, exp(h)
        local eh = k > 0 and t*t or t
        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*eh
        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

Instead of deciding settings of default 0.0 vs previous point, I added a decay factor c.
If we hit some "bad" points, it will use previous point, multiply by this decaying factor.

c=0 ≡ Default 0.0
c=1 ≡ Previous point (current default setting)

lua> function makeu(f) return function(u) return f(u*u*(3-2*u)) * 6*u*(1-u) end end
lua> uf = makeu(atanh)

lua> qthsh(uf, 0, 1, 0)       -- default 0.0
0.6931471805598393       6.9295721273908944e-012
lua> qthsh(uf, 0, 1, 1)       -- previous point
0.693147180709479         2.2281417025345338e-010

Also, limits of integration can now be swapped.

lua> qthsh(atanh, 1, 0)     -- = qthsh(atanh, 1, 0, 1, 6, 1e-9)
-0.6931471805599405      9.370022523658699e-015

---

I also removed r==0 or r==1 test.
Reaching both extreme r is very unlikely.

To hit r=0, t ≥ e^6.62, with weight underflowed to 0.0
So, even if we reached it, no harm is done.

To hit r=1, it required t=1.
This implied h must be so small that t = exp(h) = 1., or h < 2^-53.

I did some timing test on my laptop.
To reach n=18 (h=2^-18), it take about 1 second.
Each additional level take about doubled the time.

Extrapolated to n=53 (h=2^-53), it take 2^(53-18) = 2^35 sec > 1000 years ! (*)
It does not worth the code to test it.


(*) At h=2^-53, exp(h) = 1+h+h^2/2+... is slightly above half-way, thus should round-up.
However, we assume exp(x) is not "mpfr" accurate.

>>> from gmpy2 import *
>>> exp(mpfr(2)**-53)
mpfr('1.0000000000000002')