http://projecteuler.net/index.php?section=problems&id=143
∠BTC = ∠CTA = ∠ATB = 120°より、
a2 = q2 + qr + r2
b2 = p2 + pq + q2
c2 = r2 + rp + p2
が成り立ちます。個々の整数解はピタゴラス数と同じように求められます。
x = q / a
y = r / a
と置き換えて、
x2 + x y + y2 = 1
グラフを描くと、
斜めになった楕円です。(0, 1)を必ず通るので、この点から第1象限で再び楕円と交わるように直線を引きます。
y = -k x + 1
その条件は、1 / 2 < k < 1 です。楕円との交点は、
x2 + x(-k x + 1) + (-k x + 1)2 = 1
(1 - k + k2)x2 + (1 - 2k)x = 0
x = (2k - 1) / (1 - k + k2)
y = (1 - k2) / (1 - k + k2)
交点の座標が有理数 ⇔ kが有理数 です。
k = n / m(m / 2 < n < m)として、分母を払って、
(q, r, a) = (m(2n - m), m2 - n2, m2 - m n + n2)
mとnが素のとき、m + nが3の倍数だとq, r, aの全て3の倍数になります。その場合、
(m', n') = ((m + n) / 3, (2m - n) / 3)
とすれば、q, r, aがそれぞれ1/3になるので（qとrはひっくり返ります）無視できます。
さて、このように必要なq, rの組を求めたら、これをグラフにします。すなわち、q, rを点と考えて、qとrを辺で結びます。そして3点でループができていれば、その3点がp, q, rとなって題意を満たしていることになります。
from itertools import * from fractions import gcd def gen_qr(max_qr): for m in xrange(2, (max_qr + 1) / 2 + 1): g = takewhile(lambda n: n * (2 * m - n) <= max_qr, xrange(m / 2 + 1, m)) for n in ifilter(lambda n: (m + n) % 3 != 0 and gcd(m, n) == 1, g): q = m * (2 * n - m) r = m * m - n * n if q < r: q, r = r, q for l in xrange(1, max_qr / (q + r) + 1): yield l * q, l * r def make_graph(max_n): def add(q, r): if q in g: g[q].append(r) else: g[q] = [ r ] g = { } for q, r in sorted(list(gen_qr(max_n))): add(q, r) add(r, q) return g def gen_triangles(g): for q, nodes in g.iteritems(): # let q be the smallest greater_node = filter(lambda x: x > q, nodes) if len(greater_node) >= 2: for r, p in combinations(greater_node, 2): if g[r].count(p) == 1: yield p + q + r N = 120000 graph = make_graph(N) s = set(n for n in gen_triangles(graph) if n <= N) print sum(s)
