Accurate x - log(1+x)

[Albert Chan]

Let f(x) = x - log(1+x)

f(x) - f(-x) = (x - log(1+x)) - (-x - log(1-x)) = 2*(x - atanh(x))

LHS f(±x) never negative, thus possible massive cancellations.
We can go from accurate (atanh(x)-x) to accurate (log(1+x)-x), but not in reverse.
Code for RHS should be factored out.

Note that I flipped the order, to be consistent throughout, func_sub(x) = func(x) - x
I also use better pade approximation for atanh_sub_tiny(x), making it more useful by itself.

Code:

function atanh_sub_tiny(x)  -- abserr ≈ 256/2760615 * x^15
    local y = x*x           -- pade(atanh(x)-x,x,14,8)
    return x*y * (5005 - y*(5082 - y*969))
               / (15015 - y*(24255 - y*(11025 - y*1225)))
end

function log1p_sub_tiny(x)  -- abserr ≈ -1/176679360 * (x/(1+x/2))^15
    local y = x/(x+2)
    return atanh_sub_tiny(y)*2 - x*y    -- approx -x*y ≤ 0
end

function log1p_sub(x)       -- Free42 Decimal, 0.2 -> 0.01
    if abs(x) < 0.2 then return log1p_sub_tiny(x) end
    if x < -.56 or x > .89 then return log1p(x)-x end
    x = x / (sqrt(1+x)+1)   -- "half-angle" formula
    return log1p_sub(x)*2 - x*x
end

function expm1_sub(x)
    if x > .86 then return expm1(x) - x end
    if x < -.7 then return exp(x)-(1+x) end
    local z = expm1(-log1p_sub(x))
    return z + z*x
end

function ln_nr(n,s)         -- log(probability of no repetition)
    local x = -s/n
    local y = log1p_sub(x)
    return x/(12*(n-s)) - n*y + (s-0.5)*(x+y)
end

Example, expression ((1+r)^-n - 1)/r + n

lua> r, n = 1e-6, 1e3
lua> expm1(-n*log1p(r))/r + n --> error = 788 ULP
0.5003328749086222
lua> z = log1p_sub(r)
lua> (expm1_sub(-n*(z+r)) - n*z)/r --> error = -1 ULP
0.5003328749087098

lua> mapm = require'mapm'
lua> r = mapm.new(r) -- exact binary to decimal conversion
lua> (((1+r)^-n - 1)/r + n):tofixed(20)
0.50033287490870967825