Project Euler 9

Problem 9
ピタゴラス数は、a < b < cで、a2 + b2 = c2を満たす3つの自然数の組である。
(略)a + b + c = 1000を満たすピタゴラス数は唯一つ存在する。このときのabcを求めよ。
http://projecteuler.net/index.php?section=problems&id=9

もっとも素直に書いたらこうでしょうか。


N = 1000
for c in range(3, N / 2 + 1):
for b in range(2, c):
a = N - (b + c)
if a < b and a * a + b * b == c * c:
print a, b, c, a * b * c

しかし、ピタゴラス数には隈なく生成する公式があるので、これを使うのが常套手段です。

a = 2lmn b = l(m2 - n2) c = l(m2 + n2)
m > n
mnは互いに素
m + nは奇数

これによると、

a + b + c = 2lm(m + n)

なので、500を素因数分解して、積が500になる3つの自然数の組を生成して、上の条件を満たすものを挙げればよいです。



from fractions import gcd

def get_exp(n, p):
e = 0
while n % p == 0:
e += 1
n /= p
return e, n

def gen_prime_candidates(n):
yield 2
for p in range(3, n, 2):
yield p

def factorize(n):
f = [ ]
for p in gen_prime_candidates(n):
if p * p > n:
break
e, n = get_exp(n, p)
if e:
f.append( (p, e) )
if n != 1:
f.append( (n, 1) )
return f

def mul(x, y):
return x * y

def value(f):
return reduce(mul, map(lambda t: t[0] ** t[1], f), 1)

def add_fac(f, p, e):
if e == 0:
return f
else:
for k in range(len(f)):
if f[k][0] == p:
return f[:k] + (p, f[k][1] + e) + f[k+1:]
return f + ((p, e),)

def gen_divisors(f, l, k = 0):
if l > 2:
for f1, f2 in gen_divisors(f, 2):
for t in gen_divisors(f2, l - 1):
yield (f1, ) + t
else:
if k == len(f):
yield ( (), () )
else:
p, e = f[k]
for f1, f2 in gen_divisors(f, 2, k + 1):
for e1 in range(e + 1):
yield add_fac(f1, p, e1), add_fac(f2, p, e - e1)

def gen_divs(f, l):
for t in gen_divisors(f, l):
yield tuple(map(value, t))

def gen_solutions(N):
f = factorize(N / 2)
for l, m, p in gen_divs(f, 3):
n = p - m
if gcd(m, p) == 1 and 0 < n and n < m and p % 2 == 1:
a = l * (m * m - n * n)
b = l * 2 * m * n
c = l * (m * m + n * n)
if a > b:
a, b = b, a
yield a, b, c

N = 1000
for a, b, c in gen_solutions(N):
print a, b, c, a * b * c