漫ろ草 SozoroGusa: プロジェクトオイラー　問６６

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

問６６「ディオファントス方程式 x2 - Dy2 = 1 を調べ上げよ」

次の形式の, 2次のディオファントス方程式を考えよう:

x2 - Dy2 = 1

たとえば D=13 のとき, x を最小にする解は 6492 - 13×1802 = 1 である.

D が平方数(square)のとき, 正整数のなかに解は存在しないと考えられる.

D = {2, 3, 5, 6, 7} に対して x を最小にする解は次のようになる:

32 - 2×22 = 1

22 - 3×12 = 1

92 - 5×42 = 1

52 - 6×22 = 1

82 - 7×32 = 1

したがって, D ≤ 7 に対して x を最小にする解を考えると, D=5 のとき x は最大である.

D ≤ 1000 に対する x を最小にする解で, x が最大になるような D の値を見つけよ.

－－－－－

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

私の解答例は以下です。畳んでいます。

def mc(L):
	M = [len(s) for s in L ]
	if min(M)-max(M): return False
	return len(L), len(s)

def mm(L, M):
	(Lr, Lc), (Mr, Mc) = mc(L), mc(M)
	if Lc!=Mr: return False
	return [[sum([L[i][k]*M[k][j] for k in xrange(Lc)]) \
			for j in xrange(Mc)] for i in xrange(Lr)]

