Project Euler 121

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


4回の場合、全体の場合の数は2*3*4*5、4勝する場合の数は1*1*1*1、3勝する場合の数は1回目に負ける場合の数は1*1*1*1、2回目なら1*2*1*1、3回目1*1*3*1、4回目1*1*1*4なので、勝つ確率は11/120となります。
これをコード化すると、

from itertools import *

def factorial(n):
    return 1 if n == 0 else n * factorial(n - 1)

def mul(a):
    return reduce(lambda x, y: x * y, a, 1)

N = 15
print factorial(N + 1) / sum(mul(a) for n in xrange((N + 1) / 2)
                        for a in combinations(range(1, N + 1), n))

負け数で分けて組み合わせを出して場合の数を計算して和を取っています。
これはいかにも効率が悪いので、もうちょっとなんとかしましょう。こういう問題は再帰に持ち込んでメモ化と相場が決まっています。n回残っていてm回以上勝たなければならないときの場合の数をメモ化しています。

def factorial(n):
    return 1 if n == 0 else n * factorial(n - 1)

def num_cases(n, m):
    x = (n, m)
    if x in memo:
        return memo[x]
    
    if n == 0:
        y = 1
    elif m <= 0:
        y = (n + 1) * num_cases(n - 1, 0)
    elif n == m:
        y = num_cases(n - 1, m - 1)
    else:
        y = num_cases(n - 1, m - 1) + n * num_cases(n - 1, m)
    
    memo[x] = y
    return y

N = 15
memo = { }
print factorial(N + 1) / num_cases(N, N / 2 + 1)

もっとエレガントな手法があります。それは母関数を使う方法です。

G(x) = (1 + x)(1 + 2x)(1 + 3x)(1 + 4x)
= 1・1・1・1 + (1・1・1・1 + 1・2・1・1 + 1・1・3・1 + 1・1・1・4)x + ...

定数項が0敗のとき、xの係数が1敗のときの場合の数になっています。

def factorial(n):
    return 1 if n == 0 else n * factorial(n - 1)

def mul(f, g):
    h = [ 0 ] * (len(f) + len(g) - 1)
    for i in range(len(f)):
        for j in range(len(g)):
            h[i+j] += f[i] * g[j]
    return h

N = 15
f = reduce(lambda x, y: mul(x, [1, y]), xrange(1, N + 1), [1])
print factorial(N + 1) / sum(f[k] for k in xrange((N - 1) / 2 + 1))