Puzzle DE Programming / 分割数

パズルの説明

整数 n を 1 以上の自然数の和で表すことを考えます。これを「整数の分割」といいます。整数を分割するとき、同じ自然数を何回使ってもかまいませんが、並べる順序が違うだけのものは同じ分割とします。簡単な例を示しましょう。次の図を見てください。

─┬─ 6                           : 6
  │
  ├─ 5 ─ 1                      : 5 + 1
  │
  ├─ 4 ┬ 2                      : 4 + 2
  │     │
  │     └ 1 ─ 1                 : 4 + 1 + 1
  │
  ├─ 3 ┬ 3                      : 3 + 3
  │     │
  │     ├ 2 ─ 1                 : 3 + 2 + 1
  │     │
  │     └ 1 ─ 1 ─ 1            : 3 + 1 + 1 + 1
  │
  ├─ 2 ┬ 2 ┬ 2                 : 2 + 2 + 2
  │     │   │
  │     │   └ 1 ─ 1            : 2 + 2 + 1 + 1
  │     │
  │     └ 1 ─ 1 ─ 1 ─ 1       : 2 + 1 + 1 + 1 + 1
  │
  └─ 1 ─ 1 ─ 1 ─ 1 ─ 1 ─ 1  : 1 + 1 + 1 + 1 + 1 + 1

                    図 : 整数 6 の分割

6 の場合、分割の仕方は上図のように 11 通りあります。この数を「分割数」といいます。分割の仕方を列挙する場合、整数 n から k 以下の整数を選んでいくと考えてください。まず、6 から 6 を選びます。すると、残りは 0 になるので、これ以上整数を分割することはできません。次に、6 から 5 を選びます。残りは 1 になるので、1 を選ぶしか方法はありません。

次に、4 を選びます。残りは 2 になるので、2 から 2 以下の整数を分割する方法になります。2 から 2 を選ぶと残りは 0 になるので 2 が得られます。1 を選ぶと残りは 1 になるので、1 + 1 が得られます。したがって、4 + 2, 4 + 1 + 1 となります。同様に、6 から 3 を選ぶと、残りは 3 から 3 以下の整数を選ぶ方法になります。

6 から 2 以下の整数を選ぶ方法は、残り 4 から 2 以下の整数を選ぶ方法になり、そこで 2 を選ぶと 2 から 2 以下の整数を選ぶ方法になります。1 を選ぶと 4 から 1 以下の整数を選ぶ方法になりますが、これは 1 通りしかありません。最後に 6 から 1 を選びますが、これも 1 通りしかありません。これらをすべて足し合わせると 11 通りになります。

それでは問題です。

分割数を求めるプログラムを作ってください。
分割の仕方をすべて求めるプログラムを作ってください。
整数を奇数のみで分割する場合を考えます。たとえば、6 を分割すると、5 + 1, 3 + 3, 3 + 1 + 1 + 1, 1 + 1 + 1 + 1 + 1 + 1 の 4 通りになります。この制限で分割数と分割の仕方を求めるプログラムを作ってください。

本ページでは Python3 (ver 3.8.10) を使ってプログラムを作ることにします。

●解答１

整数 n を k 以下の整数で分割する総数を求める関数を \(p(n, k)\) とすると、\(p(n, k)\) は次のように定義することができます。

\(\begin{eqnarray} p(n, k) = \begin{cases} 0 & \mathrm{if} \ n \lt 0 \ \mathrm{or} \ k \lt 1 \\ 1 & \mathrm{if} \ n = 0 \ \mathrm{or} \ k = 1 \\ p(n - k, k) + p(n, k - 1) & \mathrm{others} \end{cases} \end{eqnarray}\)

たとえば、p(6, 6) は次のように計算することができます。

p(6, 6) => p(0, 6) + p(6, 5)
        => 1 + p(1, 5) + p(6, 4)
        => 1 +    1    + p(2, 4) + p(6, 3)
        => 1 + 1 + 2 + 7
        => 11

