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 ネタバレ注意の追記など

2012年11月19日月曜日

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

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


問65「e についての連分数である近似分数の100項目の分子の桁の合計を求めよ」
2の平方根は無限連分数として書くことができる。
            1
√2 = 1 + -------------------
                1
          2 + ---------------
                    1
              2 + -----------
                        1
                  2 + -------
                      2 + ...

無限連分数である√2 = [1;(2)]と書くことができるが、
(2)は2が無限に繰り返されることを示す。

同様に、√23 = [4;(1,3,1,8)]。
平方根の部分的な連分数の数列から良い有理近似が得られることが分かる。

√2の近似分数について考えよう。
    1   2
1 + - = -
    2   3

      1     7
1 + ----- = -
        1   5
    2 + -
        2

      1         17
1 + --------- = --
          1     12
    2 + -----
            1
        2 + -
            2

      1             41
1 + ------------- = --
          1         29
    2 + ---------
             1
        2 + -----
                1
            2 + -
                2

従って、√2の近似分数からなる数列の最初の10項は:
1, 3/2, 7/5, 17/12, 41/29, 99/70, 239/169, 577/408, 1393/985, 3363/2378, ...


もっとも驚くべきことに、数学的に重要な定数、
e = [2; 1,2,1, 1,4,1, 1,6,1 , ... , 1,2k,1, ...]。

eの近似分数からなる数列の最初の10項は:
2, 3, 8/3, 11/4, 19/7, 87/32, 106/39, 193/71, 1264/465, 1457/536, ...

10項目の近似分数の分子の桁を合計すると1+4+5+7=17である。
eについての連分数である近似分数の100項目の分子の桁の合計を求めよ。
-----



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





私の解答例は以下です。畳んでいます。
def gcd(a, b):
	if a<b: a, b = b, a
	if b: return gcd(b, a%b)
	else: return a

