Project Euler 122

http://projecteuler.net/index.php?section=problems&id=122

べき乗の計算はいわゆるバイナリ法を使うと速いですが、実はバイナリ法より乗算の回数が少ない方法があるという話です。
例えば5乗なら、

a -> a² -> a⁴ -> a⁵
a -> a² -> a³ -> a⁵
a -> a² -> a³ -> a⁴ -> a⁵

という計算手順があります。これらを、

(1, 2, 4, 5)
(1, 2, 3, 5)
(1, 2, 3, 4, 5)

と表します。最初の手順から一手で

a⁵ -> a⁶
a⁵ -> a⁷
a⁵ -> a⁹
a⁵ -> a¹⁰

このように全ての手順を求めて、最小の乗算の回数の和を取ります。

N = 20
a = [ [], [(1,)] ] + [ [] for k in xrange(N - 1) ]
for x in (x for n in xrange(1, N) for x in a[n]):
    for n in (n for n in (m + x[-1] for m in x) if n <= N):
        a[n].append(x + (n,))

print sum(min(len(x) - 1 for x in a[n]) for n in xrange(1, N + 1))

しかしこれは遅いです。乗算の回数が多すぎる手順まで記憶してしまうからです。ただ、どこで打ち切ればいいのかがよくわかりません。とりあえず、ここでは半分から計算しています。例えば、a^7なら、(a^3)^2 * aだからa^3の回数+2を最小とします。これでは過剰のように思えますが、正しい答えが出ます。

from itertools import *

def calc_max_m(c, n):
    if n % 2 == 0:
        return c[n/2] + 1
    else:
        return c[n/2] + 2

N = 200
M = next(n for n in count() if 1 << n >= N)
a = [ [] for k in xrange(N + 1) ]
a[1] = [ (1,) ]
c = [ 0 ]
for n in xrange(1, N + 1):
    c.append(min(len(x) - 1 for x in a[n]))
    max_m = calc_max_m(c, n)
    for x in a[n]:
        for m in x:
            n = m + x[-1]
            if n <= N and len(x) <= calc_max_m(c, n):
                a[n].append(x + (n,))

print sum(c[n] for n in xrange(1, N + 1))