プロジェクトオイラー 問９２

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

問９２「桁の2乗による連鎖」
各桁の2乗を足し合わせて新たな数を作ることを, 同じ数が現れるまで繰り返す.

例えば

44 → 32 → 13 → 10 → 1 → 1
85 → 89 → 145 → 42 → 20 → 4 → 16 → 37 → 58 → 89

のような列である. どちらも1か89で無限ループに陥っている. 
驚くことに, どの数から始めても最終的に1か89に到達する.

では, 10,000,000より小さい数で89に到達する数はいくつあるか.





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






私の回答例は以下の通りです。

def f(n):
	P = [i*i for i in xrange(10)]
	s = len(str(n))
	m = 9*9*s
	
	L = [0]+[-1]*(m-1)
	L[1], L[89] = 0, 1
	for i in xrange(1, m):
		j = i
		while L[j]<0:
			j, r = 0, j
			while r:
				j += P[r%10]
				r //= 10
			L[i] = L[j]
			
	M = [0]*m
	for i in xrange(10): M[P[i]] = 1
	for u in xrange(2, s):
		W = [0]*m
		for i in xrange(m):
			for j in xrange(10):
				if P[j]<=i: W[i] += M[i-P[j]]
		M = W[:]
		
	return sum([i*j for i, j in zip(L, M)])
	
n = 10000000
print f(n)

1から10,000,000まで、桁の2乗和の連鎖を総当りで計算すると、とても1分ルールに合いません。

試行錯誤したところ、桁2乗和は１回計算するだけで値がずいぶんちいさくなることがわかりました。
桁数sのときの桁2乗和最大値mは9がs個並んだ数で、m = 9*9*s です。
7桁ならば、m = 9*9*7 = 567 です。

そこで、
・m以下の数は桁2乗和の連鎖を最後まで計算して89ループになるか判定しておく。
・桁2乗和は１回だけ計算しm以下にしたところで数える。
ことにしました。

・89ループ判定結果リストL
1からmまでの桁2乗和の連鎖を最後まで行い、89ループになる場合は1、その他は0を設定します。
-1で初期設定しているので、while L[j]<0として未設定が続く限り連鎖をたどります。
rのwhileループで、rの10で割った余りとして、下から1桁ずつ取り出して2乗和を累積します。
2乗の計算はいちいち行わず、2乗リストP[0,1,4,9,16...]を参照します。

・桁2乗和ごとの件数リストM
初期値として1桁数値の桁2乗和の位置(1,4,9...)に1を設定します。
2桁以上の数値については、１つ少ない桁数までの累積件数の、最上位の数値の2乗分だけ小さい位置にある値を、採用し累積していきます。
このようにしてだんだん桁数を増やしていきます。

こうして、
・89ループ判定結果リストL（89ループに達する数の位置に1が立ち、他は0が設定）
・桁2乗和を１回だけ計算した値ごとの件数リストM
の同じ位置の要素同士の積を求めることで89ループの部分を取り出します。
この総和が求める件数です。


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

プロジェクトオイラー 問９１

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

問９１「整数座標における直角三角形の個数」
点P(x1, y1)と点Q(x2, y2)はともに整数係数の点であり, 原点O(0,0)と合わせてΔOPQをなす.

各係数が0と2の間にあるとき, すなわち0≦x1, y1, x2, y2≦2のとき, 直角三角形は14個できる:

では, 0≦x1, y1, x2, y2≦50のとき, 直角三角形は何個作れるか?





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






私の回答例は以下の通りです。

def f(n):
	c = 0
	L = [(i, j) for i in xrange(n+1) for j in xrange(n+1)][1:]
	for x1, y1 in L:
		for x2, y2 in L:
			if x1==x2 and y1==y2: continue
			if not (x1*(x1-x2)+y1*(y1-y2)) \
			or not (x2*(x1-x2)+y2*(y1-y2)) \
			or not (x1*x2+y1*y2): c += 1
	return c//2

n = 50
print f(n)

２点P，Qの座標を総当りで作り出します。
0からnまでの数の重複を許す組み合わせをP, Qの座標候補リストLとします。
Lは内包表記に二重ループを入れてみました。
ただし、三角形なので３点O,P,Qはそれぞれ別の座標なので、まず点O(0,0)を除きます。
P(x1,y1),Q(x2,y2)のそれぞれの座標を二重ループで取得します。
まず点P,Qが同じ座標ならば三角形にならないので次の座標へ処理を進めます。

直角三角形か判定は、ベクトルの内積が0になることを利用します。
点Pが直角：ベクトルOPとベクトルQPの内積、x1*(x1-x2)+y1*(y1-y2) = 0
点Qが直角：ベクトルOQとベクトルQPの内積、x2*(x1-x2)+y2*(y1-y2) = 0
点Oが直角：ベクトルOPとベクトルOQの内積、x1*x2+y1*y2 = 0
数字の0は論理的にFalseなので、(not 内積値)で直角三角形であることを判定します。
行末の「\」は改行してもまだ同じ行と見なす継続行であることを示します。


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