p(2, 4) => p(-2, 4) + p(2, 3)
        =>    0     + p(-1, 3) + p(2, 2)
        =>    0     +    0     + p(0, 2) + p(2, 1)
        => 0 + 0 + 1 + 1
        => 2

p(6, 3) => p(3, 3) + p(6, 2)
        => p(0, 3) + p(3, 2) + p(4, 2) + p(6, 1)
        =>    1    + p(1, 2) + p(3, 1) + p(2, 2) + p(4, 1) + 1
        =>    1    +    1    +    1    + p(0, 2) + p(2, 1) + 1 + 1
        => 1 + 1 + 1 + 1 + 1 + 1 + 1
        => 7

分割数を求める関数を partition_number() とすると、関数 p(n, k) を使って次のようにプログラムすることができます。

リスト : 分割数 (partnum.py)

def part_num(n, k):
    if n == 1 or k == 1:
        return 1
    if n < 0 or k < 1:
        return 0
    else:
        return part_num(n - k, k) + part_num(n, k - 1)

def partition_number(n):
    return part_num(n, n)

def test(func, n):
    s = time.time()
    print(func(n))
    print(time.time() - s)

>>> test(partition_number, 40)
37338
0.014701366424560547
>>> test(partition_number, 60)
966467
0.3988070487976074
>>> test(partition_number, 80)
15796476
5.382274389266968

実行環境 : Ubuntu 20.04 LTS (WSL1), Intel Core i5-6200U 2.30GHz

関数 part_num() は p(n, k) の定義をそのままプログラムしただけです。ただし、このプログラムは二重再帰で何度も同じ値を求めているため実行速度はとても遅くなります。関数 part_num() をメモ化すると高速化することができますが、大きな値を計算すると Python のスタックがオーバーフローしてしまいます。

●動的計画法による高速化

動的計画法を使うと、大きな値でも高速に計算することができます。次の図を見てください。

k 
1 : [1,  1,  1,  1,  1,  1,  1] 

2 : [1,  1,  1+1=2, 1+1=2, 2+1=3, 2+1=3, 3+1=4]
 => [1,  1,  2,  2,  3,  3,  4]

3:  [1,  1,  2,  1+2=3, 1+3=4, 2+3=5, 3+4=7]
 => [1,  1,  2,  3,  4,  5,  7]

