2012年12月23日日曜日

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

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


問68「"magic" 5-gon ring に当てはまる16桁の数字列は何か?」
下に示す図のようなものを"magic" 3-gon ringという.
これは1~6の数字を当てはめて, 各列の数字の和が9となっている.
これを例として説明する.

外側のノードのうち一番小さいものの付いた列(例では4,3,2)から時計回りに回ってそれぞれ列の数字を3つ連ねて説明する.
例えば例のものは4,3,2; 6,2,1; 5,1,3という組で説明することができる.
1~6の数字を当てはめて, 各列の数字の和が等しくなるものは次の8通りある.
合計
94,2,3; 5,3,1; 6,1,2
94,3,2; 6,2,1; 5,1,3
102,3,5; 4,5,1; 6,1,3
102,5,3; 6,3,1; 4,1,5
111,4,6; 3,6,2; 5,2,4
111,6,4; 5,4,2; 3,2,6
121,5,6; 2,6,4; 3,4,5
121,6,5; 3,5,4; 2,4,6
この組の各数字を連結して, 9桁の数字で表すことができる.
例えば, 上の図のものは4,3,2; 6,2,1; 5,1,3であるので432621513である.
さて, 下の図に1~10の数字を当てはめ, 各列の数字の和が等しくなる"magic" 5-gon ringを作って, それを表す16桁または17桁の数字のうち, 16桁のものの最大の数字を答えよ.
(注, 3つの場合の例を見ても分かる通り, 列の始まりの数字を比べた時一番小さい数字で始まる列から時計回りに繋げるという条件のもとで文字列を生成する必要があります. この条件下で最大となる数字を答えてください. )

-----


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






私の解答例は以下です。畳んでいます。
def g(M, Li, Lt):
	t = min(Lt)
	for L in M:
		if t in L: break
	L = [t]+sorted([s for s in L if s!=t], reverse=True)
	R = [L]
	while len(R)<len(M):
		for P in M:
			if sorted(P)!=sorted(L) and (L[-1] in P): break
		L = [s for s in P if s in Lt] + [L[-1]] \
		+ [s for s in P if s in Li and s!=L[-1]]
		R.append(L)
	return "".join([str(s) for L in R for s in L])

