2014年12月27日土曜日

プロジェクトオイラー 問100

プロジェクトオイラーの問題をpythonでやってみます。
日本語翻訳サイトは、プロジェクトオイラー日本語 でネット検索です。

問100「組み合わせの確率」
箱の中に15個の青い円盤と6個の赤い円盤からなる21個の色のついた円盤が入っていて, 無作為に2つ取り出すとき, 青い円盤2つを取り出す確率P(BB)は
P(BB) = (15/21) × (14/20) = 1/2
であることがわかる.

無作為に2つ取り出すとき, 青い円盤2つを取り出す確率がちょうど1/2となるような次の組み合わせは, 箱が85個の青い円盤と35個の赤い円盤からなるときである.

箱の中の円盤の合計が1012 = 1,000,000,000,000を超えるような最初の組み合わせを考える. 
箱の中の青い円盤の数を求めよ.






注意!!!
自力で解きたい人へ
以降の記述には解法に関するネタバレが含まれます。






私の回答例は以下の通りです。
def f(n):
	x, y, i = 1, 1, 0
	while i<=n:
		x, y = x+y*2, x+y
		if x*y%2: i, j = (x+1)/2, (y+1)/2
	return i,j

n = 10**12
s, b = f(n)
print "blue : ", b
print "red  : ", s-b
print "sum  :", s


問題のまま実装するととても1分では処理が終わらず、プログラミングの前に数学的に解法を整理しておかないと解けないという、プロジェクトオイラーらしい問題です。

円盤総数をi個、青円盤数をj個とすると、定義から以下が成り立ちます。
j/i * (j-1)/(i-1) = 1/2
2j(j-1) = i(i-1)
i2-i - 2(j2-j) = 0
(i - 1/2)2-(1/4) - 2(j - 1/2)2 + (1/2) = 0
(i - 1/2)2 - 2(j - 1/2)2 + 1/4 = 0
4(i- 1/2)2 - 8(j - 1/2)2 + 1 = 0
(2i-1)2 - 2(2j-1)2 = -1  ...(a)
ここで、媒介変数xとyを
x = 2i-1, y = 2j-1 ...(b)
とおくと、
x2 -  2y2 = -1 ...(c)
これは、プロジェクトオイラーではおなじみのペル方程式で、
円盤総数i=(x+1)/2が、nを超えた最初の解のときの青円盤数j=(y+1)/2が求める値です。

問94 より、ペル方程式の公式は以下の通り。
x2 - Dy2 = ±1 (Dは平方数以外。平方数ならば整数解が存在しない)
の解は、最小整数解を(x0, y0)、それ以降の解を順に、(x1, y1), (x2, y2)、...として、以下の連立方程式が成り立つ。
┌xn+1 = x0*xn + D*y0*yn
└yn+1 = y0*xn +   x0*yn ...(d)

最小整数解(x0, y0)は、式(c)から明らかで、x0=1, y0=1。また式(c)からD=2。
これらから連立方程式(d)は以下のとおり。
┌xn+1 = xn + 2yn
└yn+1 = xn + yn ...(e)
なお、式(a)より、xnとynは両方とも奇数。
これを実装します。

1.関数f(n)
・x, y, i = 1, 1, 0
 媒介変数xとyの初期値はともに1、総円盤数iの初期値は0としておきます。

・while i<=n:
 総円盤数iが指定値nを超えるまで処理を行います。

・x, y = x+y*2, x+y
 連立方程式(e)を実装します。

・if x*y%2: i, j = (x+1)/2, (y+1)/2
 媒介変数xとyの両方が奇数という判定は以下のとおりです。
 xとyの両方が奇数のときだけ、その積が奇数になるので、
 xとyの積を2で割った余りが1(論理判定ではTrue)になるかで判定します。
 そして、xとyが両方とも奇数のときに、式(b)から割り戻した値を
 円盤総数iと青円盤数jに退避しておきます。


解答はこのすぐ下の行です。文字の色を白にしてます。選択状態にすると見えます。
blue :  756872327473
red  :  313506783024
sum  : 1070379110497

0 件のコメント:

コメントを投稿