毎日が誕生日パーティー - 忘れても大丈夫

問題. 毎日が誕生日パーティー（勝手に命名）

1年が365日の世界に， $n$ 人の社員がいる．各社員の誕生日は365日の間に等確率で分布する．このとき，どの日にも少なくとも1人の社員が誕生日となる確率を求めたい（誕生日の人がいるとパーティー）．
　(1) $n$ が 365 のときの確率を求めよ
　(2) 確率が 0.5 を超える最小の社員数を求めよ

背景

Twitter上で @mikupaccho さんの次のツイートを見かけたので解いてみました．問題の解法自体は難しくないのですがC++ Boostの多倍長整数を初めて使用したので記事を書きました．

何気ない私の疑問が東大生インターンをきずつけた pic.twitter.com/Wr96JGtfjm
— Miku Nozaki (@mikupaccho) 2018年8月6日

解法. 全射の総数

社員の集合を $E = \{1, \ldots, n\}$ ，誕生日の集合を $D = \{1, \ldots, 365\}$ として， $E$ から $D$ への写像の総数を $f(n)$ ，全射の総数を $g(n)$ とすると，社員数 $n$ のときの求める確率 $P(n)$ は，
　 $P(n) = \frac{g(n)}{f(n)}$
となります．写像の総数は，各社員に対して誕生日の選び方が365通りあるので $f(n) = 365^n$ となります．一方，全射の総数は包除原理を用いると， $g(n) = \sum_{i = 1}^{365} (-1)^{365 - i} \binom{365}{i} i^n$ となります．全射の総数の詳しい求め方は全射の個数の証明とベル数 | 高校数学の美しい物語を参照してください．
以上から求める確率は，
　 $P(n) = \frac{g(n)}{f(n)} = \frac{\sum_{i = 1}^{365} (-1)^{365 - i} \binom{365}{i} i^n}{365^n}$
です．

問題(1)

問題(1)も上の式に当てはめて計算しても良いのですが，社員数が誕生日の数と等しいので，同じ誕生日の人が2人以上いると全射になりません．したがって全射の数は置換の数と等しく $g(365) = 365!$ となり， $P(365) = \frac{365!}{365^{365}} \sim 1.4549 × 10^{-157}$ となります．計算は Wolfram|Alpha: Computational Intelligence で行いました．

問題(2)

問題(2) は $P(n) \ge 0.5$ となる最小の $n$ を求める問題となります．浮動小数点数で計算しても良いのですが数値誤差がコワいので式変形をして $2 \times g(n) \ge f(n)$ となる最小の $n$ を求めました．答えは，社員数が2287人以上いれば確率0.5以上で毎日がパーティーになります．

C++のSTLには多倍長整数はサポートされていませんが，Boostにはサポートされています．詳しい説明は多倍長整数 - boostjp を参照してください．
Boostの多倍長整数にはBoost独自実装，GMPバックエンドとlibtommathバックエンドの3つの多倍長整数を使用することができます．これらを使用したときとRubyで実装したときの， $n$ を 365からインクリメントしながら上の式を計算したときの実行時間は次の通りです（1回のみの測定）．アルゴリズムの改良は出来るかも知れませんが愚直に書いています．

実装方法	実行時間
Ruby	58.546 秒
C++ Boost独自実装	58.356 秒
C++ GMPバックエンド	24.989 秒
C++ tommathバックエンド	4 分 52 秒

C++のソースコードを表示

#include <bits/stdc++.h>
#include <boost/multiprecision/cpp_int.hpp> // boost独自の実装
// #include <boost/multiprecision/gmp.hpp> // GMPバックエンドの多倍長整数（オプション： -lgmp）
// #include <boost/multiprecision/tommath.hpp> // libtommathバックエンドの多倍長整数（オプション： -ltommath）

using namespace std;
namespace mp = boost::multiprecision;

using Bint = mp::cpp_int;  // boost独自の実装
// using Bint = mp::mpz_int;  // GMPバックエンドの多倍長整数
// using Bint = mp::tom_int;  // libtommathバックエンドの多倍長整数

const Bint DAYS = 365;

Bint Combination(const Bint &n, const Bint &k) {
    Bint res = 1;
    for (Bint i = 1; i <= k; ++i)
        res = res * (n - i + 1) / i;
    return res;
}

Bint Count(const int n) {
    Bint res = 1;
    for (Bint i = 1; i <= DAYS; ++i)
        res += ((DAYS - i) % 2 == 0 ? 1 : -1) * Combination(DAYS, i) * mp::pow(i, n);
    return res;
}

int main() {
    int n = static_cast<int>(DAYS);

    for (Bint sum = mp::pow(DAYS, n); 2 * Count(n) < sum; ++n)
        sum *= DAYS;

    cout << "n = " << n << endl;

    return 0;
}

グラフ. 確率分布

社員数 $n$ を365から4485までの20刻みとしたときの確率 $P(n)$ をプロットしたグラフを下に載せます．cpp_int型同士の割り算を実数に変換するために，分子と分母が cpp_int型の有理数cpp_rational型を使用しました．
f:id:tatanaideyo:20180807231626p:plain

C++のソースコードを表示

#include <bits/stdc++.h>
#include <boost/multiprecision/cpp_int.hpp> // boost独自の実装

using namespace std;
namespace mp = boost::multiprecision;

using Bint = mp::cpp_int;  // boost独自の実装
using Bratio = mp::cpp_rational; // 有理数型

const Bint DAYS = 365;

Bint Combination(const Bint &n, const Bint &k) {
    Bint res = 1;
    for (Bint i = 1; i <= k; ++i)
        res = res * (n - i + 1) / i;
    return res;
}

Bint Count(const int n) {
    Bint res = 1;
    for (Bint i = 1; i <= DAYS; ++i)
        res += ((DAYS - i) % 2 == 0 ? 1 : -1) * Combination(DAYS, i) * mp::pow(i, n);
    return res;
}

int main() {
    int n = static_cast<int>(DAYS);

    for (Bint sum = mp::pow(DAYS, n); n <= 4500; n += 20) {
        Bratio prob(Count(n), sum);
        cout << n << '\t' << prob.convert_to<double>() << endl;
        sum *= mp::pow(DAYS, 20);
    }

    return 0;
}