Project Euler 94

Problem 94
辺の長さが整数で面積も整数の正三角形が存在しないことは簡単に示される。しかし、「ほぼ正三角形」の5-5-6という三角形は12平方単位の面積を持つ。
「ほぼ正三角形」を2つの辺の長さが等しく第三の辺はそれと1単位より違わない三角形と定義しよう。
辺の長さが整数で面積も整数で周囲の長さが10億を超えない「ほぼ正三角形」の周囲の長さの総和を求めよ。
http://projecteuler.net/index.php?section=problems&id=94

等しい辺をa、第三の辺を2bとしましょう。
いつものようにピタゴラス数の生成の式を使います。考えられるのは次の2つのパターンです。

(a, b) = (m2 + n2, m2 - n2)
(a, b) = (m2 + n2, 2mn)

上の場合、

2b - a = m2 - 3n2 = ±1

2乗して3で割って2余る整数はないので、符号は+であることがわかります。そうすると、これはペル方程式なので簡単です。周囲の長さlは、

l = 2a + 2b = 4m2

下の場合、

2b - a = 4mn - m2 - n2 = ±1
(m + n)2 - 3(m - n)2 = ±2

これもペル方程式に近いです。上と同じように符号は-であることがわかります。
さて、

x2 - 3y2 = -2

ですが、(x, y) = (1, 1)が解であることがわかります。すると、

xn - yn√3 = (1 - √3)(2 - √3)n

で与えられる(xn, yn)も解であることがわかるでしょう(共役の数を掛け合わせれば-2になる)。どうやらこれで解は尽くされているようです。
周囲の長さlは、

l = 2a + 2b = 2(m + n)2

from itertools import takewhile

def gen_Pell(m, n):
    while True:
        m, n = m * 2 + n * 3, m + n * 2
        yield m, n

N = 10 ** 9
def less_eq(p): return p <= N
print sum(takewhile(less_eq, (m * m * 4 for m, n in gen_Pell(1, 0)))) \
    + sum(takewhile(less_eq, (m * m * 2 for m, n in gen_Pell(1, 1))))