Skip to content

Instantly share code, notes, and snippets.

@ian-abbott ian-abbott/birthdays.cpp
Last active Dec 6, 2018

Embed
What would you like to do?
Determines the probability that at least 'N' people in a group of 'M' people share the same birthday.
#include <iostream>
#include <string>
#include <gmpxx.h>
static mpz_class falling_fac(unsigned long n, unsigned long k)
{
mpz_class result = 1;
if (k > n) {
k = n;
}
while (k--) {
result *= n--;
}
return result;
}
static mpz_class binom(unsigned long n, unsigned long k)
{
mpz_class result;
mpz_bin_uiui(result.get_mpz_t(), n, k);
return result;
}
static mpz_class expon(mpz_class n, unsigned long k)
{
mpz_pow_ui(n.get_mpz_t(), n.get_mpz_t(), k);
return n;
}
static mpz_class factorial(unsigned long n)
{
mpz_class result;
mpz_fac_ui(result.get_mpz_t(), n);
return result;
}
static mpz_class k_(unsigned long n, unsigned long a, unsigned long m)
{
mpz_class result = 0;
if (n > 2) {
n--;
mpz_class fac_n_exp_k = 1;
mpz_class fac_n;
if (m / n) {
fac_n = factorial(n);
} else {
// fac_n won't be used if floor(m / n) is 0,
// so don't bother calculating it
fac_n = 1;
}
for (unsigned long k = 0; k <= m / n; k++) {
result += binom(a, k) * falling_fac(m, n * k) *
k_(n, a - k, m - n * k) / fac_n_exp_k;
fac_n_exp_k *= fac_n;
}
} else {
result = falling_fac(a, m);
}
return result;
}
static mpz_class k(unsigned long n, unsigned long a, unsigned long m)
{
if (n < 2) {
/* invalid */
return 0;
}
return k_(n, a, m);
}
static mpq_class calc_prob_q(unsigned long a, unsigned long m, unsigned long n)
{
mpz_class a_exp_m;
if (a == 0 || n > m) {
return 0;
}
if (n < 2 || n < m / a) {
return 1;
}
a_exp_m = expon(a, m);
return mpq_class(a_exp_m - k(n, a, m)) / a_exp_m;
}
int main(int argc, char *argv[])
{
unsigned long n, m, a;
mpq_class prob_q;
double prob;
if (argc < 3 || argc > 4) {
std::cerr << "usage: " << argv[0] << " N M [A]\n"
"N = minimum number with same birthday\n"
"M = number of people in group\n"
"A = number of possible birthdays (default 365)\n"
"\n"
"Determines the probability that at least 'N' people in a group of 'M' people\n"
"share the same birthday when the birthdays are linearly distributed at random\n"
"over 'A' possible birthdays. 'A' defaults to 365 if not specified.\n"
"\n"
"A limitation is that there is no account for leap years or for seasonal\n"
"variation in birth rate.\n";
return 2;
}
n = std::strtoul(argv[1], NULL, 10);
m = std::strtoul(argv[2], NULL, 10);
if (argc > 3) {
a = std::strtoul(argv[3], NULL, 10);
} else {
a = 365;
}
prob_q = calc_prob_q(a, m, n);
prob = prob_q.get_d();
std::cout << "Probability of at least one sub-group of " << n <<
" or more people out of " << m <<
"\nsharing one out of " << a << " birthdays is:\n\n" << prob <<
"\n\n" << prob_q.get_str() << "\n";
return 0;
}
@ian-abbott

This comment has been minimized.

Copy link
Owner Author

ian-abbott commented Jun 21, 2017

usage: ./birthdays N M [A]
N = minimum number with same birthday
M = number of people in group
A = number of possible birthdays (default 365)

Determines the probability that at least N people in a group of M people share the same birthday when the birthdays are linearly distributed at random over A possible birthdays. A defaults to 365 if not specified.

A limitation is that there is no account for leap years or for seasonal variation in birth rate.

Based on the answer by Nicolas56 in Probability of 3 people in a room of 30 having the same birthday.

@ian-abbott

This comment has been minimized.

Copy link
Owner Author

ian-abbott commented Jun 21, 2017

Note that the program is only practical for relatively small values due to the heavy amount of recursion involved.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.