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

Let's try our new x_sub_log1p function. (note: y = x - log1p(x) ≥ 0)

Code:

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

function P_collision(n,s) return -expm1(ln_nr(n,s)) end

lua> P_collision(365, 30)
0.7063162427241915
lua> P_collision(2^128, 5e12)
3.6734198463188465e-014
lua> P_collision(2^128, 2.6153e18)
0.009999839905579181

Let's see what happen if we replace with naive implementation ...

lua> function x_sub_log1p(x) return x - log1p(x) end
lua> P_collision(365, 30) -- still OK, n is not huge
0.7063162427241918
lua> P_collision(2^128, 5e12) -- bad
7.346839692638293e-014
lua> P_collision(2^128, 2.6153e18) -- bad
0.019899683013021148