CODEWITHSUNDEEP
cprobability
John Doe
John Doe

Reputation: 211

Calculating Birthday Probability for large numbers

The probability that two people have the same birthday in a room full of n people is 1-p. Where:

p = 365! / 365^n(365 - n)!

Obviously the numbers will be too big to solve this equation, what is a creative way to go about this?

I already solved this in a different way using simulation, but I figured the formula might be more elegant.

Upvotes: 5

Views: 2015

Answers (6)

Severin Pappadeux
Severin Pappadeux

Reputation: 20130

Holly Macaroni! What a show!

Anyway, right way to compute such things with large intermediates is to log() them

p = exp(log(p))

log(p) = log(365!) - n*log(365) - log((365 - n)!)

For factorial, use Gamma function, G(n+1)=n!, and there is very handy function in C library which computes log(G(x)): lgamma(x)

No more loops, no long doubles, no bignum libraries, no overflows...

Code

#include <math.h>
#include <stdio.h>

double b(int n) {
    double l = lgamma(365.0 + 1.0) -
               (double)n * log(365.0) -
               lgamma(365.0 - (double)n + 1.0);

    return exp(l);
}

int main() {
    double p = b(20);
    printf("%e %e\n", p, 1.0 - p);

    return 0;
}

Upvotes: 5

chux
chux

Reputation: 154602

tgamma(n+1) is very close to n!. No need to loop hundred of times which can reduce accuracy as each *, / loses a faction of a bit accuracy with each iteration.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <float.h>

long double fact(int n) {
  return roundl(tgammal(n + 1));
}

double bd_prob(int n) {
  return fact(365)/(powl(365,n)*fact(365-n));
}

int main(void){
  // No problem with 365!
  printf("fact(365) %Le\n", fact(365));
  // No problem with 365 to the 365 power 
  printf("365^365 %Le\n", powl(365, 365));

  printf("prob(22) %f\n", bd_prob(22));
  exit(EXIT_SUCCESS);
}

Output

fact(365) 2.510413e+778
365^365 1.725423e+935  
prob(22) 0.524305

Upvotes: 0

mstou
mstou

Reputation: 71

You can take advantage of 365!/(365-n)! = 365 * 364 * ... * (365-(n-1))

So to calculate this term ( let it be A=365!/(365-n)! ) you can simply the above numbers like this:

unsinged double A=1; // to make sure there is no overflow
for(int i=0;i<n;i++) A*=365-i;

To take it one step further : p=A/365^n = (364*363*...*(365-(n-1)))/365^(n-1)= 364/365 * 363/365 * ... (365-(n-1))/365.

so p can be calcuated like this:

unsigned double p=1;
for(int i=0;i<n;i++) p*= (365-i)/365.0;

in linear time

I think this should work :P

Upvotes: 3

dbush
dbush

Reputation: 225827

You don't want to calculate the full factorial. Instead, calculate each term and multiply to the result.

The probability you don't share a birthday with:

  • 1 person: 364/365
  • 2 people: 364/365 * 363/365
  • 3 people: 364/365 * 363/365 * 362/365
  • ...

Given this, you calcuate p as follows.

int n = 30;
int i;
double p = 1;
for (i = 1; i < n; i++) {
    p *= (365 - i) / 365.0;
    printf("i=%d, p=%f\n", i, 1-p);
}

Upvotes: 3

magicleon94
magicleon94

Reputation: 5192

I would write a function that looks like this:

double p(int n){
    double res = 1;
    while (n>0){
        res *= (365 - (n--))/365.0;
    }
    return res;
}

Upvotes: 0

Marcos
Marcos

Reputation: 4653

Another solution (an approximation):

The probability of any two people not having the same birthday is 364/365. In a room containing n people, there are C(n, 2) = n(n − 1)/2 pairs of people. So:

p(n) = 364/365 ^ (n * (n-1)/2)

And for values bigger than n = 100, you can safely use the next table:

n   p(n)
1   0.0%
5   2.7%
10  11.7%
20  41.1%
23  50.7%
30  70.6%
40  89.1%
50  97.0%
60  99.4%
70  99.9%
100 99.99997%
200 99.9999999999999999999999999998%
300 (100 − (6×10−80))%
350 (100 − (3×10−129))%
365 (100 − (1.45×10−155))%
366 100%
367 100%

Upvotes: 0

Related Questions