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

[Albert Chan]

An improvement to previous post Q.peak().

Instead of calculating f(a+x)*x, this version use f(a+x), saving some computations.
Bonus: this may avoid overflow/underflow issues, if x is huge/tiny.

Code:

function Q.peak(f, a, d)    -- find peak of f(a+x)*x, guess=d
    a = a or 0
    d = d or 1
    local k, L, C, R = 2, f(a+d), f(a+d*2)
    while L/C > 2.353 do d=d/2; R=C; C=L; L=f(a+d) end
    if R then d=d*4 else d=d*2; R=C; C=L end
    while R/C > 0.425 do    -- stop if f(a+x)*x passes peak
        d = d+d
        L, C, R = C, R, f(a+d)
    end
    L, C = L/4, C/2         -- 3 points: (0,R), (1,C), (2,L)
    local d0, d1 = C-R, L-C -- finite difference, fit quadratic p(x)
    C = d0/(d1-d0) - 0.5    -- interpolate for peak, p'(x)=0, for -x
    R = R/(1.5*d0-0.5*d1)   -- extrapolate x-intercept, -x = R/p'(0)
    return d*2^C, d*2^R     -- peak, x-intercept
end

Code assumed peak is bell-shaped, not bumpy.

Also, code is optimized for peak close to user inputted d.
Example, first integral from qthsh.pdf Appendix C:

lua> f = function(x) return exp(-x)/x end
lua> Q.peak(f, 1) -- default d=1
0.5980678336446792       3.218395385172793

With default d=1, code does 3 function evaluations. (note: code optimized away last column)

Code:

N   x         f(1+x)                   f(1+x)*x
1   1         0.06766764161830635      0.06766764161830635
2   2         0.01659568945595465      0.0331913789119093
3   1/2       0.14875344009895322      0.07437672004947661

lua> Q.quad(Q.count(f), 1, huge), Q.n                -- d = 1 (default)
0.2193839343909888       129
lua> Q.quad(Q.count(f), 1, huge, 3.218), Q.n      -- d ≈ extrapolated x-intercept
0.21938393439552076     69

Update: code adjusted to quit early, if detected sign changes.