- 本記事はProjectEulerの「100番以下の問題の説明は記載可能」という規定に基づいて回答のヒントが書かれていますので、自分である程度考えてみてから読まれることをお勧めします。
問題 100. 確率2分の１になる配合
原文 Problem 100: Arranged probability
問題の要約：合計n個の青と赤のディスクが入った箱からランダムに2枚取り出したときに青が2枚になる確率$P(BB)$がちょうど50%に$n>10^{12}$で初めてなるような青のディスクの枚数を求めよ
ちょうど50%になる最初の2つの例は以下になります。
| ディスクの合計（n） | 青のディスク（b） | 青が2枚になる確率$P (BB)$ |
|---|---|---|
| 21 | 15 | $(15/21)×(14/20) \ = 1/2$ |
| 120 | 85 | $(85/120)×(84/119) \ = \frac{1}{2}$ |
100問の最後になる問題ということでかなり難しい問題となっていますね。
これは式にすると以下のようになります.
\large \frac{b(b-1)}{n(n-1)} = \frac{1}{2} \\
展開すると \\
\large 2b^2-2b-n^2+n=0
となり2元2次不定方程式のディオファントス方程式の整数解を求める問題になります。これは過去に投稿した「2元2次不定方程式をfind_DNでペル型方程式に変換して解く」の解法がそのまま使えます。
詳しいことはリンクを参照して頂けば良いのですが。要点は以下のようになります。
- 2元2次不定方程式をペル型方程式に変換(find_DN)
- ペル型方程式の解を生成(dion_DN)
- 元の解への変換Matrixを使って解を生成(transformation_to_DN)
上記のリンクで作った関数diophQuadがそのまま使うと、以下のような簡単なコードになります。
ansDQ = diophQuad("2*x**2 - 2*x - y**2 + y")
while True:
(x,y) = next(ansDQ)
print(x, y, 2*x**2 - 2*x - y**2 + y)
if y > 10**12: break
（開発環境：Google Colab）