Determines the probability that at least 'N' people in a group of 'M' people share the same birthday.

birthdays.cpp

#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;
}

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.

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