Project Euler: Problem 23

問題はこちら。毎回、答え合わせに使っているProject Euler: Python solutionsよりも高速にできた。
方針としては

1. 28123以下の過剰数のリスト abd を求める。
2. 1.で作ったリスト abd を使い、1から28123までの全ての数について2つの過剰数の和で書けるか調べる。
3. 2つの過剰数で書けない数の和を求める。

まずは、真の約数のリストを返すdivisorsを定義。（本当は真の約数の和を返せば十分だけど、後々使えそうなので）

def divisors(n):
    divs = [1]
    for i in xrange(2, int(math.sqrt(n)) + 1):
        if n % i == 0:
            divs.append(i)
            if n/i not in divs:
                divs.append(n/i)
    divs.sort()
    return divs

はじめ、メインの処理は以下のように書いていた。sum_of_two は後で説明するが、2つの過剰数の和で書けるかをチェックする関数。

abd = [i for i in xrange(28124) if sum(divisors(i)) > i]
print(sum([i for i in xrange(1, 28124) if not sum_of_two(i, abd)]))

sum_of_twoは自然数 n と過剰数のリスト abd を受け取って、True/False を返す。はじめ何も考えずに書いたのがこれ。

def sum_of_two(n, abd):
    for i in abd:
        if i > n / 2:
            break
        if n - i in abd:
            return True
    return False

これでも動く（はず）なのだけれど、5分待ってもメインの処理が終わらなかった。いくらなんでも遅すぎ。
実は、5行目の n - i のメンバシップ確認を 'in' で書いてしまうと、毎回 abd を順番に走査しはじめるので遅い。最悪、len(abd) 回のチェックが必要になってしまう。abd のメンバシップ確認が一番多く行われるので、ここは毎回 O(1) の手間でメンバシップ確認できるようにしておくべきだ。というわけで、メイン処理にキャッシュを作る処理を追加。

cache = [False] * 28124
for i in xrange(1, 28124):
    if sum(divisors(i)) > i:
        cache[i] = True
abd = [i for i in xrange(28124) if cache[i]]
print(sum([i for i in xrange(1, 28124) if not sum_of_two(i, abd)]))

sum_of_two もキャッシュを利用することに。

def sum_of_two(n, abd):
    for i in abd:
        if i > n / 2:
            break
        if cache[n - i]:
            return True
    return False

これだと4秒で答えの 4179871 が出て終わる。