Project Euler 44（2）

少し数学的にアプローチしましょう。
求めるDを、

D = P_l = P_k - P_j

とすると、

l(3l - 1) = k(3k - 1) - j(3j - 1)
l(3l - 1) = (k - j)(3k + 3j - 1)

k - j と 3k + 3j - 1 は l(3l - 1) の約数になります。ただし、制約がいろいろあります。3の因子はすべて k - j が受け持つことになります。また、k - j と k + j は偶奇が同じでなければなりません。よって、k - j と 3k + 3j - 1 は偶奇が反対になります。すなわち、2の因子はどちらかがすべて受け持つことになります。だから、l(3l - 1) の素因数分解をして、2と3の因子を取り去った数の約数をすべて出せばよいです。このとき、3k + 3j - 1 は k - j の3倍より大きくなければなりません。

具体的に、l = 12 で考えてみましょう。2P₁₂ = 12 • 35 = 2² • 3 • 5 • 7 なので、35の約数を出せばよいです。

約数1	約数2	k - j	3k + 3j - 1	mod3 = 2?	3倍以上?	k + j	k	j
1	35	12	35	○	×	-	-	-
		3	140	○	○	47	25	22
5	7	60	7	×	×	-	-	-
		15	28	×	×	-	-	-
7	5	84	5	×	○	-	-	-
		21	20	×	○	-	-	-
35	1	420	1	×	×	-	-	-
		105	4	×	×	-	-	-

k = 25, j = 22のみが P₁₂ = P_k - P_l を満たすことが分かります。
最後に、P_j + P_k が五角数かどうか判定します。

from itertools import takewhile, count

from math import sqrt

import time
def is_pentagonal(n):

    p = 1 + 24 * n

    q = int(sqrt(p))

    return q * q == p and q % 6 == 5
def P(n):

    return n * (3 * n - 1) / 2
def gen_prime_candidates():

    yield 2

    yield 3

    n = 5

    while True:

        yield n

        yield n + 2

        n += 6
def index(f, p):

    for k in range(len(f)):

        if f[k][0] == p:

            return k

    return -1
def pop_p(f, p):

    k = index(f, p)

    if k == -1:

        return f, 0

    else:

        return f[:k] + f[k+1:], f[k][1]
def calc_exp(n, p):

    e = 0

    while n % p == 0:

        e += 1

        n /= p

    return e, n
def factorize(n):

    f = ()

    for p in takewhile(lambda p: p * p <= n, gen_prime_candidates()):

        e, n = calc_exp(n, p)

        if e:

            f = f + ( (p, e),)

    

    if n > 1:

        f = f + ( (n, 1),)

    return f
def gen_divisors(f, k = 0):

    if k == len(f):

        yield 1, 1

    else:

        p, e0 = f[k]

        for d1, d2 in gen_divisors(f, k + 1):

            for e in range(e0 + 1):

                yield d1 * p ** e, d2 * p ** (e0 - e)
def gen_divs(l):

    f = factorize(l) + factorize(3 * l - 1) # coprime

    f, e2 = pop_p(f, 2)

    f, e3 = pop_p(f, 3)

    n2 = 2 ** e2

    n3 = 3 ** e3

    for d1, d2 in gen_divisors(f):

        if e2 == 0:

            yield n3 * d1, d2

        else:

            yield n2 * n3 * d1, d2

            yield n3 * d1, n2 * d2
def solve():

    for l in count(1):

        for d1, d2 in gen_divs(l):

            if d1 * 3 < d2 and d2 % 3 == 2:

                d3 = (d2 + 1) / 3

                n = (d1 + d3) / 2

                m = d3 - n

                if is_pentagonal(P(n) + P(m)):

                    return P(l)
t = time.clock()

print solve()

print time.clock() - t