Project Euler 94

プロジェクトオイラー
http://projecteuler.net/index.php

Q94.
それぞれの辺の長さが整数で、辺の差が1以内のものを「ほぼ正三角形」と呼ぶ。面積も整数になるもので、周囲の長さが10億以内の「ほぼ正三角形」の周囲の長さの総和

辺の長さをそれぞれ、a, b, cとすると、面積Sはヘロンの公式より、

S = √s(s-a)(s-b)(s-c)
s = (a + b + c) / 2

2辺が同じ長さで、1辺が1差だから、どちらかが奇数でどちらかが偶数だが、簡単に2辺のほうが奇数であることがわかる。
2n-1と2nのとき、s=3n-1, S=n√(3n-1)(n-1)となるから、

3n2-4n+1=k2

(3n-2)2 - 3k2 = 1

2n+1と2nのとき、s=3n+1, S=n√(3n+1)(n+1)となるから、

3n2+4n+1=k2

(3n+2)2 - 3k2 = 1

よって、どちらの場合もペル方程式、

x2 - 3y2 = 1

を解けばよい。この解は、

am + bm√3 = (2 + √3)m

x = am, y = bm

と表される。最初の解のxは3で割って2余るから、3n+2でn=0だから、辺の長さが0になる。その次からは面積は正になる。

sn = an + 1 (an ≡ 1(mod 3))

sn = an - 1 (an ≡ 2(mod 3))

N = 10 ** 9
def gen_solution():

    a = 2

    b = 1

    while True:

        a, b = a * 2 + b * 3, a + b * 2

        if a % 3 == 1:

            s = a + 1

        else:

            s = a - 1

        if s * 2 > N:

            break

        yield s * 2
print sum(gen_solution())