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

プロジェクトオイラーの問題を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