def f(n):
	import itertools
	N = n*2
	Q = []
	for Li in itertools.combinations(xrange(1, N+1), n):
		Li, Lt = list(Li), [s for s in xrange(1, N+1) if s not in Li]
		Sa = N*(N+1)//2 + sum([i for i in Li])
		LL = [list(L)+[Sa//n-sum(L)] for L in itertools.combinations(Li, 2)]
		for M in itertools.combinations(LL, n):
			La = [s for L in M for s in L]
			if sorted(La)==sorted(list(Li)+list(Li)+Lt): Q.append(g(M, Li, Lt))
	t = min([len(s) for s in Q])
	return max([int(s) for s in Q if len(s)==t])

n = 5
print f(n)

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

1.関数g(M, Li, Lt)
・枝リストから、一番外側の数字が最小である枝から時計回りに、
 全ての枝のノードの連結値のうちの最大値を返します。
・引数Mは枝(ノードのリスト)のリスト、Liは内側ノードのリスト、Ltは外側ノードのリストです。

・t = min(Lt)
 tは外側ノードの最小値

・for L in M:
  if t in L: break
 for文で枝リストMからノードリストをループ変数Lとして取り出します。
 if文で外側ノードの最小値がある枝の場合、for文を抜けます。
 ノードリストLは外側ノードの最小値を持つ枝で、最初の枝ということになります。

・L = [t]+sorted([s for s in L if s!=t], reverse=True)
 最初の枝のノードリストLを並び替えます。
 並び順は、外側ノード、その後に内側ノードの大きい順です。
 これで最初の枝のノード並びが確定します。
 sorted文で、内側ノードを大きい順に並び替えています。
 sorted文の第1引数のリストは内包表記しています。
 ノードリストLからループ変数としてノードsを取り出し、後ろのif文に回します。
 if文でsとtが異なる場合、つまり内側ノードの場合、sorted対象のリストの要素としてためます。
 そして第2引数reverse=Trueで、逆順(大きい順)に並べるように指定します。
 外側ノードtを[]でリスト化し、内側ノードリストと結合して、改めてリストLとします。

・R = [L]
 ノード並びが確定したので、最初の枝を結果リストRにためます。

・while len(R)<len(M):
 結果リストRの要素数が元の枝の数になるまで、2本目以降の枝の処理をします。

・for P in M:
  if sorted(P)!=sorted(L) and (L[-1] in P): break
 前の枝の末尾ノードは次の枝にだけ共有されていることを利用して、次の枝Pを求めます。
 「for P in M」で枝リストMから枝をループ変数Pとします。
 if文で枝P枝と1つ前の枝Lとをそれぞれソートして異なること、そして前の枝Lの末尾ノードL[-1]が存在するような枝Pになったところで処理を抜けます。

・L = [s for s in P if s in Lt] + [L[-1]] \
  + [s for s in P if s in Li and s!=L[-1]]
 枝のノードを外側からの順に並び替えます。
 並び替え後ノードリストLは、Pの次のループで前の枝として扱います。

 まず、[s for s in P if s in Lt]は外側ノード1つだけのリストです。
 内包表記です。まず枝Pからノードを1つずつ取り出しループ変数sとして
 後ろのif文に回します。そして外側ノードリストLtにある場合だけ先頭に回します。

 次に、[L[-1]]は1つ前の枝の最後の値だけのリストです。枝Pでは2番目の要素です。

 最後に、[s for s in P if s in Li and s!=L[-1]]は上記2つ以外のノードです。
 内包表記です。まず枝Pからノードを1つずつ取り出しループ変数sとして
 後ろのif文に回します。そして内側ノードリストLiにあり、先ほどの2番目の要素と異なる場合だけ採用し先頭に回します。

 この文全体で上記3要素それぞれだけのリスト3つを+演算子で1つのリストに結合します。

・R.append(L)
 枝のノードが実際の並びになったので、結果リストRに追加してためます。

・return "".join([str(s) for L in R for s in L])
 すべての枝のノードを並び直したので、全ノードを文字扱いして順に連結た値を返します。
 内包表記です。前のfor文で結果リストRから枝Lを取り出し後ろに回します。
 後ろのfor文で枝Lからノードsを取り出し先頭に回します。
 先頭で文字型に変換し、その全要素をリストにします。
 そしてこのリストをjoin関数で区切り文字なしで連結します。

2.関数f(n):
・マジックn角形リングについて、問題で求めるノードの連結値を返します。
 具体的には、マジックn角形リングの外側から内側へ各ノードの値を連結した枝のうち、一番外側の数字が最小の枝から時計回りに、全ての枝数字を連結して最小桁数の最大値を返します。
 引数nは枝の本数でもあります。

・import itertools
 組合せの関数を使用するために、標準モジュールの「itertools」を使えるようにしておきます。

・N = n*2
 Nはノード数です。マジック五角形はノードが10個あります。

・Q = []
 Qは返却する連結値候補リストです。初期クリアしておきます。

・for Li in combinations(xrange(1, N+1), n):
 Liはすべての内側ノードのリストです。
 xrange関数で1からNまでの整数を発生させ、combinations関数でそのうちn個ずつの組合せをLiとします。この時点ではLiはタプルです。

・Li, Lt = list(Li), [s for s in xrange(1, N+1) if s not in Li]
 Liは内側ノードリストで、タプルからリストに変換しておきます。
 Ltは外側ノードリストです。1からNまで整数のうち内側ノードリストに無い値のリストです。
 内包表記です。1からNまで整数をループ変数sとして後ろに回し、if文で内側ノードリストLiに無い値だけ、先頭に回してリストにためます。

・Sa = N*(N+1)//2 + sum([i for i in Li])
 Saはライン総和です。ライン総和には内側ノード2回分と外側ノード1回分が使われるので、「ノード総和+内側ノード総和」としています。
 「N*(N+1)//2」はノード総和です。1からNまでの整数の和の公式そのままです。
 //演算子は割り算で商の整数部を返します。
 「sum([i for i in Li])」は内側ノードLiの総和です。

・LL = [list(L)+[Sa//n-sum(L)] for L in combinations(Li, 2)]
 リストLLは枝候補リストです。内包表記で記述しています。
 まず内側ノードLiからノード2つの組合せタプルをループ変数Lとして先頭に回します。
 list(L)は直後にリストを連結するため、編集不可のタプルから編集可能なリストへ変換しておきます。
 「Sa//n」はライン総和Saを枝数nで割ることで枝ライン1つ分の和のことです。
 [Sa//n-sum(L)]は、この枝1つ分の和から内側ノード2つ分の和sum(L)を差し引くことで、 この枝の外側ノードの値となります。
 list(L)+[Sa//n-sum(L)の全体で、枝候補リストになります。

・for M in combinations(LL, n):
 枝候補リストLLから、枝本数分の組合せを取り出した枝候補組合せリストをループ変数Mとします。

・La = [s for L in M for s in L]
 Laは枝候補組合せリストMの全ノードリストです。内包表記しています。
 枝候補組合せリストMから枝Lを取り出し後ろに回します。
 そして取り出した枝Lからノードsを取り出し、先頭に回しそのままリストにためます。

・if sorted(La)==sorted(Li+Li+Lt): Q.append(g(M, Li, Lt))
 上記の枝候補組合せリストの全ノードについて、内側ノードLiが2回と外側ノードLtが1回使われているかチェックします。
 リストLaと、内側ノードLi2つ分と外側ノードLtの連結リストについて、
 それぞれsorted関数で小さい順に並べて、まったく同じになるかをチェックします。
 チェックOKならば、関数gに枝候補組合せリストM、内側ノードリストLi、外側ノードリストLtを
 渡してます。
 そして、関数gの戻り値である、外側ノード最小値からのノード連結値が最大になるようにを並べ替えた
 連結値を連結値候補リストQにためます。

ここまでで全ての枝候補組合せについて、その連結値候補を取得しました。

・t = min([len(s) for s in Q])
 tは連結候補の最小桁数です。
 連結値候補リストQから連結値をループ変数sとして取り出し先頭に回し、
 len関数で桁数を求めてリストにためます。
 min関数でそのリストの最小値を取得します。

・return max([int(s) for s in Q if len(s)==t])
 連結値候補リストQから連結値をループ変数sとして取り出し、後ろのif文に回します。
 そして最小桁数と一致する場合だけ先頭に回し、int関数で数値化してリストにためます。
 max関数でそのリストの最大値を取得し、呼び出し元に返却します。

上記処理は1秒程度で終わるので1分ルールは守っていますが、
処理速度をさらに向上するアイデアは以下の通りです。
・ライン総和Saがライン数nで割り切れない組は処理を飛ばす。
・枝候補リストLLでの外側ノード値が、外側ノードリストLtにない組は処理を飛ばす。
・枝候補リストLLでの外側ノード値で、外側ノードリストLtの全要素を使ってなければ処理を飛ばす。

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

(追記)
・20130120 ネタバレ注意の追記など

2012年12月10日月曜日

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

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



問67「効率的なアルゴリズムを使い三角形内の最大の和を求めよ」

以下の三角形の頂点から下まで移動するとき, その数値の合計の最大値は23になる.
3
7 4
4 6
8 5 9 3
この例では 3 + 7 + 4 + 9 = 23
100列の三角形を含んでいる15Kのテキストファイル triangle.txt (右クリックして, 『名前をつけてリンク先を保存』)の上から下まで最大合計を見つけてください.
:これは, Problem 18のずっと難しいバージョンです.
全部で299 通りの組み合わせがあるので, この問題を解決するためにすべてのルートをためすことは可能でありません!
あなたが毎秒1兆本の(1012)ルートをチェックすることができたとしても, 全てをチェックするために200億年以上かかるでしょう. 
解決するための効率的なアルゴリズムがあります. ;o)

-----


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





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

def f(u):
	#1.init
	p = []
	for i in u.split("\n"):
		if not i: continue
		p.append([int(j) for j in i.split() if j])
	n = len(p)
	A = [["" for i in xrange(n)]]
	q = p[n-1][:]

	#2.val List, LR List
	for i in xrange(n-2, -1, -1):
		L = []
		for j in xrange(i+1):
			s, t = q[j], q[j+1]
			v = max(s, t)
			L.append({s:"L", t:"R"}[v])
			q[j] = p[i][j] + v
		A.insert(0, L)
		q.pop(i+1)
 
	#3.Route Search for MaxValue
	L, M, i = [], [], 0
	#for s, t in zip(A[::-1], p):
	for s, t in zip(A, p):
		L.append(s[i])
		M.append(t[i])
		if s[i]=="R": i+=1
	return sum(M), M, L

import os, sys
sIn = os.path.join(os.path.dirname(sys.argv[0]), "problem067_triangle.txt")
u = open(sIn).read()
s, val, LR = f(u)

print "sum:", s
print "val:", val
print "LR :", LR

1個の関数を使っています。
問18の関数f(u)そのままです。
関数の外側で、提示されたファイルを読み込んで、関数f(u)を呼び出しています。
sys.argv[0]は、当pythonファイルのフルパスです。
os.path.dirname(sys.argv[0])で、当pythonファイルのフォルダのフルパスを求め、os.path.joinで、これと、当pythonファイルと同じフォルダ位置にある当問題の提示ファイル名と連結しています。

解答はこのすぐ下の行です。文字の色を白にしてます。選択状態にすると見えます。
sum: 7273
val: [59, 73, 52, 53, 87, 57, 92, 81, 81, 79, 81, 32, 86, 82, 97, 55, 97, 36, 62, 65, 90, 93, 95, 54, 71, 77, 68, 71, 94, 8, 89, 54, 42, 90, 84, 91, 31, 71, 93, 94, 53, 69, 73, 99, 89, 47, 80, 96, 81, 52, 98, 38, 91, 78, 90, 70, 61, 17, 11, 75, 74, 55, 81, 87, 89, 99, 73, 88, 95, 68, 37, 87, 73, 77, 60, 82, 87, 64, 96, 65, 47, 94, 85, 51, 87, 65, 65, 66, 91, 83, 72, 24, 98, 89, 53, 82, 57, 99, 98, 95]
LR : ['L', 'L', 'R', 'R', 'R', 'R', 'L', 'R', 'L', 'R', 'L', 'R', 'R', 'R', 'R', 'R', 'R', 'L', 'L', 'R', 'L', 'L', 'R', 'L', 'R', 'L', 'R', 'R', 'L', 'L', 'R', 'R', 'R', 'R', 'R', 'R', 'R', 'R', 'L', 'L', 'R', 'R', 'L', 'R', 'R', 'R', 'R', 'R', 'L', 'L', 'L', 'R', 'L', 'R', 'R', 'R', 'L', 'L', 'L', 'L', 'L', 'L', 'R', 'R', 'R', 'R', 'R', 'L', 'R', 'L', 'L', 'L', 'L', 'L', 'L', 'R', 'L', 'L', 'R', 'R', 'L', 'L', 'L', 'L', 'L', 'R', 'L', 'L', 'L', 'R', 'L', 'R', 'R', 'L', 'R', 'R', 'R', 'L', 'R', '']
(追記)
・20130127 ネタバレ注意の追記など

2012年12月6日木曜日

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

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

問66「ディオファントス方程式 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のときにピタリと処理が進まなくなり、
1分ルールどころではありませんでした。

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

まず、Dの1つひとつについて問題の式が成り立つ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の近似分数の分子と分母です。
連分数の近似分数は問65の問題文にあるように
> 平方根の部分的な連分数の数列から良い有理近似が得られることが分かる。
とのことです。

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

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

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

(参考文献2による式) ※文字式の文字は当ページに合わせて変えています。
┌      ┐┌      ┐      ┌      ┐ ┌           ┐
│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行列の列成分の左右が異なりますが、本質的には同じです。私は参考文献1の式を実装しました。

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

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

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

1.関数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)
 行数、列数を返します。

2.関数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のループで、上記の1行分のリストを作成します。
 さらにiのループで、上記の全ての行をリストにして、結果の行列として返却します。

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

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

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

・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]]
 参考文献1に基づいた、行列の初期値です。
 sは√Dの連分数表示での固定部分です。

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

・Q = [[0, 1], [1, i]]
 参考文献1に基づいた、ループ変数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とします。

5.関数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以下の整数を1つずつ設定します。

・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を呼び出し元に返します。

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

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

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

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


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