(32S, 32SII) Birthday paradox, hash collision probability

[Dave Britten]

The birthday paradox - called a paradox because the results seem so counterintuitive - states that in a group of just 23 people, there is approximately a 50% chance that some will share the same birthday (ignoring year), or stated another way, that the group does not have entirely unique birthdays.

This program allows you to compute the probability that there will be at least one collision when selecting S samples from N possibilities, with repetition allowed (sampling with replacement).

I won't go too deeply into the math and combinatorics behind it, as there are other resources out there that can do it more elegantly, but I'll describe the two formulas the program uses.

The simple formula, which works with smaller inputs such as birthdays, is this:

1-((N-1)/N)^((S^2-S)/2)

With large inputs, like calculating hash collisions, this formula potentially loses all of its significant digits - (N-1)/N will be rounded to 1, giving you an erroneous probability of zero.

Fortunately there's a second approximation which, while giving poor accuracy for small inputs, provides good accuracy for large ones:

1-e^(-S(S-1)/(2N))

(Reference: https://preshing.com/20110504/hash-colli...abilities/)

Of course, you may notice there's still a problem with this formula: if the x in e^x is very small, the result is 1 + an extremely small number, and the result will still be rounded to 1, giving a probability of 0. Fortunately, there's a simple way to accurately calculate e^x-1 (if you've programmed a TVM solver, you might recognize this):

e^x-1 = 2*sinh(x/2)*e^(x/2)

In order to use this transformation, we have to first restate the problem in the form of e^x-1, which is easy enough:

1-e^(-S(S-1)/(2N))
= -1*(e^(-S(S-1)/(2N))-1)
= -1*2*sinh(-S(S-1)/(2N)/2)*e^(-S(S-1)/(2N)/2)
= -2*sinh(-S(S-1)/(4N))*e^(-S(S-1)/(4N))

So this program first attempts to use the formula for small inputs, and if a result of 0 is obtained, it then applies the second formula to produce a result.

Usage

This program can either be run interactively to compute a probability for a given N and S, or via the solver to calculate N or S given the other two of N, S, and probability.

Interactive use:

XEQ B
Enter value for N, press R/S.
Enter value for S, press R/S.
Enter 0 for P (this variable is used with the solver), press R/S.

Probability will be returned in x. To see the probability expressed as "1 in x chance", press 1/X.

Solver use:

FN= B
Enter initial guesses in the X register and the variable you will be solving for (this is important, as there may be another solution to the equation with a meaningless negative value).
SOLVE, choose N S or P
Wait a couple of seconds, and results will be returned to x.

Examples:

What is the probability that at least two people in a group of 30 will share the same birthday?

XEQ B
"N?" 365 R/S
"S?" 30 R/S
"P?" 0 R/S

.6968, or about 69.68%.


What is the probability of a GUID collision if 5 trillion GUIDs are generated? (A GUID is a 128-bit globally unique identifier, with 2^128 possible values.)

XEQ B
"N?" 2 ENTER 128 y^x R/S
"S?" 5 E 12 R/S
"P?" 0 R/S

3.6734E-14, or about 1 in 27.22 trillion, i.e. not very likely!


How many GUIDs would you need to generate before there is even a 1% chance of a collision?

FN= B
2 ENTER 64 y^x
STO S (initial guess)
SOLVE S
"N?" 2 ENTER 128 y^x R/S
"P?" .01 R/S
"S=2.6153E18"

You would need to generate about 2,615,300,000,000,000,000 GUIDs.

Program Code

Code:

LBL B [CK=EB70 058.5]
INPUT N
INPUT S
INPUT P
RCL N
1/x
LASTx
1
-
*
RCL S
x^2
LASTx
-
2
/
y^x
1
x<>y
-
x!=0?
GTO C
RCL S
+/-
RCL S
1
-
*
4
/
RCL/ N
ENTER
SINH
x<>y
e^x
*
2
*
+/-

LBL C [CK=7EE2 004.5]
RCL- P
RTN

Back story: We had a data copy job at work that was failing, complaining about duplicate values in a GUID column. The table had around 5 million records in it, and I wanted to first assess whether or not there was any possibility that we legitimately had duplicate GUIDs coming from the source system. After determining that was effectively impossible, I traced the issue to some phantom reads that were happening due to a sudden increase in write activity in the source data, and changed the query's concurrency model to prevent the fake duplicate rows from appearing.

[Albert Chan]

(04-25-2020 04:49 PM)Dave Britten Wrote:  The simple formula, which works with smaller inputs such as birthdays, is this:

1-((N-1)/N)^((S^2-S)/2)

With large inputs, like calculating hash collisions, this formula potentially loses all of its significant digits - (N-1)/N will be rounded to 1, giving you an erroneous probability of zero.

Fortunately there's a second approximation which, while giving poor accuracy for small inputs, provides good accuracy for large ones:

1-e^(-S(S-1)/(2N))

(Reference: https://preshing.com/20110504/hash-colli...abilities/)

First formula can be calculated more accurately using log1p:

lua> m = require 'mathx'
lua> function p1(n,s) return -m.expm1(0.5*s*(s-1) * m.log1p(-1/n)) end
lua> function p2(n,s) return -m.expm1(0.5*s*(s-1) / -n) end
lua> function both(n,s) return p1(n,s), p2(n,s) end
lua>
lua> both(365,30)
0.6968162999539532              0.6963200177219244
lua> both(2^128, 5e12)
3.6734198463188465e-014      3.6734198463188465e-014
lua> both(2^128, 2.6153e18)
0.009999839905579181           0.009999839905579181

Probability of all different birthday, for 30 people = perm(365,30) / 365^30 ≈ 29.37%

Probability of "collision" ≈ 1 - 29.37% = 70.63%

Approximation p1() seems better than p2(), at least for test cases.

https://www.hpmuseum.org/forum/thread-11...#pid119843
https://www.hpmuseum.org/forum/thread-14...#pid126297

[Dave Britten]

(04-26-2020 07:12 PM)Albert Chan Wrote:  First formula can be calculated more accurately using log1p:

lua> m = require 'mathx'
lua> function p1(n,s) return -m.expm1(0.5*s*(s-1) * m.log1p(-1/n)) end
lua> function p2(n,s) return -m.expm1(0.5*s*(s-1) / -n) end

Very clever! And also usable on calculators that don't have a built-in ln(1+x) function:

ln(1+x) = ln(1+x)*x/(1+x-1)

[Albert Chan]

(04-26-2020 07:58 PM)Dave Britten Wrote:  Very clever! And also usable on calculators that don't have a built-in ln(1+x) function:

ln(1+x) = ln(1+x)*x/(1+x-1)

Thanks. Where do you get it ?

Note: watch out for divide-by-zero, if x is tiny.
Or, use Dieter's log1p and expm1 approximation.

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 + ...)

Code:

m = require 'mathx'
function ln_nr(n,s) -- log(probability of no repetition)
    local x = s/n   -- assumed small s, big n
    s = s - 0.5
    return -((((1/4)*x + 1/3)*x + 1/2)*x*(s-n) + s)*x
end

function p3(n,s) return -m.expm1(ln_nr(n,s)) end

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

[Dave Britten]

(04-27-2020 07:54 PM)Albert Chan Wrote:  

(04-26-2020 07:58 PM)Dave Britten Wrote:  Very clever! And also usable on calculators that don't have a built-in ln(1+x) function:

ln(1+x) = ln(1+x)*x/(1+x-1)

Thanks. Where do you get it ?

Note: watch out for divide-by-zero, if x is tiny.
Or, use Dieter's log1p and expm1 approximation.

Can't remember where it came from exactly, but it was related to TVM calculations. 15C Advanced Functions handbook maybe? And yes, if x+1-1 returns 0, then you just return x, since ln(1+x) approaches x as x gets smaller.

[Albert Chan]

(04-27-2020 07:54 PM)Albert Chan Wrote:  x + log(1-x) = x - (x + x^2/2 + x^3/3 + x^4/4 + ...) = -(x^2/2 + x^3/3 + x^4/4 + ...)

If we let y = x/(2-x) → x = 2*y - x*y

x + log(1-x) = x - 2*(y + y^3/3 + y^5/5 + ...) = -x*y - 2*(y^3/3 + y^5/5 + ...)

→ ln(P(n,s)) ≈ -s/(12*n*(n-s)) + 2*(n-s+0.5)(y^3/3 + y^5/5 + ...) - (s-1)y

I find that y^3 term already gives very good approximation.
Drop the first term, this is my revised ln_nr(n, s), shorter and more accurate.

Code:

function ln_nr(n,s)       -- log(probability of no repetition)
    local y = s/(n-s+n)   -- assumed small s, big n
    return (2/3*(n-s+0.5)*y*y - (s-1))*y
end

function P_nr(n,s)        -- return P_nr, 1 - P_nr
    local x = ln_nr(n,s)
    local y = exp(x)
    return y, 1-y+y*(log(y)-x)
end

lua> P_nr(365,30)
0.2936840560221026       0.7063159439778974
lua> P_nr(2^128,5e12)
0.9999999999999633       3.6734198463188465e-014
lua> P_nr(2^128,2.6153e18)
0.9900001600944208       0.009999839905579181

We can even drop y^3 term ...

lua> function ln_nr(n,s) return s*(s-1)/(s-n-n) end

lua> P_nr(365,30)
0.2885585859237581       0.711441414076242
lua> P_nr(2^128,5e12)
0.9999999999999633       3.6734198463188465e-014
lua> P_nr(2^128,2.6153e18)
0.9900001600944208       0.009999839905579181

Update: switched sign of y so x and y has same sign. This make it clear that x + log(1-x) ≤ 0

[Albert Chan]

(04-27-2020 07:54 PM)Albert Chan Wrote:  ln_P(n,s) := -s/(12*n*(n-s)) - (s + (n-s+0.5)*ln(1-s/n))

(04-27-2020 10:43 PM)Albert Chan Wrote:  ln(P(n,s)) ≈ -s/(12*n*(n-s)) + 2*(n-s+0.5)(y^3/3 + y^5/5 + ...) - (s-1)y

This version use the 2nd form, to avoid numerical cancellation problem of original formula.
It assumed IEEE double, ln(2^-52) ≈ -36.0, and sum y^n terms in reverse order

With all this work, we might as well add more corrections.

XCas> series(1/(e^D-1),D)
→ 1/D -1/2 +D/12 -D^3/720 +D^5/30240 +D^6*order_size(D)

Let u=n-x, f(x) = ln(n-x) = ln(u) → df = ln(n-x) dx = -ln(u) du
With x limit from 0 to s, we have u = n-x = from n to n-s

XCas> diff(-ln(u),u,3) / -720       → 1/(u^3*360)
XCas> diff(-ln(u),u,5) / 30240    → -1/(u^5*1260)

Code:

function ln_nr(n,s)         -- log(probability of no repetition)
    local c = n-s
    local y = s/(c+n)
    local y2 = y*y
    local k = floor(-36/log(y2))*2 + 3
    local t = 1/k
    for k = k-2,3,-2 do t = y2*t + 1/k end
    y2 = (c+c+1)*t*y2*y     -- = (s-1)*y - (s+(n-s+0.5)*ln(1-s/n))
    k = 1/(c*n)
    t = s*k                 -- =  1/(n-s)   - 1/n
    c = t*t + 3*k           -- = (1/(n-s)^3 - 1/n^3) / t
    k = c*(c-k) - k*k       -- = (1/(n-s)^5 - 1/n^5) / t
    return (k - 3.5*c + 105)*t*(-1/1260) + y2 - (s-1)*y
end

lua> m = require'mathx'
lua> function fac(n) return m.tgamma(n+1) end
lua> fac(69.95)/fac(69)
56.58432038352886

lua> function perm(n,s) return exp(ln_nr(n,s)) * n^s end
lua> perm(69.95, 0.95)
56.58432038352816
lua> 70*perm(69.95,-0.05) -- 69.95!/69! = 70*(69.95!/70!)
56.58432038352817
lua> 70/perm(70,0.05)       -- 69.95!/69! = 70/(70!/69.95!)
56.58432038352817

lua> ((69.95+70)/2)^0.95  -- inspired from Approximating Gamma Ratios
56.58427578530249

[Albert Chan]

Assume x ≫ a, and let \(m = (x+{1+a\over2})\), prove \(\large {(x+a)!\over x!} ≈ m^a\)

\(\large {(x+a)!\over x!}
= {m!\,/\,x! \over m!\,/\,(x+a)!}
= {perm(m, {1+a\over2}) \over perm(m, {1-a\over2})}\)

A few lines in XCas proved this ...

XCas> P_nr(n,s) := exp(-s*(s-1)/(2n))            // probability of no repetition, assumed n ≫ s
XCas> my_perm(n,s) := P_nr(n,s) * n^s         // avoid name conflict with XCas perm()
XCas> m := x+(1+a)/2
XCas> simplify(my_perm(m,(1+a)/2) / my_perm(m,(1-a)/2) - m^a)
→ 0

XCas> subst((x+a)! / x!, [x,a]=[69,0.95])
→ 56.5843203835
XCas> subst(m^a, [x,a]=[69,0.95])
→ 56.5842757853
XCas> subst(m^a, [x,a]=[70,-0.05]) * 70       // 69.95!/69! = 69.95!/70! * 70
→ 56.5843440581

Update: this ratio formula is better

XCas> ratio(x,a) := (x*x + (a+1)*(x + (a+2)/6)) ^ (a/2)
XCas> ratio(69.,0.95)
→ 56.5843203845
XCas> ratio(70., -.05) * 70                            // 69.95!/69! = 69.95!/70! * 70
→ 56.5843203829

[Albert Chan]

(05-02-2020 11:15 PM)Albert Chan Wrote:  XCas> ratio(x,a) := (x*x + (a+1)*(x + (a+2)/6)) ^ (a/2)
XCas> ratio(69.,0.95)
→ 56.5843203845
XCas> ratio(70., -.05) * 70                            // 69.95!/69! = 69.95!/70! * 70
→ 56.5843203829

We can reverse the logic, and produce a compact, accurate ln_nr(n, s)

P_nr(n, s) = perm(n, s) / n^s = ratio(n-s, s) / n^s = b^(s/2)

XCas> b := expand(ratio(n-s,s)^(2/s) / n^2)   → \(-\frac{s}{n}+\frac{s^{2}}{6\cdot n^{2}}+1+\frac{1}{n}-\frac{s}{2\cdot n^{2}}+\frac{1}{n^{2}\cdot 3}\)
XCas> factor(b-1)                                          → \(\large\frac{(-s+1) (6\cdot n-s+2)}{6\cdot n^{2}}\)

Code:

m = require 'mathx'
function ln_nr(n,s) return 0.5*s* m.log1p((s-1)*(s-2-6*n)/(6*n*n)) end
function P_collision(n,s) return -m.expm1(ln_nr(n,s)) end

lua> n, s = 365, 30
lua> x1 = ln_nr(n,s)
lua> x1, -m.expm1(x1)
-1.2252495190150705       0.706315588664643

We can improve above by halving the "interval"

lua> h = s/2
lua> x2 = ln_nr(n,h) + ln_nr(n-h,h) + h * m.log1p(-h/n)
lua> x2, -m.expm1(x2)
-1.2252516087363756       0.7063162023825731

I noticed halving the interval improved ln_nr(n,s) by about 16 times
We can extrapolate for better estimate

lua> x = x2 + (x2-x1)/15   -- (x - x1) ≈ 16*(x - x2)
lua> x, -m.expm1(x)          -- extrapolated ln_nr(n,s), collision probability
-1.2252517480511294       0.7063162432970561

Update: I found a general formula. However, convergence is not good ...

XCas> ln_nr(n,s,kmax) := ln(1+s/(n-s)) - sum(subst(sum((j/n)^k, j=1..jmax),jmax=s)/k, k=1..kmax)
XCas> map(range(1,11), k -> [k, 1 - exp(ln_nr(365.,30.,k))])

\(\left(\begin{array}{cc}
1 & 0.695232406385 \\
2 & 0.705857478651 \\
3 & 0.706293132999 \\
4 & 0.706314950579 \\
5 & 0.706316165395 \\
6 & 0.706316237864 \\
7 & 0.706316242403 \\
8 & 0.706316242698 \\
9 & 0.706316242718 \\
10 & 0.706316242719
\end{array}\right) \)

[Albert Chan]

(04-27-2020 09:16 PM)Dave Britten Wrote:  Can't remember where it [ ln(1+x) = ln(x+1)*x/(1+x-1) ] came from exactly, but it was related to TVM calculations.
15C Advanced Functions handbook maybe?

You are right !
It was proved in Theorem 4, What Every Computer Scientist Should Know About Floating-Point Arithmetic

Quote:And yes, if x+1-1 returns 0, then you just return x, since ln(1+x) approaches x as x gets smaller.

I've noticed the formula failed for IEEE double with binary exponent of 53, and odd last bit.

It is impossible to represent odd 54-bits integer x with IEEE double (53-bits mantissa).
Let "⇒" represent round-to-even rules

If x last bit is 0: (x+1)-1 ⇒ (x)-1 ⇒ x
If x last bit is 1: (x+1)-1 ⇒ (x+2)-1 ⇒ x+2

lua> log = math.log
lua> log1p = require('mathx').log1p
lua> b = 2^53
lua> for i=1,1000 do x=b+2*i; print(i, log1p(x) , log(1+x)*x/(x+1-1)) end
1       36.7368005696771       36.736800569677094
2       36.7368005696771       36.7368005696771
3       36.7368005696771       36.73680056967709
4       36.7368005696771       36.7368005696771
5       36.7368005696771       36.736800569677094
...

[Werner]

(05-03-2020 01:00 PM)Albert Chan Wrote:  

Code:

m = require 'mathx'
function ln_nr(n,s) return 0.5*s* m.log1p((s-1)*(s-2-6*n)/(6*n*n)) end
function P_collision(n,s) return -m.expm1(ln_nr(n,s)) end

Fantastic, Albert!
A 41-compatible routine, input n ENTER s:

Code:

>LBL "PCOL"
 50
 %
 DSE ST Y
 X<> ST Z    
 STO/ ST Y    
 3
 x
 STO/ ST Y
 +/-
 R^
 +
 1
 -
 x
 LN1+X
 x
 E^X-1
 +/-
 END

Cheers, Werner