def h(n):
	if n<0: return 0
	L = [2]
	for i in xrange(n-1):
		L.append({True:(i//3+1)*2, False:1}[i%3==1])
	return L

def g(L):
	b, c = L.pop(), 1
	for a in L[::-1]:
		x, y = a*b+c, b
		s = gcd(x, y)
		b, c = x//s, y//s
	return b, c

def f(n):
	L = h(n)
	x, y = g(L)
	s = str(x)
	M = [int(s[i]) for i in xrange(len(s))]
	return sum(M), x, y

n = 100
print f(n)


4つの関数を使っています。

1.関数gcd(a, b)
・「ユークリッドの互除法」を使用して、引数a,bの最大公約数を返します。
 説明は問5を参照してください。


2.関数h(n)
・n項までのネイピア数の連分数リストを返します。
 但し問題文ではセミコロンで示している区切り文字も全てカンマ区切りにします。
 e = [2; 1,2,1, 1,4,1, 1,6,1 , ... , 1,2k,1, ...]
  →[2, 1,2,1, 1,4,1, 1,6,1 , ... , 1,2k,1, ...]


・if n<0: return 0
 nが0未満は未定義なので0を返します。


・L = [2]
 初期値は連分数の固定部分なので2です。


・for i in xrange(n-1):
 ループ変数iには0以上n-1未満のn個の整数が順に設定されます。


・L.append({True:(i//3+1)*2, False:1}[i%3==1])
 「1,2k,1, ...」を発生させ、リストLに設定します。
 %演算子は割り算の余りです。
 [i%3==1]で3つごとに余り1になればTrue、他はFalseです。
 これを前の{}内の辞書要素のうち当てはまる値を設定します。
 //演算子は割り算の商です。True、つまり、3で割り余り1になる場合、
 「(i//3+1)*2」で上記式の2kにあたる部分の値を設定し、その他は1を設定します。


3.関数g(L):
・リストLとして連分数の固定部分aと循環部分b1,b2,...を受け取り、
 その近似分数x/yの分子x、分母yを返します。式で書くと以下の通りです。
      1             x
a + ------------- = -
           1        y
    b1 + --------
               1
         b2 + ---
              ...


・b, c = L.pop(), 1
 リストのpop関数はリストの要素を取り出して、それをリストから削除します。
 引数未指定の場合は最後の要素を取り出して、それをリストから削除します。
 初期値として、bに連分数の循環部分の最後、cに1を設定します。


・for a in L[::-1]:
 リストLを逆順に1つずつ取り出しループ変数aに設定します。


・x, y = a*b+c, b
    c   a*b+c   x
a + - = ----- = -
    b     b     y
上記式を実装しました。


・s = gcd(x, y)
 計算したx,yが約分できることがあるので、関数gcdで最大公約数を求めておきます。


・b, c = x//s, y//s
 先ほど求めたx,yをその最大公約数で割り、できる限り簡単な分数に直します。
 これを次のb,cとしてループ変数aのループを最後まで回します。


・return b, c
 連分数リストLから算出したその近似分数の分子と分母を返します。

4.関数f(n):
・ネイピア数のn項目までの連分数の近似分数の分子の各桁の合計、及びその分子、分母を返します。


・L = h(n)
 関数h(n)でn項までのネイピア数の連分数リストを取得しLとします。


・x, y = g(L)
 関数g(L)で上記リストLの近似分数の分子と分母を取得し、x,yとします。


・s = str(x)
 L = [int(s[i]) for i in xrange(len(s))]
 上記で求めた分子xを文字型にして、1桁ずつ数値にしたリストをMとします。


・return sum(M), x, y
 分子の各桁の合計値、分子、分母を返します。


5.関数の外
・n = 100
 print f(n)
 問題に合うように100を関数fに渡し戻り値を表示します。
 
解答はこのすぐ下の行です。文字の色を白にしてます。選択状態にすると見えます。
272,
6963524437876961749120273824619538346438023188214475670667,
2561737478789858711161539537921323010415623148113041714756

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

2012年11月11日日曜日

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

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

問64「N ≤ 10000 について奇数の周期をもつ連分数はいくつあるか?」
平方根は連分数の形で表したときに周期的であり, 以下の形で書ける:
           1
N = a0 + --------------------
                1
         a1 + ---------------
                     1
              a2 + ----------
                          1
                   a3 + -----
                         ...

例えば, √23を考えよう.
                             1              1
√23 = 4 + 23 - 4 = 4 + ------------ = 4 + ---------------
                             1                 √23 - 3
                          --------         1 + --------
                          √23 - 4                  7
となる.
この操作を続けていくと,
              1
√23 =  4 + --------------------------
                   1
             1 + ---------------------
                        1
                  3 + ----------------
                             1
                       1 + -----------
                                   1
                             8 + -----
                                  ...
を得る.
操作を纏めると以下になる:
a0 = 4, 1/(√23-4) =  (√23+4)/7  = 1 + (√23-3)/7
a1 = 1, 7/(√23-3) = 7(√23+3)/14 = 3 + (√23-3)/2
a2 = 3, 2/(√23-3) = 2(√23+3)/14 = 1 + (√23-4)/7
a3 = 1, 7/(√23-4) = 7(√23+4)/7  = 8 + (√23-4)
a4 = 8, 1/(√23-4) =  (√23+4)/7  = 1 + (√23-3)/7
a5 = 1, 7/(√23-3) = 7(√23+3)/14 = 3 + (√23-3)/2
a6 = 3, 2/(√23-3) = 2(√23+3)/14 = 1 + (√23-4)/7
a7 = 1, 7/(√23-4) = 7(√23+4)/7  = 8 + (√23-4)

よって, この操作は繰り返しになることが分かる.
表記を簡潔にするために,
√23 = [4;(1,3,1,8)] と表す.
(1,3,1,8)のブロックは無限に繰り返される項を表している.

最初の10個の無理数である平方根を連分数で表すと以下になる.
√2=[1;(2)], period=1
√3=[1;(1,2)], period=2
√5=[2;(4)], period=1
√6=[2;(2,4)], period=2
√7=[2;(1,1,1,4)], period=4
√8=[2;(1,4)], period=2
√10=[3;(6)], period=1
√11=[3;(3,6)], period=2
√12= [3;(2,6)], period=2
√13=[3;(1,1,1,1,6)], period=5
N ≦ 13で奇数の周期をもつ平方根は丁度4つある.

N ≦ 10000 について奇数の周期をもつ平方根が何個あるか答えよ.
-----



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






私の解答例は以下です。畳んでいます。
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 f(n):
	t = 0
	for i in xrange(n+1):
		s = g(i, int(i**(1.0/2)), 1, [], [])
		if len(s)%2: t += 1
	return t

n = 10000
print f(n)


以下のように考えました。
  c                     √n-B
------ ・・・(A),   A + ----- ・・・(B)
√n-b                     C
として、(A)を連分数(B)にすることを考えます。
なお、(B)で計算したAの値を循環部分として保存し、
B, Cの値を次世代のb, cとして再帰的に計算します。
まず(A)を有理化します。
  c          c*(√n+b)     c*(√n+b)
------ = --------------- = ----------- ・・・(A)'
√n-b     (√n-b)*(√n+b)    (n-b*b)
連分数展開後の(B)の形に合わせるようにて分母Cを決めます。
    (n-b*b)
C = ------- ・・・(C)
       c
(A)'に(C)を当てはめ、
(√n+b)
------- ・・・(D)
   C
(D)から、連分数(B)のAに相当する部分をくくりだします。
そのために、√nの整数部分tとして、
t = int(√n) ・・・(E)
(b+t)   √n-t
----- + ------ ・・・(F)
  C       C
ですが、連分数(B)のAは整数のため、
Aは(F)の第1項の整数部分で、余りuは第2項に含めます。
A = (b+t)//C ・・・(G)
u = (b+t)%C  ・・・(H)
(G)(H)を(F)に反映して
 (b+t)        √n-t +u
{----- - u} + ---------
   C               C
             √n - (t-u)
= (b+t)//c + -----------
                 C
以上より、A,B,Cをまとめておくと以下の通りです。
(E)(H)のt,uを使って、
A = (b+t)//C ・・・(G)
B = t-u      ・・・(I)
C = (n-b*b)/c・・・(C)

2つの関数を使っています。

1.関数g(n, b, c, L, M)
  c                        √n-B
------ を連分数展開して、A + ------ (分数部分の絶対値は1未満)
√n-b                         C
に式変形したときの、循環部分のAのリストLを取得します。
リストMにはbとcの組、タプル(b,c)を貯めます。
例)√23 = [4;(1,3,1,8)]の場合、[1,3,1,8]を返します。
・関数内に自分自身を呼び出す記述のある再帰関数です。

・if (b, c) in M: return L
 まず再帰の終了条件です。
 引数b,cの組が今までにでてきていれば、循環部分がひと回りしたことになるので処理を終了し循環部分のAのリストLを呼び出し元に戻します。

・if n==b*b: return L
 この条件の場合、入力値の分母が0でゼロ割り算エラーのため連分数にならないのでこの時点でのリストLを返却します。

・M.append((b, c))
 ここまでくれば、引数b,cは連分数の循環しない部分なので、タプル(b,c)をリストMに貯めます。

・C, t = (n-b*b)/c, int(n**(1.0/2))
 上記説明の(C)(E)を実装しました。
 **演算子は累乗で、平方根を求める部分を1/2乗として記述しました。
 整数で1/2とするとint型のため値が0になるので、1.0とfloat型にすることで
 1/2がfloat型の0.5になります。分母の2を2.0にしても同じことです。
 共通モジュールのmathのsqrt関数を使っても同じことです。
 import math
 t = math.sqrt(n)

・u = (b+t)%C
 上記説明の(H)を実装しました。

・L.append((b+t)//C)
 循環部分リストLには、上記説明の(G)を実装しました。

・return g(n, t-u, C, L, M)
 連分数展開を続けます。
 現世代のB,Cが次世代のb,cなので、上記説明の(C)(I)を実装します。


2.関数f(n)
・問題に合う、n以下で奇数の周期を持つ平方根の個数を返します。

・t = 0
 個数カウンタtを初期クリアしておきます

・for i in xrange(n+1):
 ループ変数iとして、0以上n+1未満(n以下)の整数を発生させます。

・s = g(i, int(i**(1.0/2)), 1, [], [])
 関数gで連分数展開したときの循環部分リストを取得しsとします。
 第1引数は連分数展開する平方根の元の値なので、ここではループ変数iです。
 第2引数は√iからくくりだせる最大整数なので、ここでは√iの整数部分です。
 第3引数には1、第4第5引数には空のリストを設定します。

・if len(s)%2: t += 1
 %演算子は割り算の余りです。
 len(s)は関数gで取得した連分数循環部分リストsの要素数です。
 この要素数が奇数ということは2で割った余りが1で、1は論理的にTrueなので、上記if文が成り立ちます。このとき、個数カウンタtを1つカウントアップします。

・return t
 変数iのループが終了したら、個数カウンタを呼び出し元に返します。

3.関数の外
・n = 10000
 print f(n)
 問題に合うように10000を関数fに渡し戻り値を表示します。
 
解答はこのすぐ下の行です。文字の色を白にしてます。選択状態にすると見えます。
1322


(追記)
・20121208 式(C)の分子分母が逆だったので修正。説明文の中に2箇所あり。
 誤)C = c/(n-b*b)・・・(C)
 ↓
 正)C = (n-b*b)/c・・・(C)
・20130127 ネタバレ注意の追記など

2012年11月2日金曜日

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

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

問63 「n 乗して得られる n 桁の正整数は何個あるか?」
5桁の数 16807 = 75は自然数を5乗した数である.
同様に9桁の数 134217728 = 89も自然数を9乗した数である.

自然数をn乗して得られるn桁の正整数は何個あるか?

-----

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






私の解答例は以下です。畳んでいます。
def f():
	import math
	return [(a,b,a**b) for a in xrange(1,10) for b in xrange(1, int(1/(1-math.log(a, 10)))+1)]

L = f()
print len(L),L


自然数aをb乗してb桁の正整数ということは、桁数の区切りは10の累乗値なので以下の式が成り立ちます。
0 ≦ 10(b-1) ≦ ab < 10b ・・・(A)
a>1、b>1 ・・・(B)

まず、aの変域を求めます。(A)の右部分の2式より、
1 ≦ a < 10 ・・・(C)

次にbの変域を求めます。
(A)の左部分をbについて解きます。
まず、左部分の2式の対数をとって、
(b-1)*log10 ≦ b*(log a)

次にbを左辺に寄せて、
(b-1)/b ≦ (log a)/log10

対数の底を10にして、
1-(1/b) ≦ log10a

1-log10a ≦ 1/b

b ≦ 1/(1-log10a)

(B)よりb>0なので、bの変域は
1 ≦ b ≦ 1/(1-log10a) ・・・(D)

上記(C)(D)の範囲で二重ループさせます。

2つの関数を使っています。

1.関数f()
・import math
 数学用関数モジュールmathをインポートして使えるようにしておきます。

・return [(a,b,a**b) for a in xrange(1,10) for b in xrange(1, int(1/(1-math.log(a, 10)))+1)]
 内包表記内に二重ループを入れています。
 まず前の「for a in xrange(1,10)」で、aが1以上10未満の整数でループします。
 そしてループ変数aを後ろのfor文へ渡します。
 後ろの「for b in xrange(1, int(1/(1-math.log(a, 10)))+1)」で、
 bが1以上1/(1-math.log10a)未満の整数でループします。
 math.log()関数は対数関数です。第2引数は対数の底です。
 int関数で切り捨てて整数にします。
 このa,bの値を先頭の式、(a,b,a**b)に渡します。
 これは問題に合うa, b, aのb乗の値のタプルです。
 このタプルを内包表記でリストにため、そのまま呼び出し元に返します。

2.関数の外
・L = f()
 print len(L),L
 関数f()で問題に合うa, b, aのb乗のタプルを格納したリストを取得しLとします。
 求める個数、とそのリストを表示します。

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

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

2012年10月31日水曜日

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

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

問62 「立方数になるような桁の置換をちょうど5つもつ最小の立方数を求めよ」
立方数 41063625 (3453) は, 桁の順番を入れ替えると2つの立方数になる:
56623104 (3843) と 66430125 (4053) である.
41063625は, 立方数になるような桁の置換をちょうど3つもつ最小の立方数である.
立方数になるような桁の置換をちょうど5つもつ最小の立方数を求めよ.


-----





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





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

def f(n):
	j = 0
	while True:
		j += 1
		L = [i*i*i for i in xrange(10**j)]
		M = [[str(s).count(str(t)) for t in xrange(10)] for s in L]
		for s in M:
			t = M.count(s)
			if t>=n:
				P = [i for i, u in enumerate(M) if s==u]
				N = [L[i] for i in P]
				return N[0], N, P

n = 5
print f(n)

最初は立方数の1つひとつの順列を作り、それが立方数かチェックしようとしましたが、1分では処理が終わりません。
そこで、立方数の1つひとつを構成する数字別個数リストを作って、同じリストになるものの数が指定個数になるかチェックするようにしました。

関数を1つ使っています。

1.関数f(n)
・j = 0
 変数jは10の累乗ごとに区切ってチェックするのでその累乗値カウンタです。初期クリアします。

・while True:
 無限ループです。

・j += 1
 累乗カウントをカウントアップします。

・L = [i*i*i for i in xrange(10**j)]
 Lは立方数リストです。10**jは10のj乗で、0からこの値までを1区切りにして処理を続けます。

・M = [[str(s).count(str(t)) for t in xrange(10)] for s in L]
 Mは数字別個数リストのリストです。二重の内包表記です。
 まず外側の内包表記の「for s in L」で先ほどの立方数リストLから1つずつ取り出しfor文の前の内側の内包表記に渡します。
 次に内側の内包表記の「for t in xrange(10)」で、ループ変数tに0以上10未満の整数が1つずつ設定されfor文の前に渡されます。
 「str(s).count(str(t))」で、立方数sに含まれる0以上10未満のいずれかのtの個数を数えます。
 例えば、s=1006012008の場合、内側の内包表記は[5,2,1,0,0,0,1,0,1,0]です。
 このような内側の内包表記が数字別個数リストで、外側の内包表記でそれをさらにリストにためています。

・for s in M:
 Mは数字別個数リストのリストから1つずつ取り出した数字別個数リストをループ変数sに設定します。

・t = M.count(s)
 リストMにある、ループ変数sの数字別個数リストと同じものの件数を数えtとします。

・if t>=n:
 上記tが指定個数n以上ならば、問題に合います。

・P = [i for i, u in enumerate(M) if s==u]
 Pは問題にあう立方数の元の数のリストです。
 enumerate関数で0からの連番を発生させます。このインデックス値が立方数の元の数です。

・N = [L[i] for i in P]
 リストPの値をインデックス値に立方数リストLから値を集め、リストNとします。

・return N[0], N, P
 問題に合う最小の立方数、桁を置換した立方数のリスト、その元の数のリストの3つを呼び出し元に返します。
 
2.関数の外
・n = 5
 print f(n)
 問題に合うように指定個数を5個として関数fで関連情報リストも含めて取得し、表示します。

解答はこのすぐ下の行です。文字の色を白にしてます。選択状態にすると見えます。(127035954683L,
[127035954683L, 352045367981L, 373559126408L, 569310543872L, 589323567104L],
[5027, 7061, 7202, 8288, 8384])

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

2012年10月30日火曜日

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

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

問61 「6つの巡回する4桁の図形数となる唯一の集合の和を求めよ」
三角数, 四角数, 五角数, 六角数, 七角数, 八角数は多角数であり, それぞれ以下の式で生成される.

三角数P3,n=n(n+1)/2 1, 3, 6, 10, 15, ...
四角数P4,n=n2 1, 4, 9, 16, 25, ...
五角数P5,n=n(3n-1)/2 1, 5, 12, 22, 35, ...
六角数P6,n=n(2n-1) 1, 6, 15, 28, 45, ...
七角数P7,n=n(5n-3)/2 1, 7, 18, 34, 55, ...
八角数P8,n=n(3n-2) 1, 8, 21, 40, 65, ...

3つの4桁の数の順番付きの集合 (8128, 2882, 8281) は以下の面白い性質を持つ.
  1. この集合は巡回的である. 最後の数も含めて, 各数の後半2桁は次の数の前半2桁と一致する
  2. それぞれ多角数である: 三角数 (P3,127=8128), 四角数 (P4,91=8281), 五角数 (P5,44=2882) がそれぞれ別の数字で集合に含まれている
  3. 4桁の数の組で上の2つの性質をもつはこの組だけである.
同じように, 6つの4桁の数からなる順番付きの集合で,
  1. 集合は巡回的である
  2. 三角数, 四角数, 五角数, 六角数, 七角数, 八角数が全て表れる唯一の集合を求め, その和を求めよ.

-----

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






私の解答例は以下です。たたんでいます。
def p3(n): return n*(n+1)//2
def p4(n): return n*n
def p5(n): return n*(3*n-1)//2
def p6(n): return n*(2*n-1)
def p7(n): return n*(5*n-3)//2
def p8(n): return n*(3*n-2)

def h(f, a, b):
	m = __import__(__name__)
	L = []
	for i in xrange(b):
		s = getattr(m, f)(i)
		if a<=s<b: L.append(s)
		elif b<=s: break
	return L

def g(L, M=[], n=2):
	if len(L)<=len(M): return M
	for s in L[len(M)]:
		if M and str(M[-1])[-n:]!=str(s)[:n]: continue
		if len(L)-len(M)>1:
			M = g(L, M[:]+[s], n)
			if len(L)<=len(M): break
			if s==L[len(M)][-1]: break
		elif str(s)[-n:]==str(M[0])[:n]: return M+[s]	
	else:
		M.pop()
	return M

def f(a, b, n=4, m=2):
	import itertools
	P = [h("p"+str(i), 10**(n-1), 10**n) for i in xrange(a, b+1)]
	for t in itertools.permutations(xrange(b+1-a)):
		L = [P[int(i)] for i in t]
		M = g(L, n=m)
		if M: break
	return sum(M), M, [str(i+a)+"kakusu" for i in t]

if __name__=="__main__":
	s, M, t = f(a=3, b=8, n=4, m=2)
	print s, M
	print t

角数の関数6個と他3個の関数を使っています。

1.関数p3(n)~p8(n)
・三角数から八角数の定義式のとおりです。

2.関数h(f, a, b)
・当pythonソースファイル内の関数fにb未満の値を渡し、a以上b未満の戻り値のリストを返します。

・m = __import__(__name__)
 __import__関数は、指定したモジュールを実行時に動的にインポートする関数です。
 _name__変数には、このプログラムのpythonソースのモジュール名(ファイル名の拡張子よりも前の部分)が設定されています。
 通常のimport関数はインポートするモジュール名を事前に固定文字列で指定しておきますが、自分のファイル自身をインポートする場合はファイル名変更しても影響を受けないように、動的にインポートします。
 後で関数を動的にインポートするための準備です。
 変数mには動的にインポートしたモジュールオブジェクトが設定されます。
 
・for i in xrange(b):
 ループ変数iは0以上b未満の範囲で1ずつ増加します。

・s = getattr(m, f)(i)
 getattr(m, f)関数は、モジュールmの関数fを呼び出します。
 変数iは関数fの引数です。
 当pythonソースファイル中の関数fに先ほどのループ変数iを渡し、その戻り値を変数sとします。

・if a<=s<b: L.append(s)
 elif b<=s: break
 上記の値sが、a以上b未満ならリストLに追加し、b以上なら終了します。

3.関数g(L, M=[], n=2)
・リストLは三角数などの角数リストが格納されたリストです。
 リストMはリストLの要素である角数リストを順にn桁で巡回的連結した値を格納したのリストです。
 リストLの全部の要素リストをn桁で巡回的連結した場合、その連結した値のリストを返します。

・if len(L)<=len(M): return M
 自分自身を呼び出す再帰関数なので、まず再帰呼び出し終了条件です。
 リストLの要素リストの個数以上に、順に巡回的連結した値の個数があれば順に巡回的に連結できた値のリストMを返し処理終了です。

・for s in L[len(M)]:
 len(M)は巡回的連結済みの個数であり、リストLのインデックス値としては次の巡回的連結要素を探すリストのインデックス値です。
 リストLから次に巡回的に連結する値を探す対象の角数リストを選び、ループ変数sとします。

・if M and str(M[-1])[-n:]!=str(s)[:n]: continue
 リストMに要素があり、その最後の要素がループ変数sにn桁で巡回的連結しなければsは対象外なので、次のループ変数sにループを進めます。
 リストMの最後の要素とループ変数sがn桁で巡回的連結している場合は以下の処理を続けます。

・if len(L)-len(M)>1:
  M = g(L, M[:]+[s], n)
  if len(L)<=len(M): break
 巡回的連結値の個数とリストLの要素数の差が1より大きい、つまり、リストLの最後までいっていない場合、さらに自分自身である関数g()を再帰呼び出しして、巡回的連結可能な値を探します。
 この再帰呼び出しのときの引数はリストLと巡回判定桁数nはそのままの値で、巡回的連結値リストMにはsを追加した状態で次を探します。
 再帰呼び出しから戻ってきたら、最初の行と同じ終了判定をします。

・elif str(s)[-n:]==str(M[0])[:n]: return M+[s]
 巡回的連結値の個数よりも、リストLの要素数が多くないということは、これらの個数が等しく、リストLの要素である角数リストすべてから巡回的連結値を探し出した状態なので、最後のチェックとして先頭要素との巡回的連結で連結の輪が閉じれるかを判定します。まだリストMに追加前のループ変数sからリストMの最後の要素にn桁で巡回的連結できればリストMが完成しこれを呼び出し元に返します。
 M+[s]は、リストMと要素sだけのリストの連結です。appendして返すことと同じです。

・else:
 M.pop()
 このelseはループ変数sのfor文のelseです。
 for文が途中でbreakせず、ループ変数の最後の要素まで処理し終った後に行う処理です。
 ここでは呼び出し元に返すこともなくfor文が終わったということなので、現在にリストMの最後の要素が巡回的連結先が無いことが確定しているので、このリストMの最後の要素をpop()関数で取り除きます。
 
・return M
 最後まで処理が来た場合は、リストMをそのまま返します。

4.f(a, b, n=4, m=2)
・n桁のa角数からb角数を1つずつ選びm桁で巡回的につながる組を返します。

・import itertools
 標準モジュールitertoolsを使えるようにしておきます。

・P = [h("p"+str(i), 10**(n-1), 10**n) for i in xrange(a, b+1)]
 リストPは、関数h()による角数リストをためたリストです。内包表記で書いています。
 for文でループ変数iにa以上b以下の整数を1つずつ設定します。
 ループ変数iは前に送られ、関数h()を呼び出しその戻り値を設定します。
 関数hの第1引数で、"p"+str(i)、つまり三角数などの関数、「p3」などを指定します。第2引数と第3引数で戻り値がn桁になるようします。例えばn=3の場合、100、1000を引き渡します。
 
・for t in itertools.permutations(xrange(b+1-a)):
 ループ変数tは、後ででてくるリストLのindex番号の組合せに使います。
 permutations()は順列を返す関数です。
 例えばrange(3)ならば、012, 021, 102, 120, 201, 210と順に返します。
 xrange(b+1-a)で、角数リストの個数分だけ0からの連番を発生させます。
 
・L = [P[int(i)] for i in t]
 上記tで先ほどのリストPの中を組み替えて、6つの角数リストを任意の順に組み替えます。
・M = g(L, n=m)
 if M: break
 関数g()を呼び出し、空のリスト(論理的にFalse)ならば、ループを続け次の
 順列の角数リストで巡回的連結可能な値を探しにいきます。

・return sum(M), M, [str(i+a)+"kakusu" for i in t]
 巡回的連結がした値の合計値、その値リスト、およびどの角数リストからとってきたかを呼び出し元に返します。

5.関数の外
・s, M, t = f(a=3, b=8, n=4, m=2)
 print s, M
 print t
 三角数から八角数まで各角数のうち4桁の値を1つずつ2桁で循環的連結できる値の和、そのときの値のリスト、どの角数リストからとってきたかのリストの3つを受け取り表示します。

解答はこのすぐ下の行です。文字の色を白にしてます。選択状態にすると見えます。
28684 [8256, 5625, 2512, 1281, 8128, 2882]
['3kakusu', '4kakusu', '7kakusu', '8kakusu', '6kakusu', '5kakusu']


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

2012年10月9日火曜日

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

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

問60 「任意の2つの素数を繋げたときに別の素数が生成される, 5つの素数の組を求めよ」
素数3, 7, 109, 673は非凡な性質を持っている.
任意の2つの素数を任意の順で繋げると, また素数になっている.
例えば, 7と109を用いると, 7109と1097の両方が素数である.
これら4つの素数の和は792である.
れは, このような性質をもつ4つの素数の組の和の中で最小である.

任意の2つの素数を繋げたときに別の素数が生成される, 5つの素数の組の和の中で最小のものを求めよ.


-----

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





最初にコーディングしたときにfor文の五重ループになってしまったので、どうしても再帰関数に改良したかったことと、処理速度の向上とで、1ヶ月近くかかってしまいました。
上限値を設けてやっと1分ルールを満たしているので、まだまだロジックが甘いのですが、行き詰まってしまったこともあり、現在のプログラムのまま提示しておきます。

答えは合っていますが、問題では「最小のもの」を求めるところ、題意を満たす「最初のもの」を求めるプログラムです。

私の解答例は以下です。畳んでいます。
def h(n):
	if n==2: return True
	if n<2 or (not n%2): return False
	for i in xrange(3, n, 2):
		if i*i>n: break
		if not n%i: return False
	return True

def e(n):
	i = n+1+int(n%2)
	while not h(i): i += 2
	return i

def g(L, t):
	t = str(t)
	for s in L:
		s = str(s)
		if not h(int(s+t)): return False
		if not h(int(t+s)): return False
	return True

def f(n, L=[], p=0):
	if len(L)>=n: return L
	while p<10**(n-1):
		p = e(p)
		if g(L, p):
			L = f(n, L+[p], p)
			if len(L)>=n: break
	else:
		L.pop()
	return L

n = 5
s = f(n)
print sum(s), s


4つの関数を使っています。

1.関数h(n)
・素数判定です。問58で作成したものを流用しました。
 素数ならTrue、素数でなければFalseを返します。

2.関数e(n)
・引数を超える最小の素数を返します。n>1のとき有効です。

・i = n+1+int(n%2)
 変数iが素数候補です。
 %演算子は割り算の余りです。割り切れれば0、割り切れないと1なので、
 式全体では偶数なら1足し、奇数なら2足すことになります。
 後続処理の準備として、素数候補iを奇数にしておきます。

・while not h(i): i += 2
 関数h()で素数判定をして、素数でなければ2ずつ増やしてチェックします。
 前準備でこの処理にくるときには素数候補iは奇数にしているので、
 iが偶数で無限ループになることはありません。

・return i
 ここまでの処理でnを超える最初の素数としてiを返却します。

3.関数g(L, t)
・引数リストLの各要素と引数tを任意の順に2つ連結した値が、すべて素数かどうかを判定します。

・t = str(t)
 後で連結するので、引数tを文字型にします。

・for s in L:
 引数リストLから順にループ変数sに設定します。

・s = str(s)
 後で連結するので、ループ変数sを文字型にします。

・if not h(int(s+t)): return False
 if not h(int(t+s)): return False
 上記の文字列sと文字列tを前後とその入れ替えの2パターンで連結し、
 int関数で数値化して、関数hで素数判定します。
 素数でなければ即終了でFalseを返します。

・return True
 引数リストLと引数tでここまでで処理を抜けていなければ、そのすべての組合せが素数なので、Trueを返します。

4.関数f(n, L=[], p=0)
・問題に合う値のリストを返します。
 自分で自分を呼び出す記述のある再帰関数です。
 引数nは求める素数の個数です。
 引数Lは素数を格納するリストで初期値は空リストです。
 引数pは最後にチェックした素数で、初期値は0とします。

・if len(L)>=n: return L
 再帰関数の終了条件です。
 リストLに求める個数n個以上の要素がたまったら処理終了で、リストLを返します。

・while p<10**(n-1):
 対象とする素数が青天井では処理が終わらないので、とりあえず求める素数の個数分の桁数、10**(n-1)、未満の素数で、処理することにしました。
 求める素数が2個なら10以内の素数の組で最小組[3,7]が見つかり、
 同じく3個なら100以内の素数の組で最小組[3, 37, 67]が見つかり、
 同じく4個なら問題文にあるように1000以内の素数の組が見つかります。

・p = e(p)
 得られている最後の素数pを関数eで処理し、次に大きい素数を改めてpとします。

・if g(L, p):
 関数gで引数リストLの各要素と引数pを任意の順に2つ連結した値が、すべて素数かどうかを判定します。Falseの場合はwhileループの次の値へ進みます。

・L = f(n, L+[p], p)
 1つ前の関数g()ですべて素数の場合、次に大きい素数pを候補として、さらに題意にある素数を探して追加するため、当関数f()を再帰呼び出しします。
 なお、第2引数の「L+[p]」は、再帰呼び出しで引数Lで受けるので、お試しでリストLの末尾に追加することになります。

・if len(L)>=n: break
 再帰呼び出しから戻り、リストLの要素数がn以上ならば探索終了でwhileループから抜けます。

・else:
 L.pop()
 このelse節はインデント(行頭のずらし位置)から明らかなように、for文のelseです。直前のif文のelse節ではありません。
 for文のelse節の処理タイミングは、forループの最後の値まで処理が行われた場合のループ終了直後です。
 ここまでの処理でbreakで抜けない場合、再帰呼び出しの際にお試しで追加したpが求める値ではないということになります。
 このため、関数pop()で末尾の要素を削除しておきます。

・return L
 処理の最後にリストLを呼び出し元に戻します。

5.関数の外
・n = 5
 s = f(n)
 print sum(s), s
 問題に合うように引数の値を5として関数f()を呼び出します。
 返されたリストをsとして、そのリスト要素合計、及び参考としてリスト要素を表示します。


 
解答はこのすぐ下の行です。文字の色を白にしてます。選択状態にすると見えます。
26033 [13, 5197, 5701, 6733, 8389]


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

2012年9月12日水曜日

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

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

問59「XOR暗号化された暗号文を解読できるか?」
各文字はそれぞれ一意のコードに割り当てられている.よく使われる標準としてASCII (American Standard Code for Information Interchange) がある.ASCIIでは, 大文字A = 65, アスタリスク (*) = 42, 小文字k = 107というふうに割り当てられている.

モダンな暗号化の方法として, テキストファイルの各バイトをASCIIに変換し, 秘密鍵から計算された値とXORを取るという手法がある.
XOR関数の良い点は, 暗号化に用いたのと同じ暗号化鍵でXORを取ると平文を復号できる点である.
65 XOR 42 = 107であり, 107 XOR 42 = 65である.

破られない暗号化のためには, 鍵は平文と同じ長さのランダムなバイト列でなければならない.ユーザーは暗号文と暗号化鍵を別々の場所に保存する必要がある.また, もし一方が失われると, 暗号文を復号することは不可能になる.
悲しいかな, この手法はほとんどのユーザーにとって非現実的である.

そこで, 鍵の変わりにパスワードを用いる手法が用いられる.パスワードが平文より短ければ (よくあることだが), パスワードは鍵として繰り返し用いられる.
この手法では, 安全性を保つために十分長いパスワードを用いる必要があるが, 記憶するためにはある程度短くないといけない.

この問題での課題は簡単になっている.
暗号化鍵は3文字の小文字である.
cipher1.txtは暗号化されたASCIIのコードを含んでいる.
また, 平文はよく用いられる英単語を含んでいる.

この暗号文を復号し, 平文のASCIIでの値の和を求めよ.

-----


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






私の解答例は以下です。畳んでいます。
def g(s, n):
	import itertools
	for t in itertools.product(s, repeat=n): yield "".join(t)

def f(sIn, a, n):
	L = [int(t) for t in open(sIn).read().split(",")]
	for s in g(a, n):
		M = ([ord(t) for t in s]*len(L))[:len(L)]
		N = [i^j for i,j in zip(L, M)]
		R = [chr(i) for i in N]
		t = "".join(R)
		if " the " in t: return sum(N), s, t

import os
sIn = os.path.join(os.path.dirname(__file__), "cipher1.txt")
a = "abcdefghijklmnopqrstuvwxyz"
n = 3
print f(sIn, a, n)


2つの関数を使っています。

1.関数g(s, n)
・文字列sからn個選んだ直積を連結した文字列を返します。
 直積とは、重複を許し順番が違えば別のものとして扱う組み合わせです。
 例えば、abcから2つ選んだ直積は、aa,ab,ac,ba,bb,bc,ca,cb,ccです。

・import itertools
 標準モジュールitertoolsを使えるようにしておきます。

・for t in itertools.product(s, repeat=n): yield "".join(t)
 標準モジュールitertoolsにあるproduct関数で、文字列sからをn個選んだ直積の構成要素のタプルを生成し、for文のループ変数tに設定します。
 forループの中で、タプルtを区切り文字無しで連結し、yield文で返します。
 yield文で返すと、返した時点で処理内容が一旦凍結され、再度呼び出されたときに処理を再開して続きの処理による値を返します。
 このような動きをするyield文で返す関数をジェネレータ関数といいます。
 なお、単発呼び出しで終了するreturn文で返す関数はイテレータ関数といいます。
 また、product関数はpython2.6以上です。
 python2.5にもitertoolsモジュールはありますが、product関数が含まれていません。

2.関数f(sIn, a, n)
・主処理です。
・入力ファイルsInの中の文字列を、指定した文字列aからn個選んだ秘密鍵sによって順にxor演算で復号した数値列N、及びこれをasciiコード扱いして得られる復号文字列tがあるとき、数値列Nの合計値、キー文字列s、復号文字列tを返します。

・引数sInは入力ファイルのフルパス、文字列aは秘密鍵の構成文字候補の文字列、nは秘密鍵の桁数です。

・L = [int(t) for t in open(sIn).read().split(",")]
 暗号文の各数値をもつリストLを作成します。
 入力ファイルsInを開き、その内容をすべて読み込んで、カンマで区切った文字列を
 ループ変数tに設定します。
 このtをfor文の前に送り、int関数で数値化してリストLとします。
 なお、開いたファイルオブジェクトはpythonがタイミングをみて閉じくれます。

・for s in g(a, n):
 関数gで文字列aからn個選んだ秘密鍵候補をループ変数sとします。

・M = ([ord(t) for t in s]*len(L))[:len(L)]
 Mは、秘密鍵候補sの構成要素の数値を暗号文の長さまで繰り返した復号用鍵のタプルです。
 内包表記のループ変数tとして秘密鍵候補sの文字1つずつを、for文の前に送り、ord関数でasciiコードとして対応する数値にしてリスト化します。
 このリストをlen(L)倍、つまり暗号文の長さ倍して、秘密鍵を順に繰り返し連結した状態にしてから、[:len(L)]で暗号文の長さの分だけ取り出しておきます。
 
・N = [i^j for i,j in zip(L, M)]
 Nは復号化した数値のリストです。
 zip関数で上記のリストLとリストMから1つずつ組み合わせて取り出し、ループ変数i,jとし、内包表記のfor文の前に送り、^演算子で数値j,jのxorを計算しリストにため、リストNとします。

・R = [chr(i) for i in N]
 復号数値リストNから1つずと取り出しchr関数でascii文字化してリストにため、リストRとします。
 
・t = "".join(R)
 上記の文字列リストRの要素を区切り文字無しで連結し、復号文tとします。

・if " the " in t: return sum(N), s, t
 復号が正しいか判定します。
 判定は「よく用いられる英単語」があるという問題文に従い、半角空白ではさんだ「 the 」としてみました。結果的にこれで判定できました。
「 the 」があれば、復号化した数値リストNの要素の合計、秘密鍵、復号文tを返します。

3.関数の外
・import os
 ファイルを扱うので、標準モジュールosを使えるように取り込んでおきます。

・sIn = os.path.join(os.path.dirname(__file__), "cipher1.txt")
 暗号文のフルパスです。
 「__file__」は実行するpythonファイルのフルパスです。
 os.path.dirname関数で親フォルダまで取得し、ファイル名とともに、
 os.path.joinでパス表示を連結しました。
 つまり、暗号文は実行するpythonファイルと同じフォルダに置いて実行するように記述しています。
 
・a = "abcdefghijklmnopqrstuvwxyz"
 秘密鍵の候補です。問題文から小文字のアルファベットです。

・n = 3
 秘密鍵の桁数です。問題文から3としました。

 

解答はこのすぐ下の行です。文字の色を白にしてます。選択状態にすると見えます。
107359
秘密鍵はgod。平文の内容は、聖書ヨハネによる福音書の第1章第1節から第14節。

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

2012年9月9日日曜日

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

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

問58「渦巻きに配置された格子の対角線上にある素数の数を調べ上げよ」
1から始めて, 以下のように反時計回りに数字を並べていくと,
辺の長さが7の渦巻きが形成される.


37 36 35 34 33 32 31
38 17 16 15 14 13 30
39 18  5  4  3 12 29
40 19  6  1  2 11 28
41 20  7  8  9 10 27
42 21 22 23 24 25 26
43 44 45 46 47 48 49


面白いことに, 奇平方数が右下の対角線上に出現する.
もっと面白いことには, 対角線上の13個の数字のうち, 8個が素数である.
ここで割合は8/13 ≒ 62%である.


渦巻きに新しい層を付け加えよう. すると辺の長さが9の渦巻きが出来る.
以下, この操作を繰り返していく.


対角線上の素数の割合が10%未満に落ちる最初の辺の長さを求めよ.

-----


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





私の解答例は以下です。畳んでいます。
def h(n):
	if n<2 or (not n%2): return False
	for i in xrange(3, n, 2):
		if i*i>n: break
		if not n%i: return False
	return True

def g(n):
	if n<1 or (not n%2): return []
	elif n==1: return [1]
	s, t = n-1, n-2
	u = t*t+s
	return [u, u+s, u+s*2, u+s*3]

def f(p):
	i, u = 1, 0
	while not (0<float(u)/(4*(i//2)+1)<p):
		i += 2
		u += len([j for j in g(i) if h(j)])
	return i

p = 0.1
print f(p)


3つの関数を使っています。

1.関数h(n)
・素数判定です。2と、3以上の奇数で割り切れるか判定します。

・if n<2 or (not n%2): return False
 2以下、または2で割り切れる場合、素数ではありません。
 %演算子は割り算の余り、余り0は論理的にFalseです。

・for i in xrange(3, n, 2):
 素因数候補として、3以上の奇数を順にループ変数iに設定します。
 割り切り判定の割る数として使用します。
 終了値をnにしていますが、後続ロジックでnになる前に必ずループを終了します。

・if i*i>n: break
 引数nを構成する素因数の限界値は√nです。これを超える素因数はあり得ません。そこで、素因数候補がこの限界値まで達したら、ループを抜けます。

・if not n%i: return False
 割り切り判定です。割り切れたら、素数でないのでFalseを返します。
 %演算子は割り算の余り、余り0は論理的にFalseです。

・return True
 割り切れる素因数候補が無く、素数なのでTrueを返します。

2.関数g(n)
 問28では同じ渦巻きの四隅和を返す関数を作りましたので、
 これの返り値を四隅和から四隅値のリストに修正しました。

3.関数f(n)
・問題文の渦巻きについて、対角線上の値のうち素数の割合が引数p未満になる最初の辺の長さを返します。

・i, u = 1, 0
 初期設定です。iは渦巻きの辺の長さ、uは対角線上の素数の個数です。

・while not (0<float(u)/(4*(i//2)+1)<p):
 「float(u)/(4*(i//2)+1)」は、対角線上の素数の割合です。
 この値が0と引数pの間に入ったら、ループ終了です。
 uは対角線上の素数の個数で、割り算で小数点以下まで求めるので、float関数で浮動小数点型にしています。
 「(4*(i//2)+1)」は対角線上の値の個数です。辺の長さiにより決まります。

・i += 2
 辺の長さは渦巻きを1つ巻くごとに、2ずつ増えます。

・u += len([j for j in g(i) if h(j)])
 内包表記です。
 for文で関数g(i)にて取得した辺の長さiの四隅値の1つひとつをループ変数jに設定し、後ろのif文に渡し、関数h(j)にて素数判定し素数だけをfor文の前に渡します。
 for文の前はjをリスト化し、len関数でその要素数を数えて、uに累積します。

・return i
 whileループ終了時点の辺の長さiを呼び出し元に返します。

4.関数の外
・p = 0.1
 print f(p)
 問題に合うように10%=0.1を関数fに渡し、その返り値を表示します。
 
解答はこのすぐ下の行です。文字の色を白にしてます。選択状態にすると見えます。
26241

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