Project Euler 151から（3）

元のコードはこうでしたが、

from itertools import izip, count

def next(s, k):
    return s[:k] + (s[k] - 1,) + tuple(n + 1 for n in s[k+1:])

def E(s):
    num_papers = sum(s)
    if num_papers == 0:
        return 0.0
    
    return (1 if num_papers == 1 else 0) \
            + sum(n * E(next(s, k)) / num_papers
                    for k, n in izip(count(), s) if n > 0)

N = 5
print "%.6f" % (E((1,) + (0,) * (N - 1)) - 2)

E(s)は何回通るでしょうか。実行時間はこれによると思われます。
上のコードで通る回数を数えればよいですが、これでは上のコードと時間がかかります。なんとかならないでしょうか。
例えば、状態(0, 0, 1)は、(0, 1, 0)と(0, 0, 2)から移行します。だから、(0, 1, 0)と(0, 0, 2)が通った回数を足せば、(0, 0, 1)が通った回数になります。その際に次のコードのgen_sのようにすると、封筒を取り出した回数順に状態を出すことができるので、容易に通った回数を計算することができます。

from itertools import chain

def add(dic, key, value):
    if key in dic:
        dic[key] += value
    else:
        dic[key] = value

def next(s, k):
    return s[:k] + (s[k] - 1,) + tuple(n + 1 for n in s[k+1:])

def gen_s(set_s):
    if len(set_s) > 0:
        return chain(set_s, gen_s(set(next(s, k) for s in set_s
                                for k in xrange(len(s)) if s[k] > 0)))
    else:
        return set_s

def num_pass(N):
    s0 = (1,) + (0,) * (N - 1)
    dic_num = { s0 : 1 }
    for s in gen_s(set((s0,))):
        n = dic_num[s]
        for k in (k for k in xrange(len(s)) if s[k] > 0):
            add(dic_num, next(s, k), n)
    return sum(dic_num.itervalues())

for N in xrange(2, 9):
    print N, num_pass(N)

これを実行すると、

2 3
3 8
4 125
5 191170
6 7087287305719
7 404634508398297036890780530296
8 151348294904659598414501102210976841112270718960049871705594229385

ものすごい勢いで増加しますね。どうりで最初のコードではN=6のときに結果が返ってこないわけです。