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

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

問９５「友愛鎖」
ある数の真の約数とは, それ自身を除く約数すべてである. 例えば, 28 の真の約数は 1, 2, 4, 7, 14 である. これらの約数の和は 28 に等しいため, これを完全数と呼ぶ.
面白いことに, 220 の真の約数の和は 284 で, 284 の真の約数の和は 220 となっており, 二つの数が鎖をなしている. このため, 220 と 284 は友愛数と呼ばれる.

さらに長い鎖はあまり知られていないだろう. 例えば, 12496 から始めると, 5 つの数の鎖をなす.
12496 → 14288 → 15472 → 14536 → 14264 (→ 12496 → ...)
この鎖は出発点に戻っているため, 友愛鎖と呼ばれる.

いずれの要素も 1,000,000 を超えない最長の友愛鎖の最小のメンバーを求めよ.






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






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

def d(n):
	c, m = 1, int(n**.5)
	for i in xrange(2, m+1):
		if not n%i:
<			c += i
<			j = n//i
<			if i<j: c += j
	return c

def f(n):
	L, N = [d(i) for i in xrange(n+1)], []
	for i in xrange(n+1):
		j, M = L[i], []
		while (j not in M) and (j<=n):
<			M.append(j)
<			j = L[j]
		if M and i==M[-1] and len(N)<len(M): N = [i]+M[:-1]
	return len(N), min(N), N

n = 1000000
n,m,L = f(n)
print "min:",m
print "num:",n
print "chain:",L

１．関数d(n)
・関数d(n)はnの友愛数を返します。問２１を改良しました。
・数値nの「真の約数」の上限は、その定義からnの正の平方根です。
　n**(0.5)　が、0.5乗ということで、nの正の平方根です。
・変数cに真の約数を累積していきます。
　//演算は割り算の商（整数）で、%演算は割り算の余りです。
　変数iのforループで、nがiで割り切れるとき、変数cにiを累積します。
　割り切れるときは余り0で、%演算の結果がFalseになるので、not演算子でこのような場合を捉えます。
　またこのときの商j=n//iも約数なので、i<jならこれも変数cに累計します。
　nが平方数でiがその平方根の場合、i=jなのでiだけ累積しjは累積しません。

２．関数f(n)
・関数f(n)は、n以下の友愛数リストを取得し、その友愛数をたどり、最長の友愛鎖を返します。
・リストLは、その位置のオフセット値の友愛数を持つリストです。例）L[220]=284
・リストNには、最長の友愛鎖を格納します。
・ループ変数iは、友愛鎖をたどる最初の数値です。
・変数jは次の友愛数が格納されます。
・リストMには、数値iから始まる友愛鎖を格納します。
・while文のループ条件は、
　友愛数jが友愛鎖リストMになく、上限値以下であることとしています。
・友愛鎖リストMを最長友愛鎖リストNに保存する条件は、
　リストMに値があり、友愛鎖たどりの最初の数値iが友愛鎖の最後の値と一致し、
　リストNよりも長い（最長の友愛鎖）こととしています。

解答はこのすぐ下の行です。文字の色を白にしてます。選択状態にすると見えます。
min: 14316
num: 28
chain: [19116, 31704, 47616, 83328, 177792, 295488, 629072, 589786, 294896, 358336, 418904, 366556, 274924, 275444, 243760, 376736, 381028, 285778, 152990, 122410, 97946, 48976, 45946, 22976, 22744, 19916, 17716, 14316]

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

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

問９４「擬正三角形」
一辺の長さが整数の正三角形は面積が整数にならないことを示すのは簡単である. 
しかし, 5-5-6の辺を持つ殆ど正三角形に近い擬正三角形 (almost equilateral triangle) は面積が12で整数である.
以降, 二等辺三角形で, 3つめの辺の長さが他と1つしか違わないもの (5-5-6, 5-5-4等) を, 擬正三角形と呼ぶ.

さて, 周囲の長さが1,000,000,000以下の面積が整数になる擬正三角形を考え, その周囲の長さの総和を求めよ.






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






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

def f0(n):
	s, c = n//3, 0
	d = {}; dd = {}
	for i in xrange(1, s):
		ii = i*i
		d[i] = ii; dd[ii] = 0
	for i in xrange(3, s, 2):
		t = i//2; u = d[i]
		if (u-d[t]  ) in dd: c += i*3-1
		if (u-d[t+1]) in dd: c += i*3+1
	return c

n = 10**9
print f0(n)

１．関数f(n)
・関数f(n)は周囲の長さがn以下で面積が整数になる擬正三角形の周囲の長さの和を返します。
・変数sは１辺の最大値、変数cは周囲の長さの累積です。
・擬正三角形の二等辺に挟まれた角からもう１辺に垂線hを下ろすと、
　垂線の足が直角である、同じ面積の直角三角形２つに分割できます。
　（以下、半分直角三角形と呼びます）
　ここで半分直角三角形の３辺は以下のように置きます。
　h：擬正三角形の二等辺で挟まれた角からの垂線。
　a：擬正三角形の二等辺であり、半分直角三角形の斜辺。
　b：hと直角を構成する。擬正三角形の底辺の半分の長さ。

