まえがき
いきなりですが、複数人でランダムにプレゼント交換を行うことを考えましょう。
プレゼント交換を行うというのですから、自分のプレゼントが自分に帰ってきてしまっては嫌ですよね。プレゼント交換に参加した人全員についてそのような状況が起こらない確率を$p$としましょう。
実は、プレゼント交換に参加する人数を増やしていくと、$p\rightarrow\frac{1}{e} (eはネイピア数)$に収束することが知られています。
これは攪乱順列(完全順列)と呼ばれる問題で、漸化式または包除原理を用いると$p$の一般項に
$e$のマクローリン展開の形が現れることから分かります。詳しい説明や導出が気になる方は以下の記事を参考にしてください。
受験の月
高校数学の美しい物語
ユニマス
実際にシミュレーションしてみる
実際にプログラムを動かしてプレゼント交換のシミュレーションを行うことによって$p$を求めます。人数が十分多い場合には$p\approx\frac{1}{e}$だったので、$\frac{1}{p}\approx e$となることからネイピア数$e$(の近似値)を求めてみましょう。
C++のメルセンヌ・ツイスタという疑似乱数生成アルゴリズムを用いて実装しました。なお、メルセンヌ・ツイスタを用いる理由としては高品質かつ高速な疑似乱数生成が可能だからです。
プレゼント交換のシミュレーション
#include <bits/stdc++.h>
#include <random>
using namespace std;
int main() {
random_device seed_gen;
mt19937 engine(seed_gen());
int n=10;
double cnt=0, success=0;
while(cnt<300000000) {
vector<int> presents(n); //人iが今持ってるプレゼント
iota(presents.begin(),presents.end(),0); //最初人iはプレゼントiを持っている
for(int i=n-1; i>=0; i--) {
int r=engine()%(i+1); //人iはプレゼント交換する相手を人0~iから選ぶ
swap(presents[i],presents[r]); //人iと人rがプレゼント交換
if (presents[i]==i) break; //人iにプレゼントiが帰ってきた
if (i==0) success++; //誰も自分のプレゼントが自分に帰ってこなかった
}
cnt++;
}
double p=success/cnt;
cout << fixed << setprecision(7) << "1/e≈" << p << endl;
cout << fixed << setprecision(7) << "e≈" << 1./p << endl;
}
このコードを実行した結果は以下のようになりました。
ネイピア数$e=2.71828 18…$と比較すると、少数第3位まで一致しました。
補足
プレゼント交換を行う人数が$n=10$であることについて軽く説明します。まず、$e$のマクローリン展開は以下のようになります。ちなみにマクローリン展開とは、一般の関数を多項式で近似する手法のことです。
e^x=\sum_{k=0}^\infty\frac{x^k}{k!}
上式に$x=1$を代入すれば$e$の値が求まります。ここで、$x=1$としたとき各$k$に対する分母は$k!$ですから、早めに切り上げても問題なさそうです。実際に、$k=10$で足し合わせを打ち切っても少数第7位まで一致する精度でネイピア数の近似値が求まります。
以上より、プレゼント交換に参加する人数は10人で抑え、試行回数を増やしました。
順列全探索で求める
今回はシミュレーションしましたが、$n=10$程度なら順列全探索の計算量が$O(n!)$なので十分高速に動作します。(C++では大体$O(10^7)$から$O(10^8)$くらいまでなら$2s$以内に実行が終了し、いま$10!=3628800=O(10^6)$です。)
つまり、全パターンを探索して理論的な確率を求めればより正確に$e$の近似値を求めることが出来ます。(もちろんシミュレーションをするより正確に...)
以下にそのコードを示します。
プレゼント交換のシミュレーション
#include <bits/stdc++.h>
using namespace std;
int main()
{
int n=10;
vector<int> presents(n);
iota(presents.begin(), presents.end(), 0);
double count=0, success=0;
do{
count++;
for(int i=0; i<n; i++){
if(presents[i]==i) break;
if (i==n-1) success++;
}
}while(next_permutation(presents.begin(), presents.end()));
double p=success/count;
cout << fixed << setprecision(10) << "1/e≈" << p << endl;
cout << fixed << setprecision(10) << "e≈" << 1./p << endl;
}
少数第6位まで$e$の理論値と一致しました。
あとがき
今回はプレゼント交換を行うことによってネイピア数を求めました。確率とかに$e$や$\pi$が出てきて、シミュレーションで近似値を求める系が好きなのでやってみました。
時間がある方は実際にプレゼント交換を繰り返してネイピア数を求めてみてください。