はじめに
「部分和問題 (subset-sum problem)」は、要素が数値の集合 S において、要素の総和が M となる部分集合があるか判定する問題です。たとえば、集合 {2, 3, 5, 8} の場合、総和が 10 となる部分集合は {2, 3, 5} と {2, 8} がありますが、14 となる部分集合はありません。部分集合の総数は、要素数を N とすると 2N 個になるので、n が大きくなるとナイーブな方法では時間がかかってしまいます。実際には、動的計画法などの手法を使うことで、現実的な時間で部分和問題を解くことができます。
●べき集合
最初にナイーブな方法で部分和問題を解いてみましょう。今回は要素を正整数に限定します。部分和問題は「べき集合」を生成する高階関数 power_set() を作ると簡単です。次のリストを見てください。
リスト : べき集合 (subsetsum.py)
def power_set(func, xs):
def iter(i, a):
if i == len(xs):
func(a)
else:
iter(i + 1, a)
a.append(xs[i])
iter(i + 1, a)
a.pop()
#
iter(0, [])
引数 func が関数で、xs が集合を表すリストです。実際の処理は局所関数 iter() で行います。引数 i が添字で、a が累積変数 (リスト) です。処理内容は簡単で、i が len(xs) と等しい場合、部分集合がひとつできたので func(a) を評価します。あとは、i 番目の要素を選ばない場合は iter() を再帰呼び出しし、選ぶ場合は a に xs[i] を追加してから iter() を再帰呼び出しします。戻ってきたら追加した要素を削除します。これですべての部分集合を生成することができます。
簡単な実行例を示します。
>>> from subsetsum import * >>> power_set(print, [1,2,3,4]) [] [4] [3] [3, 4] [2] [2, 4] [2, 3] [2, 3, 4] [1] [1, 4] [1, 3] [1, 3, 4] [1, 2] [1, 2, 4] [1, 2, 3] [1, 2, 3, 4]
部分集合は空集合 [ ] を含めて 16 通りあります。
●ナイーブな解法
この power_set() を使うと部分和問題のプログラムは次のようになります。
リスト : 部分和問題 (subsetsum.py)
def subset_sum0(n, xs):
def check(ys):
if sum(ys) == n: print(ys)
power_set(check, xs)
部分集合 ys の総和を sum で求め、n と等しい場合は print(で出力します。それでは実行してみましょう。
>>> subset_sum0(10, [2,3,5,8]) [2, 8] [2, 3, 5] >>> subset_sum0(14, [2,3,5,8]) >>>
とても簡単ですね。ただし、集合の要素数が多くなると、実行時間がかかるようになります。次のテストプログラムを見てください。
リスト : 簡単なテスト
nums = [ 1, 2, 3, 5, 8, 13, 21, 34, 55, 89,
144, 233, 377, 610, 987, 1597, 2584, 4181, 6765, 10946,
17711, 28657, 46368, 75025, 121393]
def test(func):
for x in range(20, 26):
xs = nums[0:x]
s = time.time()
print(func(sum(xs) - 1, xs))
print(time.time() - s)
配列 nums はフィボナッチ数列になっています。要素の総和を M とすると、1 から M までの整数は、要素を組み合わせて必ず作ることができます。これはフィボナッチ数列の面白い特徴です。
テストは数列の長さを 20 から一つずつ増やしながら、総和 - 1 となる組み合わせを subset_sum0() で求め、その実行時間を計測します。結果は次のようになりました。
>>> from subsetsum import * >>> test(subset_sum0) [2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765, 10946] None 0.6480607986450195 [2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765, 10946, 17711] None 1.3853185176849365 [2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765, 10946, 17711, 28657] None 2.5871338844299316 [2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765, 10946, 17711, 28657, 46368] None 5.3003764152526855 [2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765, 10946, 17711, 28657, 46368, 75025] None 11.62532901763916 [2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765, 10946, 17711, 28657, 46368, 75025, 121393] None 21.400745630264282 実行環境 : Python 3.8.10, Ubuntu 20.04 LTS (WSL1), Intel Core i5-6200U 2.30GHz
要素がひとつ増えると実行時間は約 2 倍になっていることがわかります。要素数を N とすると、subset_sum0() の実行時間は 2N に比例する遅いプログラムなのです。
●深さ優先探索＋枝刈り
次は深さ優先探索に枝刈りを加えた解法プログラムを作ってみましょう。今回の部分和問題は要素を正整数値に限定しているので、二種類の枝刈りを考えることができます。ひとつは部分集合の総和が求める値 N を超えた場合です。残りの要素は正整数なので、これ以上要素を追加しても解を得られないのは明白ですね。もうひとつは、部分集合の総和に残りの要素をすべて足しても N に満たない場合です。これも解を得られないのは明白です。
部分集合の総和を S, 残りの要素の総和を R, 求める値を N とすると、探索を行う条件は次の式で表すことができます。
S < N and S + R >= N
これをプログラムすると次のようになります。
リスト : 部分和問題 (dfs + cut)
def subset_sum1(n, xs):
def iter(i, a, s, r):
if s == n:
print(a)
elif s < n and s + r >= n:
iter(i + 1, a, s, r - xs[i])
a.append(xs[i])
iter(i + 1, a, s + xs[i], r - xs[i])
a.pop()
#
iter(0, [], 0, sum(xs))
実際の処理は局所関数 iter() で行います。引数 s は求めている部分集合の総和、r は残りの要素の総和を表します。s < n かつ s + r >= n ならば条件を満たすので、iter() を再帰呼び出しして探索を続行します。とても簡単ですね。
それでは試してみましょう。
>>> test(subset_sum1) [2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765, 10946] None 0.0007145404815673828 [2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765, 10946, 17711] None 0.0008144378662109375 [2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765, 10946, 17711, 28657] None 0.0007948875427246094 [2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765, 10946, 17711, 28657, 46368] None 0.0005850791931152344 [2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765, 10946, 17711, 28657, 46368, 75025] None 0.0007252693176269531 [2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765, 10946, 17711, 28657, 46368, 75025, 121393] None 0.00042247772216796875
とても速くなりましたね。ただし、これは枝刈りがうまくいった場合であり、データによっては枝刈りが機能しない場合もありえます。たとえば、xs の最後尾の要素から 1 を引いた値 (xs[-1] - 1) を判定してみましょう。結果は次のようになりました。
>>> test1(subset_sum1) [1, 3, 8, 21, 55, 144, 377, 987, 2584, 6765] None 0.5068469047546387 [2, 5, 13, 34, 89, 233, 610, 1597, 4181, 10946] None 1.01094388961792 [1, 3, 8, 21, 55, 144, 377, 987, 2584, 6765, 17711] None 1.967355489730835 [2, 5, 13, 34, 89, 233, 610, 1597, 4181, 10946, 28657] None 4.194068431854248 [1, 3, 8, 21, 55, 144, 377, 987, 2584, 6765, 17711, 46368] None 8.148783683776855 [2, 5, 13, 34, 89, 233, 610, 1597, 4181, 10946, 28657, 75025] None 15.926688432693481
実行時間はナイーブな方法より速くなりましたが、かなり遅いですね。ところが xs を降順に並べておくと、実行時間はとても速くなります。
>>> test1(subset_sum1) [6765, 2584, 987, 377, 144, 55, 21, 8, 3, 1] None 0.0004088878631591797 [10946, 4181, 1597, 610, 233, 89, 34, 13, 5, 2] None 0.00023365020751953125 [17711, 6765, 2584, 987, 377, 144, 55, 21, 8, 3, 1] None 0.00042176246643066406 [28657, 10946, 4181, 1597, 610, 233, 89, 34, 13, 5, 2] None 0.0004284381866455078 [46368, 17711, 6765, 2584, 987, 377, 144, 55, 21, 8, 3, 1] None 0.0004036426544189453 [75025, 28657, 10946, 4181, 1597, 610, 233, 89, 34, 13, 5, 2] None 0.00039124488830566406
ちなみに、関数 random.shuffle() でシャッフルしてみたところ、実行時間は次のようになりました。
>>> test1(subset_sum1) [21, 6765, 144, 987, 1, 2584, 55, 8, 377, 3] None 0.06557774543762207 [13, 89, 10946, 2, 1597, 34, 233, 5, 4181, 610] None 0.02338862419128418 [144, 987, 1, 2584, 21, 17711, 3, 6765, 377, 55, 8] None 0.05066370964050293 [1597, 5, 233, 34, 28657, 610, 13, 89, 2, 4181, 10946] None 0.6216027736663818 [2584, 6765, 21, 377, 3, 46368, 144, 55, 8, 987, 17711, 1] None 1.3327007293701172 [13, 2, 233, 34, 1597, 610, 28657, 75025, 5, 10946, 4181, 89] None 3.507883071899414
今回のプログラムと問題では、大きな値から試した方が枝刈りの効果は高いようです。興味のある方はいろいろ試してみてください。
●動的計画法
次は「動的計画法」で部分和問題を解いてみましょう。動的計画法の説明は拙作のページ 動的計画法 をお読みください。部分和問題の場合、総和が N となる部分集合があるか判定するだけでよければ、動的計画法で解くことが可能です。
動的計画法で部分和問題を解く場合、要素をひとつずつ追加しながら、総和 N となる部分集合があるか判定します。簡単な例を示しましょう。次の図を見てください。
上図は xs = {2, 3, 5, 8} で N = 10 の部分集合があるか判定する場合です。最初に N + 1 の配列 Ai を用意します。空集合の総和は 0 なので A0[0] に○をセットします。次に、要素 2 を追加します。部分集合は { } と {2} になります。A1[0] と A1[2] に○をセットします。その次に要素 3 を追加します。追加される部分集合は {3} と {2, 3} になるので、A2[0], A2[2], A2[3] と A2[5] に○をセットします。
つまり、i 番目の要素 x を追加する場合、Ai-1 で○が付いている位置を y とすると、Ai[y] と Ai[x + y] に○をセットすればいいわけです。添字 y は部分集合の総和を表しています。Ai[y] に○をセットすることは、その部分集合に x を加えないことを意味し、Ai[x + y] に○をセットすることは、その部分集合に x を追加することを意味するわけです。
次に 5 を追加します。A2 の○の位置は 0, 2, 3, 5 なので、これに 5 を足した 5, 7, 8, 10 の位置に○をセットします。最後に 8 を追加します。A3 の○の位置は 0, 2, 3, 5, 7, 8, 10 なので、これに 8 を足した 8, 10 の位置に○をセットします。A4[10] の値が○になので、部分和が 10 となる部分集合があることがわかります。
もうひとつ簡単な例を示しましょう。今度は総和が 14 となる部分集合があるか判定します。
3 番目で○の位置は 0, 2, 3, 5, 7, 8, 10 です。次は 8 を追加しますが、総和 14 より大きい値は不要なので、8, 10, 11, 13 の位置に○を追加します。14 の位置は×なので、総和が 14 となる部分集合は無いことがわかります。
プログラムは次のようになります。
リスト : 部分和問題 (動的計画法)
def subset_sum2(n, xs):
a = [False] * (n + 1)
a[0] = True
for x in xs:
for i in range(n - x, -1, -1):
if a[i]: a[x + i] = True
return a[n]
○を True で、×を False で表します。配列をひとつで済ますため、配列の後ろから True の位置を検索していることに注意してください。また、検索の開始位置を n - x とすることで、True をセットするときの範囲チェックを省略しています。今回のプログラムでは xs の要素をすべてチェックしていますが、x + i が n と等しくなったら return で True を返してもかまいません。あとは特に難しいところはないと思います。
それでは実際に試してみましょう。
>>> test(subset_sum2) True 0.02694845199584961 True 0.041280508041381836 True 0.08095121383666992 True 0.10300326347351074 True 0.1945786476135254 True 0.30707287788391113
ナイーブな方法よりも高速になりましたが、dfs + cut にはかないませんでした。集合の要素数を M, 総和を N とすると、今回のプログラムの実行速度は N * M に比例します。たとえば、N の値を (xs[-1] - 1) とすると、実行結果は次のようになります。
>>> test1(subset_sum2) True 0.009999275207519531 True 0.016277790069580078 True 0.023107051849365234 True 0.05161857604980469 True 0.07762503623962402 True 0.12096428871154785
N が小さくなったので、実行時間も速くなりました。このように、動的計画法では N が大きくなると、どうしても時間がかかるようになります。そこで、次は配列を使わずにプログラムを作ってみましょう。
●プログラムの改良
今までのプログラムは総和 N に対応する配列を用意しました。N が大きくなると、配列から True を検索する処理に時間がかかるようになります。そこで、配列のかわりにセット (set) を使って部分集合の総和を管理することにします。プログラムは次のようになります。
リスト : 部分和問題 (動的計画法 set 版)
def subset_sum3(n, xs):
a = {0}
for x in xs:
b = []
for y in a:
if x + y <= n: b.append(x + y)
a.update(b)
return n in a
最初に部分集合の総和を管理するセット a を {0} に初期化します。セット a から要素を取り出している間、セット a は更新できないので、新しく追加する総和の値を配列 b にセットしておきます。a の要素をすべてチェックしてから a.update(b) で新しい総和を a に追加します。このとき、重複要素のチェックが行われるので、同じ値が複数追加されることはありません。
それでは実行結果を示します。
>>> test(subset_sum3) True 0.013036012649536133 True 0.018460512161254883 True 0.023808956146240234 True 0.04677009582519531 True 0.08202719688415527 True 0.14128828048706055
subset_sum2() よりも速くなりましたが、それでも dfs + cut には及びません。そこで、dfs + cut と同様の枝刈りを行うことにします。次のリストを見てください。
リスト : 部分和問題 (動的計画法 set + cut 版)
def subset_sum31(n, xs):
a = {0}
r = sum(xs)
for x in xs:
b = []
r -= x
for y in a:
if x + y <= n and x + y + r >= n: b.append(x + y)
a.update(b)
return n in a
残りの要素の総和を変数 r で管理します。x + y + r が n 未満であれば、総和 x + y は n にならないのは明白です。x + y を集合 a を追加する必要はありません。
実行結果は次のようになりました。
>>> test(subset_sum31) True 0.0003428459167480469 True 0.00048041343688964844 True 0.0004889965057373047 True 0.0005033016204833984 True 0.0004897117614746094 True 0.0005037784576416016
dfs + cut と同等以上の実行速度になりました。
●簡単なテスト
もうひとつ、簡単なテストを行ってみましょう。次のリストを見てください。
リスト : 簡単なテスト (2)
nums2 = [93057, 89728, 88489, 84353, 83190, 73846, 63272, 63121, 62842, 60169,
54123, 51409, 46751, 39788, 35209, 30222, 29706, 24978, 14548, 10283]
def test2():
for x in [16,17,18,19,20]:
xs = nums2[0:x]
n = sum(xs) // 2
s = time.time()
subset_sum1(n, xs)
print(time.time() - s)
s = time.time()
print(subset_sum31(n, xs)
print(time.time() - s)
print("-----")
nums1 の整数値は乱数で生成して、それを降順に並べたものです。subset_sum1() と subset_sum31() で実行時間を計測したところ、結果は次のようになりました。
>>> test2() 0.008139848709106445 False 0.0045795440673828125 ----- [93057, 89728, 88489, 73846, 63272, 46751, 39788, 29706] 0.02464604377746582 True 0.008201122283935547 ----- 0.04703068733215332 False 0.010396718978881836 ----- [89728, 84353, 83190, 63272, 62842, 46751, 39788, 30222, 29706, 14548] [89728, 88489, 63272, 63121, 62842, 46751, 39788, 35209, 30222, 24978] 0.051488399505615234 True 0.02222299575805664 ----- [88489, 84353, 83190, 62842, 60169, 54123, 51409, 29706, 24978, 10283] [89728, 63272, 63121, 62842, 60169, 54123, 46751, 39788, 30222, 24978, 14548] [93057, 88489, 84353, 83190, 73846, 51409, 35209, 29706, 10283] [93057, 89728, 73846, 63272, 63121, 46751, 39788, 35209, 30222, 14548] 0.08587002754211426 True 0.0312504768371582 -----
dfs + cut よりも動的計画法の方が速くなりました。ただし、今回のプログラム (subset_sum31) は集合 xs の要素を増やすとメモリを大量に消費する場合があります。興味のある方は集合の要素数を増やすなど、いろいろなデータで試してみてください。
●擬似完全数と不思議数
自然数 n の約数において、n 以外の約数の集合 xs を考えます。xs の総和が n に等しい数を「完全数」といい、xs の部分和 (部分集合の要素の和) が n に等しい数を「擬似完全数」といいます。たとえば、6 の約数は {1, 2, 3, 6} で、6 以外の約数 {1, 2, 3} の総和は 6 になるので完全数です。12 の約数は {1, 2, 3, 4, 6, 12} ですが、12 以外の約数の総和は 16 になるので完全数ではありません。ところが、部分集合 {1, 2, 3, 6} の和は 12 になるので、12 は擬似完全数になります。なお、完全数は擬似完全数でもあります。
自然数 n の約数の総和が 2 * n より大きい数のことを「過剰数」といいます。擬似完全数は過剰数になります。過剰数だが擬似完全数ではない数のことを「不思議数」といいます。たとえば、70 の約数は {1, 2, 5, 7, 10, 14, 35, 70} ですが、70 を除いた約数の集合の部分和は 70 にならないので不思議数になります。
それでは問題です。n <= 10000 の範囲で、擬似完全数の個数とすべての不思議数を求めてみましょう。約数の求め方は拙作のページ Puzzle DE Programming 約数 で説明しています。そのとき作成した関数 divisor() を使うと、プログラムは次のようになります。
リスト : 擬似完全数と不思議数
# 擬似完全数
def semiperfect_number(n):
count = 0
for x in range(2, n + 1):
xs = divisor(x)
if sum(xs) >= 2 * x and subset_sum31(x, xs[:-1]):
count += 1
return count
# 不思議数
def weird_number(n):
for x in range(2, n + 1):
xs = divisor(x)
if sum(xs) >= 2 * x and not subset_sum31(x, xs[:-1]):
print(x, end=' ')
print()
s = time.time()
print(semiperfect_number(10000))
print(time.time() - s)
s = time.time()
weird_number(10000)
print(time.time() -s)
2 つの関数に分けましたが、一つにまとめてもかまいません。どちらの関数も約数を求める関数 divisor() を呼び出して x の約数 xs を求めます。その総和が 2 * n 以上であれば、部分和問題を解く subset_sum31() を呼び出して、x と等しくなる部分集合があるかチェックします。このとき、xs から末尾の要素 (x の値) を削除することをお忘れなく。
実行結果は次のようになりました。
$ python3 divisor.py 2485 5.665281772613525 70 836 4030 5830 7192 7912 9272 5.583245277404785 $ pypy3 divisor.py 2485 1.6282453536987305 70 836 4030 5830 7192 7912 9272 1.5397768020629883 実行環境 : Ubuntu 20.04 LTS (WSL1), Intel Core i5-6200U 2.30GHz
擬似完全数の個数は 2485 個で、不思議数の個数は 7 個になりました。実行時間は Python3 で約 6 秒、divisor() だけならば 0.1 秒程度なので、部分和問題を求める subset_sum31() に時間がかかるようです。PyPy3 を使うと約 1.6 秒になりました。興味のある方はいろいろ試してみてください。
●プログラムリスト１
#
# subsetsum.py : 部分和問題
#
# Copyright (C) 2022 Makoto Hiroi
#
import random, time
def power_set(func, xs):
def iter(i, a):
if i == len(xs):
func(a)
else:
iter(i + 1, a)
a.append(xs[i])
iter(i + 1, a)
a.pop()
#
iter(0, [])
def subset_sum0(n, xs):
def check(ys):
if sum(ys) == n: print(ys)
power_set(check, xs)
def subset_sum1(n, xs):
def iter(i, a, s, r):
if s == n:
print(a)
elif s < n and s + r >= n:
iter(i + 1, a, s, r - xs[i])
a.append(xs[i])
iter(i + 1, a, s + xs[i], r - xs[i])
a.pop()
#
iter(0, [], 0, sum(xs))
def subset_sum2(n, xs):
a = [False] * (n + 1)
a[0] = True
for x in xs:
for i in range(n - x, -1, -1):
if a[i]: a[x + i] = True
return a[n]
def subset_sum3(n, xs):
a = set([0])
for x in xs:
b = []
for y in a:
if x + y <= n: b.append(x + y)
a.update(b)
return n in a
def subset_sum31(n, xs):
a = set([0])
r = sum(xs)
for x in xs:
b = []
r -= x
for y in a:
if x + y <= n and x + y + r >= n: b.append(x + y)
a.update(b)
return n in a
nums = [ 1, 2, 3, 5, 8, 13, 21, 34, 55, 89,
144, 233, 377, 610, 987, 1597, 2584, 4181, 6765, 10946,
17711, 28657, 46368, 75025, 121393]
def test(func):
for x in range(20, 26):
xs = nums[0:x]
s = time.time()
print(func(sum(xs) - 1, xs))
print(time.time() - s)
def test1(func):
for x in range(20, 26):
xs = nums[0:x]
n = xs[-1] -1
# xs.sort(reverse = True)
# random.shuffle(xs)
s = time.time()
print(func(n, xs))
print(time.time() - s)
nums2 = [93057, 89728, 88489, 84353, 83190, 73846, 63272, 63121, 62842, 60169,
54123, 51409, 46751, 39788, 35209, 30222, 29706, 24978, 14548, 10283]
def test2():
for x in [16,17,18,19,20]:
xs = nums2[0:x]
n = sum(xs) // 2
s = time.time()
subset_sum1(n, xs)
print(time.time() - s)
s = time.time()
print(subset_sum31(n, xs))
print(time.time() - s)
print("-----")
●プログラムリスト２
#
# divisor.py : 約数、擬似完全数、不思議数
#
# Copyright (C) 2022 Makoto Hiroi
#
import time
from subsetsum import *
#
# 素因数分解 (単純な試し割り)
#
def factor_sub(n, m):
c = 0
while n % m == 0:
c += 1
n //= m
return c, n
def factorization(n):
buff = []
c, m = factor_sub(n, 2)
if c > 0: buff.append((2, c))
c, m = factor_sub(m, 3)
if c > 0: buff.append((3, c))
x = 5
while m >= x * x:
c, m = factor_sub(m, x)
if c > 0: buff.append((x, c))
if x % 6 == 5:
x += 2
else:
x += 4
if m > 1: buff.append((m, 1))
return buff
#
# 約数を求める
#
# p^q の約数を求める
def divisor_sub(p, q):
a = []
for i in range(0, q + 1):
a.append(p ** i)
return a
def divisor(n):
xs = factorization(n)
ys = divisor_sub(xs[0][0], xs[0][1])
for p, q in xs[1:]:
ys = [x * y for x in divisor_sub(p, q) for y in ys]
return sorted(ys)
# 擬似完全数
def semiperfect_number(n):
count = 0
for x in range(2, n + 1):
xs = divisor(x)
if sum(xs) >= 2 * x and subset_sum31(x, xs[:-1]):
count += 1
return count
# 不思議数
def weird_number(n):
for x in range(2, n + 1):
xs = divisor(x)
if sum(xs) >= 2 * x and not subset_sum31(x, xs[:-1]):
print(x, end=' ')
print()
s = time.time()
print(semiperfect_number(10000))
print(time.time() - s)
s = time.time()
weird_number(10000)
print(time.time() -s)