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(m2 - n2) b = 2l m n c = l(m2 + n2)
として、
c / (a - b) = k
となる整数kがある。
m2 + n2 = k(m2 - n2 - 2l m n)
(k - 1)m2 - 2k m n - (k + 1)n2 = 0
m / n = (k ± √(2k2 - 1)) / (k - 1)
平方根の中は平方数でなければならないので、
2k2 - 1 = p2
p2 - 2k2 = -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()))