Project Euler 45（2）

P_m = H_n を式変形すると、

m(3m - 1) = 2n(2n - 1)
(36m² - 12m) / 12 = (16n² - 8n) / 4
( (6m - 1)² - 1) / 12 = ( (4n - 1)² - 1) / 4
(6m - 1)² - 3(4n - 1)² = -2

結局、次の方程式に帰着されました。

x² - 3y² = -2

これはペル方程式に似ています。
(x, y) = (1, 1) が最小解であることはすぐにわかります。

p_n - √3q_n = (2 - √3)ⁿ

としたとき、

x_n - √3y_n = (p_n - √3q_n)(1 - √3)

とすると、(x_n, y_n) が元の方程式の解であることが分かります。

x_n² - 3y_n²
= (x_n - √3y_n)(x_n + √3y_n)
= (p_n - √3q_n)(1 - √3)(p_n + √3q_n)(1 + √3)
= -2(2 - √3)ⁿ(2 + √3)ⁿ
= -2

ただし、一般的にはほかに解があるかもしれません。今回の問題については無いようです。ですから、x_nが6で割って余り5、y_nが4で割って余り3になる解を選べばよいわけです。

from itertools import count, islice
def is_square(n):

    m = int(sqrt(n))

    return m * m == n
def gen_Pell_solution():

    n, m = 1, 1

    while True:

        yield n, m

        n, m = 2 * n + 3 * m, n + 2 * m
def gen_both():

    for x, y in gen_Pell_solution():

        if x % 6 == 5 and y % 4 == 3:

            m = (x + 1) / 6

            n = (y + 1) / 4

            yield m * (3 * m - 1) / 2
print sum(islice(gen_both(), 2, 3))