Project Euler 139

http://projecteuler.net/index.php?section=problems&id=139

単にピタゴラス数で条件に合う組合せを数えても1分くらいでした。

from itertools import *
from fractions import gcd

def gen_primitive_Pythagorean(L):
    for m in takewhile(lambda m: m * (2 * m + 1) <= L, count(2)):
        n0 = m + 2 if m % 2 == 1 else m + 1
        for n_ in (n for n in xrange(n0, 2 * m, 2) if gcd(m, n) == 1):
            a = n_ * (m * 2 - n_)
            b = m * (n_ - m) * 2
            c = m * m * 2 - a
            yield (a, b, c)

L = 10 ** 8
counter = 0
for a, b, c in gen_primitive_Pythagorean(L - 1):
    if c % (a - b) == 0:
        print a, b, c, (L - 1) / (a + b + c)
        counter += (L - 1) / (a + b + c)
print counter

もう少しなんとかならないでしょうか。

a = l(m² - n²) b = 2l m n c = l(m² + n²)

として、

c / (a - b) = k

となる整数kがある。

m² + n² = k(m² - n² - 2l m n)
(k - 1)m² - 2k m n - (k + 1)n² = 0
m / n = (k ± √(2k² - 1)) / (k - 1)

平方根の中は平方数でなければならないので、

2k² - 1 = p²
p² - 2k² = -1

これもペル方程式ですね。ただし、kが負の場合も考慮する必要があります。また、m / nは通分できる可能性があるので、必ずしも直角三角形が小さいほうから並ぶとは限らないはずです。しかし、実際にやってみると小さい順に並んでいます。通分での除数もだんだん大きくなります。

from itertools import *
from fractions import Fraction

def unfold(f, s):
    while True:
        x, s = f(s)
        yield x

def gen_Pells():
    def mul(p, q):
        x1, y1 = p
        x2, y2 = q
        return (x1 * x2 + 2 * y1 * y2, x1 * y2 + y1 * x2)
    
    a = (1, 1)
    r = mul(a, a)
    return unfold(lambda s: (s, mul(s, r)), a)

def gen_perimeters():
    def gen_mn():
        for p, k in gen_Pells():
            if k != 1:
                f = Fraction(p + k, k - 1)
                yield f.numerator, f.denominator
            f = Fraction(p + k, k + 1)
            yield f.numerator, f.denominator
    
    return (2 * m * (m + n) for m, n in gen_mn() if (m + n) % 2 == 1)

N = 10 ** 20
print sum((N - 1) / l for l in 
        takewhile(lambda l: l < N * 2, gen_perimeters()))