4:  [1,  1,  2,  3,  1+4=4, 1+5=6, 2+7=9]
 => [1,  1,  2,  3,  5,  6,  9

5:  [1,  1,  2,  3,  5,  1+6=7, 1+9=10]
 => [1,  1,  2,  3,  5,  7,  10]

6:  [1,  1,  2,  3,  5,  7,  10+1=11]
 => [1,  1,  2,  3,  5,  7,  11]

大きさ n + 1 のベクタを用意します。ベクタの添字が n を表していて、p(n, 1) から順番に値を求めていきます。p(n, 1) の値は 1 ですから、ベクタの要素は 1 に初期化します。次に、p(n, 2) の値を求めます。定義により p(n, 2) = p(n - 2, 2) + p(n, 1) なので、2 番目以降の要素に n - 2 番目の要素を加算すれば求めることができます。あとは、k の値をひとつずつ増やして同様の計算を行えば p(n, n) の値を求めることができます。

プログラムは次のようになります。

リスト : 分割数 (動的計画法)

def partition_number2(n):
    table = [1] * (n + 1)
    for k in range(2, n + 1):
        for m in range(k, n + 1):
            table[m] += table[m - k]
    return table[n]

説明をそのままプログラムしただけなので、とくに難しいところはないと思います。

それでは実際に試してみましょう。1000, 2000, 4000 の分割数を求めてみました。

>>> test(partition_number2, 1000)
24061467864032622473692149727991
0.06793832778930664
>>> test(partition_number2, 2000)
4720819175619413888601432406799959512200344166
0.272350549697876
>>> test(partition_number2, 4000)
1024150064776551375119256307915896842122498030313150910234889093895
0.9898276329040527

動的計画法の効果はとても高いですね。

●さらなる高速化

ところで、数がもっと大きくなると動的計画法を使ったプログラムでも遅くなります。実際に 5000, 6000, 7000 の分割数を求めてみましょう。

>>> test(partition_number2, 5000)
169820168825442121851975101689306431361757683049829233322203824652329144349
1.5292284488677979
>>> test(partition_number2, 6000)
467172753197020909297102464397369064336462915327003703385660552892507240534
9246129
2.2825496196746826
>>> test(partition_number2, 7000)
328569308034406157862809256359241668619501515745322406596990321574322363943
74450791229199
3.158601760864258

参考 URL 『分割数 - Wikipedia』によると、次の漸化式を使うと分割数を高速に求めることができるそうです。

\( p(k) = p(k - 1) + p(k - 2) - p(k - 5) - p(k - 7) + p(k - 12) + p(k - 15) - p(k - 22) - \cdots \)

漸化式の説明を参考 URL 『分割数 - Wikipedia』より引用します。

『ここで p(0) = 1 および負の整数 k に対して p(k) = 0 とし、和は (1/2)n(3n - 1) の形（ただし n は正または負の整数全体を走る）の一般五角数全体にわたってとるものとする（順に n = 1, -1, 2, -2, 3, -3, 4, -4 ..., とすると、値として 1, 2, 5, 7, 12, 15, 22, 26, 35, 40, 51, ... が得られる）。和における符号は交互に +, +, -, -, +, +, ... と続く。』

分割数 p(k) は k - 1 以下の分割数がわかれば求めることができます。この漸化式も動的計画法を使えば簡単にプログラムできます。次のリストを見てください。

リスト : 分割数 (オイラーの五角数定理)

# 五角数
def pentagon(n): return n * (3 * n - 1) // 2

def partition_number3(n):
    p = [0] * (n + 1)
    p[0] = 1
    for i in range(1, n + 1):
        j = 1
        s = 1
        while True:
            k = pentagon(j)
            if i < k: break
            p[i] += p[i - k] * s
            k = pentagon(-j)
            if i < k: break
            p[i] += p[i - k] * s
            j += 1
            s *= -1
    return p[n]

リスト p は分割数 p(k) を記憶するために使います。p[0] を 1 に初期化したあと、for ループで 1 から n までの分割数を順番に求めていきます。あとは、漸化式をそのままプログラムするだけです。変数 s は符号 (+. -) を表していて、j が奇数のとき s は 1 になり、j が偶数のときは -1 になります。

それでは実際に 5000, 10000, 20000 の分割数を求めてみましょう。

>>> test(partition_number3, 5000)
169820168825442121851975101689306431361757683049829233322203824652329144349
0.15317678451538086
>>> test(partition_number3, 10000)
3616725132563629398882047189095369549501603033931565042208186860588795256875406
6420592310556052906916435144
0.4225122928619385
>>> test(partition_number3, 20000)
2521148138125296979166195332304704522813289496018115934368503141080342844238015
64956623970731689824369192324789351994903016411826230578166735959242113097
1.2280046939849854

20000 の分割数を 1 秒ちょっとで求めることができました。PyPy3 を使うともっと速くなります。

$ pypy3
Python 3.6.9 (7.3.1+dfsg-4, Apr 22 2020, 05:15:29)
[PyPy 7.3.1 with GCC 9.3.0] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>>> from partnum import *
>>>> test(partition_number3, 5000)
169820168825442121851975101689306431361757683049829233322203824652329144349
0.14496374130249023
>>>> test(partition_number3, 5000)
169820168825442121851975101689306431361757683049829233322203824652329144349
0.026314973831176758
>>>> test(partition_number3, 5000)
169820168825442121851975101689306431361757683049829233322203824652329144349
0.02671670913696289
>>>> test(partition_number3, 10000)
3616725132563629398882047189095369549501603033931565042208186860588795256875406
6420592310556052906916435144
0.07793188095092773
>>>> test(partition_number3, 20000)
2521148138125296979166195332304704522813289496018115934368503141080342844238015
64956623970731689824369192324789351994903016411826230578166735959242113097
0.26520705223083496

20000 の分割数を 1 秒かからずに求めることができました。

●解答２

リスト : 整数の分割

def part_int_sub(n, k, a):
    if n == 0:   print(a)
    elif n == 1: print(a + [1])
    elif k == 1: print(a + [1] * n)
    else:
        if n >= k:
            part_int_sub(n - k, k, a + [k])
        part_int_sub(n, k - 1, a)

def partition_of_int(n): part_int_sub(n, n, [])

基本的な考え方は partition_number() と同じです。関数 part_int_sub() は選んだ数値を累積変数 a のリストに格納していくだけです。n が 0 の場合は a を出力し、n が 1 の場合は a に [1] を追加してから出力します。k が 1 の場合は [1] * n で要素が 1 で長さが n のリストを作成し、それを a と連結してから出力します。

5, 6, 7 の分割の仕方は次のようになります。

>>> partition_of_int(5)
[5]
[4, 1]
[3, 2]
[3, 1, 1]
[2, 2, 1]
[2, 1, 1, 1]
[1, 1, 1, 1, 1]
>>> partition_of_int(6)
[6]
[5, 1]
[4, 2]
[4, 1, 1]
[3, 3]
[3, 2, 1]
[3, 1, 1, 1]
[2, 2, 2]
[2, 2, 1, 1]
[2, 1, 1, 1, 1]
[1, 1, 1, 1, 1, 1]
>>> partition_of_int(7)
[7]
[6, 1]
[5, 2]
[5, 1, 1]
[4, 3]
[4, 2, 1]
[4, 1, 1, 1]
[3, 3, 1]
[3, 2, 2]
[3, 2, 1, 1]
[3, 1, 1, 1, 1]
[2, 2, 2, 1]
[2, 2, 1, 1, 1]
[2, 1, 1, 1, 1, 1]
[1, 1, 1, 1, 1, 1, 1]

●解答３

整数を奇数で分割する場合、関数 p(n, k) の k を奇数に限定するだけで求めることができます。漸化式は次のようになります。

\(\begin{eqnarray} p(n, k) = \begin{cases} 0 & \mathrm{if} \ n \lt 0 \ \mathrm{or} \ k \lt 1 \\ 1 & \mathrm{if} \ n = 0 \ \mathrm{or} \ k = 1 \\ p(n - k, k) + p(n, k - 2) & \mathrm{if} \ k \equiv 1 \pmod 2 \end{cases} \end{eqnarray}\)

k は奇数なので、p(n - k, k) + p(n, k - 1) の k - 1 を k - 2 に変更するだけで OK です。これを単純にプログラムすると次のようになります。

リスト : 奇数で分割したときの分割数

def part_num_odd(n, k):
    if n == 1 or k == 1:
        return 1
    if n < 0 or k < 1:
        return 0
    else:
        return part_num_odd(n - k, k) + part_num_odd(n, k - 2)

def partition_number_odd(n):
    if n % 2 == 0:
        return part_num_odd(n, n - 1)
    return part_num_odd(n, n)

>>> for x in range(1, 21): print(partition_number_odd(x))
...
1
1
2
2
3
4
5
6
8
10
12
15
18
22
27
32
38
46
54
64
>>> partition_number_odd(30)
296
>>> partition_number_odd(40)
1113
>>> partition_number_odd(50)
3658

動的計画法を使ったプログラムは次のようになります。

リスト : 奇数で分割したときの分割数 (2)

def partition_number_odd2(n):
    table = [1] * (n + 1)
    for k in range(3, n + 1, 2):
        for m in range(k, n + 1):
            table[m] += table[m - k]
    return table[n]

変数 k を奇数に限定するだけです。それでは、1000, 2000, 4000 の分割数を実際に求めてみましょう。

>>> test(partition_number_odd2, 1000)
8635565795744155161506
0.03764915466308594
>>> test(partition_number_odd2, 2000)
106972734349914451123354464808960
0.13586044311523438
>>> test(partition_number_odd2, 4000)
24884290037681681235695209792703824727967596678
0.5035812854766846

5000 程度の数であれば、Python3 でも高速に求めることができます。ただし、これよりも大きな数の分割数を求めるには、さらなる高速化が必要になりますが、M.Hiroi の勉強不足で漸化式がわかりません。これは今後の研究課題にしたいと思います。

分割の仕方を列挙するプログラムも簡単です。次のリストを見てください。

リスト : 奇数で分割する

def part_int_odd_sub(n, k, a):
    if n == 0:   print(a)
    elif n == 1: print(a + [1])
    elif k == 1: print(a + [1] * n)
    else:
        if n >= k:
            part_int_odd_sub(n - k, k, a + [k])
        part_int_odd_sub(n, k - 2, a)

def partition_of_int_odd(n):
    if n % 2 == 0:
        part_int_odd_sub(n, n - 1, [])
    else:
        part_int_odd_sub(n, n, [])

関数 part_int_odd_sub() の引数 k を奇数に限定するだけです。6, 7, 8 の分割の仕方は次のようになります。

>>> partition_of_int_odd(6)
[5, 1]
[3, 3]
[3, 1, 1, 1]
[1, 1, 1, 1, 1, 1]
>>> partition_of_int_odd(7)
[7]
[5, 1, 1]
[3, 3, 1]
[3, 1, 1, 1, 1]
[1, 1, 1, 1, 1, 1, 1]
>>> partition_of_int_odd(8)
[7, 1]
[5, 3]
[5, 1, 1, 1]
[3, 3, 1, 1]
[3, 1, 1, 1, 1, 1]
[1, 1, 1, 1, 1, 1, 1, 1]

●参考 URL

分割数 - Wikipedia
自然数の分割 - Wikipedia
分割数 - 桃の天然水, (inamori さん)

●プログラムリスト

#
# partnum.py : 分割数
#
#              Copyright (C) 2015-2022 Makoto Hiroi
#
import time

def part_num(n, k):
    if n == 1 or k == 1:
        return 1
    if n < 0 or k < 1:
        return 0
    else:
        return part_num(n - k, k) + part_num(n, k - 1)

def partition_number(n):
    return part_num(n, n)

def partition_number2(n):
    table = [1] * (n + 1)
    for k in range(2, n + 1):
        for m in range(k, n + 1):
            table[m] += table[m - k]
    return table[n]

# 五角数
def pentagon(n): return n * (3 * n - 1) // 2

def partition_number3(n):
    p = [0] * (n + 1)
    p[0] = 1
    for i in range(1, n + 1):
        j = 1
        s = 1
        while True:
            k = pentagon(j)
            if i < k: break
            p[i] += p[i - k] * s
            k = pentagon(-j)
            if i < k: break
            p[i] += p[i - k] * s
            j += 1
            s *= -1
    return p[n]

def part_int_sub(n, k, a):
    if n == 0:   print(a)
    elif n == 1: print(a + [1])
    elif k == 1: print(a + [1] * n)
    else:
        if n >= k:
            part_int_sub(n - k, k, a + [k])
        part_int_sub(n, k - 1, a)

def partition_of_int(n): part_int_sub(n, n, [])

def part_num_odd(n, k):
    if n == 1 or k == 1:
        return 1
    if n < 0 or k < 1:
        return 0
    else:
        return part_num_odd(n - k, k) + part_num_odd(n, k - 2)

def partition_number_odd(n):
    if n % 2 == 0:
        return part_num_odd(n, n - 1)
    return part_num_odd(n, n)

def partition_number_odd2(n):
    table = [1] * (n + 1)
    for k in range(3, n + 1, 2):
        for m in range(k, n + 1):
            table[m] += table[m - k]
    return table[n]

def part_int_odd_sub(n, k, a):
    if n == 0:   print(a)
    elif n == 1: print(a + [1])
    elif k == 1: print(a + [1] * n)
    else:
        if n >= k:
            part_int_odd_sub(n - k, k, a + [k])
        part_int_odd_sub(n, k - 2, a)

def partition_of_int_odd(n):
    if n % 2 == 0:
        part_int_odd_sub(n, n - 1, [])
    else:
        part_int_odd_sub(n, n, [])

def test(func, n):
    s = time.time()
    print(func(n))
    print(time.time() - s)

Puzzle DE Programming

数字のパズル : 分割数