cp's OEIS Frontend

A380129 Strong Birthday Problem: Number of people needed so that probability of everyone sharing a birthday out of n possible days is at least 1/2.

%I A380129 #42 Jan 24 2025 16:28:57
%S A380129 2,4,8,12,16,21,26,31,36,41,47,52,58,64,69,75,81,87,93,100,106,112,
%T A380129 119,125,131,138,144,151,158,164,171,178,184,191,198,205,212,219,226,
%U A380129 233,240,247,254,261,268,275,283,290,297,304,312,319,326,334,341,348,356
%N A380129 Strong Birthday Problem: Number of people needed so that probability of everyone sharing a birthday out of n possible days is at least 1/2.
%C A380129 The answer to the strong birthday problem is 3064, that is, a(365) = 3064. This is the number of people that need to be gathered together before there is a 50% chance that everyone in the gathering shares their birthday with at least one other person.
%H A380129 Chai Wah Wu, <a href="/A380129/b380129.txt">Table of n, a(n) for n = 1..1000</a>
%H A380129 (Note: OEIS uses definition of m and n which are switched from original references)
%H A380129 Mario Cortina Borja, <a href="https://doi.org/10.1111/j.1740-9713.2013.00705.x">The strong birthday problem</a>, Significance 10.6 (2013): 18-20.
%H A380129 Anirban DasGupta, <a href="https://www.math.ucdavis.edu/~tracy/courses/math135A/UsefullCourseMaterial/birthday.pdf">The matching, birthday and the strong birthday problem: a contemporary review</a>, Journal of Statistical Planning and Inference 130.1-2 (2005): 377-389.
%F A380129 a(n) ~ -n LambertW(-1, -Log(2)/n).
%t A380129 (* a(n)=m, n number of days, m number of people *)
%t A380129 p[n_, m_] :=Sum[((-1)^i*n^(-m)*((-i + n)^((-i) + m))*n!*m!)/(i!*((-i + n)!)*((-i + m)!)), {i, 0, m}];
%t A380129 a[n_] := Module[{m = n + 1, prob = 0}, While[prob < 0.5, prob = p[n, m]; m++;]; m - 1];
%t A380129 Table[a[n], {n, 1, 20}]
%o A380129 (Python)
%o A380129 from math import comb, factorial
%o A380129 def A380129(n):
%o A380129     def p(m): return sum((-1 if i&1 else 1)*comb(m,i)*comb(n,i)*factorial(i)*(n-i)**(m-i) for i in range(m+1))<<1
%o A380129     kmin, kmax = n, n+1
%o A380129     while p(kmax) < n**kmax: kmax<<=1
%o A380129     while kmax-kmin > 1:
%o A380129         kmid = kmax+kmin>>1
%o A380129         if p(kmid) >= n**kmid:
%o A380129             kmax = kmid
%o A380129         else:
%o A380129             kmin = kmid
%o A380129     return kmax # _Chai Wah Wu_, Jan 21 2025
%o A380129 (PARI) p(n,m) = sum(i=0, n, (-1)^i*(n-i)^(m-i)*binomial(n,i)*m!/(m-i)!)/n^m;
%o A380129 a(n) = my(ma, mb=n+1, md); while(2*p(n,mb)<1, mb<<=1); ma=mb\2; while(mb>ma, md=(ma+mb)\2; if(2*p(n,md)<1, ma=md+1, mb=md)); ma; \\ _Jinyuan Wang_, Jan 24 2025
%Y A380129 Cf. A014088, A033810.
%K A380129 nonn
%O A380129 1,1
%A A380129 _Mike Sheppard_, Jan 13 2025