def g(n, b, c, L, M):
	if (b, c) in M: return L
	if n==b*b: return L
	M.append((b, c))
	C, t = (n-b*b)/c, int(n**(1.0/2))
	u = (b+t)%C
	L.append((b+t)//C)
	return g(n, t-u, C, L, M)

def h(D):
	import math
	s = int(math.sqrt(D))
	L = g(D, s, 1, [], [])
	if not L: return [False, False]
	P = [[0, 1], [1, s]]
	for i in L+L:
		Q = [[0, 1], [1, i]]
		R = mm(P, Q)
		x, y = R[1][1], R[0][1]
		if x*x-D*y*y==1: return [x, y]
		P = R[:]

def f(n):
	r = [0, 0, 0]
	for D in xrange(n+1):
		x, y = h(D)
		if x and r[1]<x: r = [D, x, y]
	return r

n = 1000
print f(n)

最初に考えた方法は、「x2 - Dy2 = 1」が成り立つようにyをループさせながら
チェックするという以下の方法です。

def g(D):
	for y in xrange(2, 1000000000):
		x = math.sqrt(D*y*y+1)
		if round(x)==x: return int(x), y
	return [0, 0]

ですが、この方法ではD=61のときにピタリと処理が進まなくなり、
１分ルールどころではありませんでした。

まったく解法の糸口が見えず、ネット検索したところ、問６４,問６５でやってきた連分数を使えばできそうなことがわかり、本などを参考にしました。
末尾に参考文献として記載しました。

まず、Dの１つひとつについて問題の式が成り立つx,yを求めます。
問題の式をDで解くと、

x2 - Dy2 = 1
D*y2 = x2 - 1
x2 - 1 x 2 1
D = -------- = (---) - ----
y2 y y2

x,yが大きな値になる場合、1/y2は極めて小さい値になり、
x 2
D → (---)
y
√D → ±(x/y)
となり、√Dが±(x/y)に近づいていきます。
xとyは√Dの近似分数の分子と分母です。
連分数の近似分数は問６５の問題文にあるように
> 平方根の部分的な連分数の数列から良い有理近似が得られることが分かる。
とのことです。

そこで、連分数によって近似分数を求め、その分子分母を候補として
問題の式にあてはめて、これが成り立つところまで近似精度を上げていきます。

連分数の入れ子具合の階層の深さ、近似分数の計算の深さをnとして、
近似分数の分子をx(n)、分母をy(n)とします。
このとき以下の式が成り立ちます。
（証明は数学的帰納法でできますがここでは省略します。）

x(n) 1
---- = a(0) + ------------------
y(n) 1
a(1) + -------------
・・・・ 1
a(n-1) + ----
a(n)
において、
（参考文献１による式）　※文字式の文字は当ページに合わせて変えています。
┌ ┐┌ ┐ ┌ ┐ ┌ ┐
│0 1 ││0 1 │ │0 1 │ │y(n-1) y(n)│
│ ││ │・・・│ │=│ │
│1 a(0)││1 a(1)│ │1 a(n)│ │x(n-1) x(n)│
└ ┘└ ┘ └ ┘ └ ┘

（参考文献２による式）　※文字式の文字は当ページに合わせて変えています。
┌ ┐┌ ┐ ┌ ┐ ┌ ┐
│a(0) 1││a(1) 1│ │a(n) 1│ │x(n) x(n-1)│
│ ││ │・・・ │ │=│ │
│ 1 0││ 1 0│ │ 1 0│ │y(n) x(n-1)│
└ ┘└ ┘ └ ┘ └ ┘

x,yが分子分母のどちらかということと、結果の2×2行列の列成分の左右が異なりますが、本質的には同じです。私は参考文献１の式を実装しました。

なお、連分数の循環部分は１回では足りず２回繰り返す必要があります。
近似分数は真の値よりやや大きい、やや小さいを繰り返しながら真の値に近づきますので、
当問題の提示式、
x2 - D*y2 = 1　に対して、
x2 - D*y2 = -1　が成り立つ場合があり得ます。
この場合、さらにもう１回分循環部分を繰り返すなかで
当問題の提示式が成り立つ組が見つかります。

なお、行列計算がでてきますが、これにはpython数値計算ライブラリnumpyを
別途、ダウンロード、インストールして使用しても可能です。

５つの関数を使用しています。

１．関数mc(L)
・行列のリストから行数と列数を返します。
　行列は[[a, b],[c, d]]のように行ごとのリストが列の数だけあるものが対象です。

・M = [len(s) for s in L]
　リストMは行列Lの行要素sの要素数です。つまり、各行ごとの列数のことです。

・if min(M)-max(M): return False
　各行の列数がすべての行で同じであるかチェックします。
　具体的にはリストMの要素がすべて同じこと、最大値と最小値が一致することを
　チェックし、不一致ならばFalseを返します。

・return len(L), len(s)
　行数、列数を返します。

２．関数mm(L, M)
・行列リストLとMの積を返します。
　行列は[[a, b],[c, d]]のように行ごとのリストが列の数だけあるものが対象です。

・(Lr, Lc), (Mr, Mc) = mc(L), mc(M)
　行列LとMのそれぞれの行数、列数を取得します。rが行(row),cが列(column)です。

・if Lc!=Mr: return False
　行列Lの列数と行列Mの行数が一致しない場合、
　行列の積は計算できないのでFalseを返します。

・return [[sum([L[i][k]*M[k][j] for k in xrange(Lc)]) \
for j in xrange(Mc)] for i in xrange(Lr)]
　ループ変数iは結果の行列の行数、jは結果の行列の列数、
　kは結果の行列の各行の中での列番号です。

　まずkのループがある、[L[i][k]*M[k][j] for k in xrange(Lc)]は
　行列Lと行列Mの要素から各要素の積の配列を求めます。
　sum関数でこの和を求め、結果の行列の各要素とします。
　行列計算の要素の定義そのものです。

　次にjのループで、上記の１行分のリストを作成します。
　さらにiのループで、上記の全ての行をリストにして、結果の行列として返却します。

（参考）行列(wikipedia)　
　　上記ページの中ほどの「行列の積」。添え字のi,j,kも当関数と一緒です。

３．関数g(n, b, c, L, M)
・連分数の循環部分リストを取得します。
　問６４の関数gと同じです。
　例）√23 = [4;(1,3,1,8)]について、[1,3,1,8]を返します。

４．h(D)
・「x2 - Dy2 = 1」を満たす最小のx,yを返します。Dは平方根を求めたい整数です。
・参考文献１による式を実装しました。

・import math
　s = int(math.sqrt(D))
　sは√Dの連分数表示の固定部分、つまり√Dの値の整数部分です。
　モジュールmathをインポートして求めます。

・L = g(D, s, 1, [], [])
　Lは√Dの連分数表示の循環部分です。関数gで求めます。

・if not L: return [False, False]
　連分数の循環部分が無い場合、つまり関数gで平方根を求めるために設定した引数Dが
　小数部分を持たない平方数であった場合は、以降の計算ができません。
　問題文の
「Dが平方数(square)のとき, 正整数のなかに解は存在しないと考えられる.」
　に相当するケースです。
　この場合は、呼び出し元に分子分母ともFalseとして返します。

・P = [[0, 1], [1, s]]
　参考文献１に基づいた、行列の初期値です。
　sは√Dの連分数表示での固定部分です。

・for i in L+L:
　ループ変数iには、√Dの連分数表示での循環部分の値を１つずつ設定します。

・Q = [[0, 1], [1, i]]
　参考文献１に基づいた、ループ変数iを使った次の順番の行列です。

・R = mm(P, Q)
　関数mmで行列Pと行列Qの行列積を求め、行列Rに設定します。

・x, y = R[1][1], R[0][1]
　行列の積の結果行列Rから分子分母に相当する要素を取り出しそれぞれx,yに設定します。

・if x*x-D*y*y==1: return [x, y]
　問題の提示式が成り立つときは、そのときのx,yを呼び出し元に返却します。

・P = R[:]
　結果行列をループの次の行列積の準備として、行列Pとします。

５．関数f(n)
・「x2 - Dy2 = 1」において、n以下の正の整数Dに対するxを最小にする解で、
　xが最大になるようなDの値、及びそのときのx、yを返します。

・r = [0, 0, 0]
　リストrは、返却用リストです。初期クリアしておきます。

・for D in xrange(n+1):
　ループ変数Dに0以上n以下の整数を１つずつ設定します。

・x, y = h(D)
　関数hで、Dの場合の「x2 - Dy2 = 1」を満たす最小の値を求め、x, yとします。

・if x and r[1]<x: r = [D, x, y]
　xに値があり、かつそれが返却用リストのx値よりも大きい場合、
　現在のD,x,yを返却用リストrに保持します。

・return r
　ループ終了後、返却用リストrを呼び出し元に返します。

６．関数の外
・n = 1000
　print f(n)
　問題に合うように1000を関数fに渡し戻り値を表示します。

（参考文献１）
木村俊一,連分数のふしぎ,講談社,ブルーバックス1770,2012/05/20
「６　行列と連分数」　p304、他。

（参考文献２）
佐藤篤,連分数論入門（東北大学理学部数学科のサイト、pdfファイル）
「７　Pell方程式」　p49、他。

解答はこのすぐ下の行です。文字の色を白にしてます。選択状態にすると見えます。
661,
16421658242965910275055840472270471049L,
638728478116949861246791167518480580L

（追記）
・20121208　参考文献２のURLを直接表示から文献名のリンクに修正。
・20130127　ネタバレ注意の追記など

漫ろ草 SozoroGusa

2012年12月6日木曜日

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

0 件のコメント:

コメントを投稿

リンク

2012年12月6日木曜日

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

0 件のコメント:

コメントを投稿

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