分割数 - inamori’s diary

分割数というのは、例えば5を自然数で分割することを考えます。同じ数を何度も使ってもよいですが、順序違いは区別しません。そうすると、

5 = 5
  = 4 + 1
  = 3 + 2
  = 3 + 1 + 1
  = 2 + 2 + 1
  = 2 + 1 + 1 + 1
  = 1 + 1 + 1 + 1 + 1

の7通りあります。これを

p(5) = 7

と書いて5の分割数と呼びます。

漸化式

これをどう計算しましょうか。漸化式になりそうですね。最初に5を取ったらそこで終わり、4を取ったらp(1)、3を取ったらp(2)。2を取ったときは気をつけなければなりません。残りが3ですから、p(3)とすればよさそうですが、これだと次に3を取る場合も含まれてしまい、最初に3を取る場合と重複してしまいます。

5 = 3 + 2 = 2 + 3

そのため、2を取ったら次以降は2以下しか取れないなどとします。そうするとp(3)ではなく、3の2以下による分割数になります。これをp(3, 2)と書くことにしましょう。そうすると漸化式は次のように書けます。

p(5) = p(0, 5) + p(1, 4) + p(2, 3) + p(3, 2) + p(4, 1)

一般的には、

p(0, m) = 1
p(n, m) = p(n, n) (n < m)
p(n, m) = p(n - m, m) + p(n - m + 1, m - 1) + … + p(n - 1, 1)

となります。2番目の式は、例えばp(2, 3)は3個一度には取れないということです。これでコードを書いてみましょう。

memo = { }
def p(n, m1 = 0):
    m = n if m1 == 0 else min(n, m1)
    
    if n == 0:
        return 1
    elif (n, m) in memo:
        return memo[(n,m)]
    else:
        s = sum(p(n - d, d) for d in xrange(1, m + 1))
        memo[(n,m)] = s
        return s

メモ化のコストを除いてO(n³)程度ですね。少し工夫すると速くなります。

p(5, 5) = 1 + p(5, 4)

5を分割するのに、5だけか4以下のみ使って分割するかに分けています。一般的には、

p(n, 1) = 1
p(n, m) = p(n - m, m) + p(n, m - 1)

ですね。これでO(n²)になっています。

母関数

分割数の母関数は次のようになります。

G(x) = (1 + x + x² + …)(1 + x² + x⁴ + …)(1 + x³ + x⁶ + …)…

展開して、xⁿの係数がp(n)になります。こう書いたらわかりやすいかもしれません。

G(x) = (1 + x¹ + x¹⁺¹ + x¹⁺¹⁺¹ + …)(1 + x² + x²⁺² + …)…

例えば、2 + 1 + 1 + 1は

x¹⁺¹⁺¹x²

に対応します。

def mul_poly(f, g):
	h = [ 0 ] * (N + 1)
	for e1, c1 in enumerate(f):
		for e2 in g:
			e = e1 + e2
			if e > N:
				break
			h[e] += c1
	return h

N = 500
fs = [ range(0, N + 1, n) for n in xrange(1, N + 1) ]
p = reduce(mul_poly, fs, [ 1 ])
print p[N]

オイラーの五角数定理

1 + x + x² + …は1/(1 - x)と書けますね。だから、

1/G(x) = (1 - x)(1 - x²)(1 - x³)…

と書けます。この右辺はどうなるでしょうか。

N = 20
fs = [ [ 1 ] + [ 0 ] * (n - 1) + [ -1 ] for n in xrange(1, N + 1) ]
p = reduce(mul_poly, fs, [ 1 ])
print p

[1, -1, -1, 0, 0, 1, 0, 1, 0, 0, 0, 0, -1, 0, 0, -1, 0, 0, 0, 0, 0]

やたら0が多くてあとは1か-1ですね。実はオイラーの五角数定理というものがあって、

$\prod_{n=1}^{\infty}{(1-x^n)} = \sum_{n=-\infty}^{\infty}{(-1)^nx^{n(3n-1)/2}}$

が成り立ちます。あとはこの逆数を取るだけです。

G(x) = p(0) + p(1)x + p(2)x² + …

として、

G(x)(1 - x - x² + x⁵ + x⁷ - x¹2 - …) = 1

だから、

p(n) = p(n - 1) + p(n - 2) - p(n - 5) - p(n - 7) + p(n - 12) + …

となります。1/G(x)のnまでの0でない項の個数はO(√n)だから、p(n)の計算量はO(n^1.5)になります。

from itertools import *
import time

t0 = time.clock()
N = 10000
p = [ 0 ] * (N + 1)
p[0] = 1
for n in xrange(1, N + 1):
	for k in count(1):
		m1 = k * (3 * k - 1) / 2
		m2 = k * (3 * k + 1) / 2
		if n < m1:
			break
		p[n] += p[n-m1] if k % 2 == 1 else -p[n-m1]
		if n < m2:
			break
		p[n] += p[n-m2] if k % 2 == 1 else -p[n-m2]
print p[N]
print time.clock() - t0

最後に、4通りのnの分割数を求める方法の速度を測ってみました。

n	漸化式	漸化式2	母関数	五角数定理
500	1.60s	80ms	140ms	7ms
1000	14.8s	422ms	626ms	18ms
2000	-	2.03s	3.13s	51ms
5000	-	73.6s	22.3s	205ms
10000	-	-	-	605ms

Project Eulerに必要な用語集