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

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

問７５「直角三角形を作るときに1通りの折り曲げ方しか存在しない鉄線の長さの総数」
ある長さの鉄線を折り曲げたときに1通りの直角三角形を作る最短の長さは12 cmである.
他にも沢山の例が挙げられる.

12 cm: (3,4,5)
24 cm: (6,8,10)
30 cm: (5,12,13)
36 cm: (9,12,15)
40 cm: (8,15,17)
48 cm: (12,16,20)

それとは対照的に, ある長さの鉄線 (例えば20cm) は整数長さで折ることができない.
また2つ以上の折り曲げ方があるものもある. 2つ以上ある例としては, 
120cmの長さの鉄線を用いた場合で, 3通りの折り曲げ方がある.

120 cm: (30,40,50), (20,48,52), (24,45,51)

Lを鉄線の長さとする. 
直角三角形を作るときに1通りの折り曲げ方しか存在しないような
L ≦ 1,500,000 の総数を答えよ.

－－－－－


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

def gcd(a, b):
	#if a<b: a, b = b, a
	while b: a, b = b, a%b
	return a

def f(n):
	M = [0 for i in xrange(1, n+1)]
	for i in xrange(1, n):
		for j in xrange(1, i):
			if 2*i*(i+j)>n: break
			a, b, c = i*i-j*j, 2*i*j, i*i+j*j
			if gcd(c, a)>1: continue
			if gcd(c, b)>1: continue
			s = a+b+c
			for k in xrange(s, n, s): M[k] += 1
	L = [i for i, j in enumerate(M) if j==1]
	return len(L)

n = 1500000
print f(n)

直角三角形の直角を挟む辺をa、b、斜辺をcとして、
ピタゴラスの定理「a2 + b2 = c2」が成り立つかを総当りでチェックすると、
１分ルールに全然間に合いません。

まず考えたことは、最大公約数が1の既約なa,b,cの組を探し、この整数倍の組を消しこむ方法です。
例えば(3,4,5)に基づいて(6,8,10)などをチェックし探すことから外します。

次に考えたことは偶数奇数の絞込みです。
・aもbも偶数の場合
　(2m)2 + (2n)2 = c2から、c2は偶数、つまりcも偶数になってしまい、
　a,b,cは既約でないので、これは成立しません。
・aもbも奇数の場合
　(2m+1)2 + (2n+1)2
= 4m2 + 4m + 4n2 + 4n + 2
= 4(m2 + n2 + m + n) + 2  = c2
ここでc2が4で割ったときに余り2になるかというと、
　cが偶数ならば、c=2s、c2 = (2s)2 = 4s2、なのでc2は4で割りきれ余り0
　cが奇数ならば、c=2s+1、c2 = (2s+1)2 = 4(s2+s)+1、なのでc2は４で割ると余り1
つまり、c2は4で割って余り2になることはありえないので、これも成立しません。
よって、a,bは奇数と偶数の組です。偶数同士や奇数同士ではありません。
そこで、a2 + b2 = c2 より、偶数の２乗と奇数の２乗の和を求めると、
　(2m)2 + (2n+1)2
= 4m2 + 4n2 + 4n + 1
= 4(m2 + n2 +n) + 1  = c2
２乗して奇数ということは、cは必ず奇数です。

以上より直角を挟む辺a,bは奇数と偶数の組で、斜辺cは奇数です。
a,bは入れ替えても一般性を失わないので、a奇数、b偶数で固定します。

さらに、a2 + b2 = c2より、
b2 = c2 - a2
b2 = (c+a)(c-a) ...(1)

a, b, cは既約で、(1)の左辺が平方数なので、
c+a、c-aは両方とも平方数を因子に持ちます。 ...(2)

また、bは偶数なので、b=2rとおくことができて
b2 = (2r)2 = 4r2 なので、
b2は4の倍数 ...(3)

a, cは両方とも奇数なので、c+aとc-aは両方とも偶数。
c+a=2p, c-a=2q ...(4)

上記(2), (3), (4)を総合すると、0<j<iの条件で
b2 = 4i2j2
c+a = 2i2
c-a = 2j2
とおくことができる。よって、
b = 2ij
(c+a)+(c-a) = 2i2 + 2j2
         2c = 2(i2 + j2)
          c = i2 + j2
(c+a)-(c-a) = 2i2 - 2j2
         2a = 2(i2 - j2)
          a = i2 - j2
以上より、
(a, b, c) = (i2-j2, 2ij, i2+j2), 0<j<i, cが斜辺。...(5)


１．関数gcd(a, b)
・ユークリッドの互除法でa,bの最大公約数を返します。
・問５では自分で自分自身を呼び出す再帰呼び出しをしていましたが、while文を使うことで再帰呼び出しなしにしました。
・必ずa>bとなるように呼び出しすことで大小を入れ替える部分をコメントアウトしました。

２．関数f(n)
・長さn以下の針金で各辺が整数の直角三角形を作るときに1通りの折り曲げ方しか存在しない長さの個数を返します。
・折り曲げパターン数リストMは、i番目の要素に周がiで各辺が整数の直角三角形が何通りあるかを格納します。初期値はすべての要素が0です。

・for i in xrange(1, n):
for j in xrange(1, i):
if 2*i*(i+j)>n: break
a, b, c = i*i-j*j, 2*i*j, i*i+j*j
　上記(5)のまま、プログラミングしました。
　停止条件は３辺の和がnを超えたときなので、
　a+b+c = (i2-j2)+2ij+(i2+j2) = 2i2+2ij = 2i(i+j) > n です。

・if gcd(c, a)>1: continue
　if gcd(c, b)>1: continue
　a, b, cが既約でない、つまり最大公約数が1より大きい場合は、次の候補へ。

・s = a+b+c
　for k in xrange(s, n, s): M[k] += 1
　a, b, cの３辺和をsとして、sからsごとの値kが、sの整数倍の長さの針金であり、直角三角形の折り曲げ方が1つあるということで、折り曲げパターン数リストMのk番目の値を+1します。
　
・L = [i for i, j in enumerate(M) if j==1]
　return len(L)
　折り曲げパターン数リストMは、針金長ごとの折り曲げパターン数を格納してあるので、その要素が1のものだけの場合、リストLに格納します。内包表記で記述しています。
　そして、len関数でリストLの要素数を求め呼び出し元へ返します。


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