Project Euler 12

PE 数学 プログラミング

Problem 12
（略）約数の個数が500を超えるはじめての三角数を求めよ。
http://projecteuler.net/index.php?section=problems&id=12

約数の個数は、素因数分解をすれば簡単に求められますが、個別に素因数分解していては非常に遅くなります。エラトステネスのふるいを使えば比較的速く求められますが、ある程度まとまっていないと（あるいは等間隔で並んでいないと）いけません。ところが、

Tn = n(n + 1) / 2

とバラバラになるので、その手は使えません。しかし、まとめて素因数分解して約数の個数を求めて、隣同士でそれを掛け合わせば、n(n + 1)の約数の個数になります（nとn + 1は互いに素、すなわち同じ素因数を持たないので、この掛け算が成り立ちます）。偶数はあらかじめ2で割っておくと、Tnの約数の個数が求められます。T9までの約数の個数を求めておきましょう。



1 1 3 2 5 3 7 4 9 5

1 1 1 1 1 1 1 1 1 1

上の段が元の数で、下の段が約数の個数です。上は1〜10で、偶数なら2で割っておきます。上の段の隣同士を掛け合わせると三角数になることを確認しておきましょう。



1 3 6 10 15 21 28 36 45

まず、2で割れるだけ割ります。例えば、4は、2で2回割れるので、下の段を(2+1)倍します。



1 1 3 1 5 3 7 1 9 5

1 1 1 2 1 1 1 3 1 1

同様に3で割ります。



1 1 1 1 5 1 7 1 1 5

1 1 2 2 1 2 1 3 3 1

上の段で1でないのは素数なので、下の段を2倍します。



1 1 1 1 1 1 1 1 1 1

1 1 2 2 2 2 2 3 3 2

下の段の隣同士を掛ければ、Tnの約数の個数となります。



1 2 4 4 4 4 6 9 6

いつ約数の個数が500より大きくなるかわからないので、約数の個数はある程度まとめて計算して溜めておき、出し尽くしたらまたまとめて計算する、というのを繰り返します。ここでは1000個ずつまとめて計算しています。



from itertools import count, takewhile

from math import sqrt
def make_primes(s, e):

    def is_prime(n):

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

            if n % p == 0:

                return False

        return True

    

    if s == 2:

        s = 3

    for n in range(s, e, 2):

        if is_prime(n):

            primes.append(n)
def sum_pow(n, p):

    e = 0

    while n % p == 0:

        e += 1

        n /= p

    return e, n
def int_sqrt(n):

    return int(sqrt(n + 0.5))
def gen_sum_divs():

    for m in count():

        s = N * m + 1

        e = N * (m + 1) + 1

        make_primes(int_sqrt(s) + 1, int_sqrt(e) + 1)

        a = [ n if n & 1 else n / 2 for n in range(s, e) ]

        b = [ 1 ] * N

        for p in takewhile(lambda p: p * p <= e, primes):

            for n in xrange( (s + p - 1) / p * p, e, p):

                l = n - s

                q, a[l] = sum_pow(a[l], p)

                b[l] *= q + 1

        

        for k in xrange(N):

            if a[k] != 1:

                b[k] *= 2

            n = k + s

            yield n, b[k]
M = 500

N = 10 ** 3

primes = [ 2 ]

prev_s = 1

for n, s in gen_sum_divs():

    if prev_s * s > M:

        print n * (n - 1) / 2

        break

    prev_s = s


ところで、この問題は逆からは攻められないのでしょうか。すなわち、約数の個数が500より大きい自然数を小さい順から羅列し、それが三角数であるか判定する、という方法です。これがなかなかうまくいかないのです。単に一番小さい自然数を求めよ、というのなら比較的簡単です。小さい自然数は、

2e13e25e3…
e1 ≥ e2 ≥ e3…

という形になるからです。



def gen_young(n, m0 = 1, d0 = 0):

    if d0 == 0:

        d0 = n

    

    for d in range(2, d0 + 1):

        m = m0 * d

        if n > m:

            for t in gen_young(n, m, d):

                yield (d, ) + t

        else:

            yield (d, )

        if n <= m:

            break
def value(e):

    return reduce(lambda x, y: x * y,

            map(lambda k: primes[k] ** (e[k] - 1), range(len(e))))
N = 500

primes = [ 2, 3, 5, 7, 11, 13, 17, 19, 23, 29 ]

print min(map(value, gen_young(N))) # 14414400 = 2^6*3^4*5^2*7*11*13

しかし、小さい順に列挙するのは難しいです。