Accurate x - log(1+x)

[Albert Chan]

(04-27-2020 07:54 PM)Albert Chan Wrote:  I tried to re-use my formula, for log(probability of no repetition)

ln_P(n,s) := -s/(12*n*(n-s)) - (s + (n-s+0.5)*ln(1-s/n))

30 people birthday "collision" ≈ -expm1(ln_P(365,30)) = 0.706316242724
Using exact formula:           1 - perm(365,30)/365^30 = 0.706316242719

For big n, small x = s/n, we can drop the first term.
Second term has a problem, even if we have accurate log1p

(s + n*ln(1 - s/n)) / n = x + log(1-x) ≈ x - x = 0

There might be a better way. For now, I use:

x + log(1-x) = x - (x + x^2/2 + x^3/3 + x^4/4 + ...) = -(x^2/2 + x^3/3 + x^4/4 + ...)

We have correct formula for ln_P(n,s), but numerically very hard to evaluate with big n.
It would be nice if we can accurately calculate x - log1p(x)

For small x, we can sum taylor/pade series. log1p(x) = 2*atanh(x/(x+2))

Code:

function tinyx_sub_log1p(x) --   x - log1p(x)
    local y = x/(x+2)       -- = x - 2*atanh(y)
    local z = y*y           -- = y*(x+2) - 2*(y + y^3/3 + y^5/5 + ...)
    z = z*(1/3 + z*(1/5 + z*(1/7 + z/9 / (1-9/11*z))))
    return y * (x-2*z)      -- relerr O(x^11/805376)
end

y = x/(x+2) ≈ x/(0+2) = x/2
abs_err(z) ≈ (1/13 - 9/121) * y^12 ≈ (4/1573) * (x/2)^12 = x^12/1610752
(x-2z) ≈ x - 2/3*y² ≈ x - 2/3*(x²/4) = x - x²/6 ≈ x

relerr = abs_err(2z)/(x-2z) ≈ x^11/805376

relerr(x=0.1000) ≈ 1.24e-17, below Free42 Binary machine epsilon
relerr(x=0.0025) ≈ 2.96e-35, below Free42 Decimal machine epsilon

For big x, we use "half-angle" formula for log1p
Keep reducing log1p argument, until we can use tinyx_sub_log1p

log1p(x) = 2 * log1p(√(1+x) - 1) = 2 * log1p(x / (√(1+x) + 1))

Let y=sqrt(1+x), z = (y-1) = x/(y+1)
Let f(x) = x - log1p(x)

f(x) = z*(y+1) - 2*log1p(z) = z*(y-1) + 2*f(z) = z^2 + 2*f(z)

Both RHS terms are non-negative, without cancellation errors.

Code:

function x_sub_log1p(x)
    if abs(x) < 0.1 then return tinyx_sub_log1p(x) end
    x = x / (sqrt(1+x)+1)   -- "half-angle" formula
    return x_sub_log1p(x)*2 + x*x
end