整数を分割する場合の数

いわゆる分割数ではありません。分割して0も認めます。また順番が変わっただけの分割も別とみなします。すなわち、3を3分割なら、

3 = 3 + 0 + 0 = 2 + 1 + 0 = 2 + 0 + 1 = 1 + 2 + 0 = 1 + 1 + 1 = 1 + 0 + 2
  = 0 + 3 + 0 = 0 + 2 + 1 = 0 + 1 + 2 = 0 + 0 + 3

全部で10通りあります。
これは母関数を使うと簡単に求められます。まず、1分割はどの整数でも1通りとみなします。各項の係数をその指数の分割の場合の数として、

P₁(x) = 1 + x + x² + x³ + ...

2分割はこの積で求められます。

P₂(x) = P₁(x)P₁(x) = 1 + 2x + 3x² + 4x³ + ...

例えば、3を分割するとき、3 = 0 + 3 = 1 + 2 = 2 + 1 = 3 + 0ですが、

x⁰x³ + x¹x² + x²x¹ + x³x⁰

に対応しています。
このようにすると、100の10分割も簡単に求められます。

def mul(f, g):
    h = [ 0 ] * (N + 1)
    for k in xrange(N + 1):
        for l in xrange(N - k + 1):
            h[k+l] += f[k] * g[l]
    return h

def pow(f, e):
    if e == 1:
        return f
    
    g = pow(f, e / 2)
    if e % 2 == 1:
        return mul(mul(g, g), f)
    else:
        return mul(g, g)

N = 100
M = 10
f = [ 1 ] * (N + 1)
print pow(f, M)[N]

しかし、例えば10000の100分割の場合の数はこのコードではなかなか返ってきません。
このようなときは、Problem 31の方法が使えます。すなわち、n分割のときの母関数の係数は指数のn-1次関数になるので、0〜n-1の分割の場合の数を上の方法で求め、それを用いて母関数の係数の関数の形を求めます。そうすれば10000だろうが1000000だろうが簡単に分割の場合の数を求めることができます。

from itertools import *
from fractions import Fraction

def mul(f, g):
    h = [ 0 ] * (M + 1)
    for k in xrange(M + 1):
        for l in xrange(M - k + 1):
            h[k+l] = (h[k+l] + f[k] * g[l])
    return h

def add(f, g):
    h = f[:]
    for k in xrange(M + 1):
        h[k] += g[k]
    return h

def pow(f, e):
    if e == 1:
        return f
    
    g = pow(f, e / 2)
    if e % 2 == 1:
        return mul(mul(g, g), f)
    else:
        return mul(g, g)

def unit():
    return [ 1 if k == 0 else 0 for k in xrange(M + 1) ]

def value(f, x):
    return reduce(lambda p, q: p * x + q, reversed(f), 0) % D

def calc_coeffs(f):
    def calc_num(k):    # numerator
        g = reduce(mul, imap(lambda l: [ -l, 1 ] + [ 0 ] * (M - 1),
                                ifilter(lambda l: l != k, xrange(M))))
        return map(lambda x: x * f[k], g)
    
    def calc_den(k):    # denominator
        return reduce(lambda x, l: x * (k - l),
                            ifilter(lambda l: l != k, xrange(M)), 1)
    
    def calc_fraction(k):
        den = calc_den(k)
        return map(lambda x: Fraction(x, den), calc_num(k))
    
    return reduce(add, imap(calc_fraction, xrange(M)))

N = 10000
M = 100
D = 10 ** 8
f = [ 1 ] * (M + 1)
fM = pow(f, M)
cs = calc_coeffs(fM)
print value(cs, N)