・最初のforループで二乗値の準備をします。
　辞書dには二乗値を格納し、辞書ddには整数二乗値をキーにその正の平方根を格納します。

・２つ目のforループは、ループ変数iが擬正三角形の二等辺の長さです。
　これは半分直角三角形の斜辺でもあり、問７５より、３辺が整数の直角三角形は、
> 直角を挟む辺a,bは奇数と偶数の組で、斜辺cは奇数です。
　なので、ループ変数iは3から１辺最大値sまで2飛びで値を増加させます。
・変数tは二等辺の長さを2で割って切り捨てることで、
擬正三角形の３辺がi,i,i-1の場合の半分直角三角形の直角を挟む１辺の長さです。
そして、t+1が擬正三角形の３辺がi,i,i+1の場合の半分直角三角形の直角を挟む１辺の長さになります。

・二乗値辞書dから取得した、上記iとtの二乗値の差は、三平方の定理より
　擬正三角形の垂線hの長さの二乗値と等しく、この値が辞書ddのキーにあれば整数直角三角形なので、擬正三角形の周囲の長さを変数cに累積していきます。

・・・と、ここまできて、１分ルールが守れず長考に入ってしまいました・・・。

def f2(n):
	r, s, t, x, y = 0, 0, 0, 1, 0
	while s<=n or t<=n:
		x, y = 2*x+3*y, x+2*y
		
		s = x*x*4
		if s<=n: r += s
		
		t = x+3*y
		t = 2*t*t
		if t<=n: r += t
		
	return r

n = 10**9
print f2(n)

改めて、問７５より、３辺が整数の直角三角形では、整数i,j（0<j<i）を使って、
> 斜辺がi2+j2、直角を挟む２辺はi2-j2, 2ij

よって、上記の半分直角三角形に当てはめると、
a = i2+j2　：(A)
(b, h) = (i2-j2, 2ij) or (2ij, i2-j2)　：(B)
擬正三角形の周長L = 2(a+b)　：(C)
擬正三角形の定義から、a = 2b±1　：(D)

・ケース１：b = i2-j2の場合
(D)に(A)と上記を代入して
i2 + j2 = 2(i2 - j2)±1
i2 + j2 = 2i2 - 2j2 ±1
i2 - 3j2 = ±1　：(E)
また、(C)に(A)と上記を代入して
L = 2((i2+j2)+(i2-j2))
  = 4i2　：(F)

・ケース２：次に上記(B)で、b = 2ijの場合
(D)に(A)と上記を代入して
i2 + j2 = 2・2ij ±1
i2 + j2 = 4ij ±1
(i - 2j)2 - 3j2 = ±1　：(G)
また、(C)に(A)と上記を代入して
L = 2(i2+j2 +2ij)
  = 2(i+j)2　：(H)
ここで、t=i-2jと置くと、(G)、(H)は、
t2 - 3j2 = ±1　：(G')
L = 2((t+2j) +j)2
  = 2(t+3j)2　：(H')
　
・上記の２つの場合をまとめると、
x2 - 3y2 = ±1　：(J)
L = 4x2, 2(x+3y)2　：(K)

・さてここで、（J)の「±1」の「-1」の場合を考えてみます。
x2 - 3y2 = -1
x2 = 3y2 -1
xの二乗は、３の倍数に１足りないので、
xの二乗を３で割ると２余ることになります。
x = 3m, 3m+1, 3m+2　として確認していきます。
x = 3mの場合、x2 = (3m)2 = 6m、3で割ると余り0
x = 3m+1の場合、x2 = (3m+1)2 = 9m2+6m+1、3で割ると余り1
x = 3m+2の場合、x2 = (3m+2)2 = 9m2+12m+4、3で割ると余り1
上記より、二乗して３で割ったときに２余る整数は存在せず、式(J)は以下に改めます。
x2 - 3y2 = 1　：(J')

式(j')はペル方程式と呼ばれる式で、以下の公式があります。
x2 - Dy2 = 1　（Dは平方数以外。平方数の場合、整数解が存在しない）
の解は、最小整数解を(x0, y0)、それ以降、(x1, y1),(x2, y2), ...として、
┌xn+1 = x0xn + y0Dyn
└yn+1 = y0xn + x0yn

当問題に当てはめると、D=3、また、、
問題にある5-5-6の擬正三角形が最小整数解で、このときx0=2, y0=1 なので、
┌xn+1 = 2*xn + 3*yn
└yn+1 = 1*xn + 2*yn

上記式をループ内に実装し、上記(K)の擬正三角形の周の長さL = 4x2, 2(x+3y)2　を計算し累積します。




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