多元线性方程整数解的计数
本文基于 SICP 1.2 换零钱方式的统计。
问题描述
统计
S
=
a
1
x
1
+
a
2
x
2
+
.
.
.
+
a
n
x
n
(
a
>
0
)
S=a_1x_1+a_2x_2+...+a_nx_n(a>0)
S=a1x1+a2x2+...+anxn(a>0)
的非负整数解的个数,列出这些解。
枚举
静态地看,单个
x
i
x_i
xi的范围容易确定,即
x
i
∈
[
0
,
S
a
i
]
x_i\in[0,\frac{S}{a_i}]
xi∈[0,aiS],记
X
=
[
0
,
S
a
1
]
×
[
0
,
S
a
2
]
×
.
.
.
×
[
0
,
S
a
n
]
X=[0,\frac{S}{a_1}]\times[0,\frac{S}{a_2}]\times ...\times[0,\frac{S}{a_n}]
X=[0,a1S]×[0,a2S]×...×[0,anS]
这就是解空间的可能范围,枚举法代码如下:
from itertools import product
def cc(s, a):
def p(x):
return sum([x * a for x, a in zip(x, a)])
return [x for x in product(*[range(s // ai + 1) for ai in a]) if p(x) == s]
枚举法是否做了一些无用功呢?答案是肯定的。
递归
动态地看,已决定的 x i x_i xi对未决定的 x i x_i xi的范围有影响。例如,在已经决定 x 1 x_1 x1的情况下, x 2 x_2 x2的上限从 S a 2 \frac{S}{a_2} a2S变为 S − a 1 x 1 a 2 \frac{S-a_1x_1}{a_2} a2S−a1x1。
问题的规模分布在两个维度上:
- 变量的个数
- 总和 S S S
一个成功的递归算法势必要以某种方式缩减需要考虑的变量个数以及总和 S S S。SICP 提供了一种分解思路,可以将方程的解分为两组:
- x n = 0 x_n=0 xn=0,此时相当于解 S = ∑ i = 1 n − 1 a i x i S=\sum\limits_{i=1}^{n-1}a_ix_i S=i=1∑n−1aixi
- x n ≥ 1 x_n\geq1 xn≥1,此时相当于解 S − a n = ∑ i = 1 n a i x i S-a_n=\sum\limits_{i=1}^{n}a_ix_i S−an=i=1∑naixi
这种划分是不重复、无遗漏的,且产生的子问题分别在变量个数与总和上比原问题小。分解不断进行下去就会触及边界,边界条件是:
- n = 0 n=0 n=0,不是解
- S < 0 S<0 S<0,不是解
- S = 0 S=0 S=0,是解
照此思路的 python 实现如下:
def solve(s, a):
def recursion(s, a):
if len(a) == 0:
return
elif s < 0:
return
elif s == 0:
xs.append(tuple(x))
return
else:
recursion(s, a[:-1]) # 减小变量个数
x[len(a) - 1] += 1
recursion(s - a[-1], a) # 减小S
x[len(a) - 1] -= 1
return
x = [0] * len(a)
xs = []
recursion(s, a)
return xs
对比
枚举空间的大小是各个变量范围的积,递归空间的大小可以通过在代码中设置计数变量来统计。将递归过程叶子结点的个数与枚举法的 X X X中元素个数做对比,发现递归空间的大小只有枚举空间大小的 1 40 \frac{1}{40} 401。测试运行时间:
%timeit cc(100, [1, 5, 10, 25, 50])
> 574 ms ± 47.9 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
%timeit solve(100, [1, 5, 10, 25, 50])
> 8.1 ms ± 616 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
574 / 8.1
> 71
上述实验显示,枚举耗时为递归的40倍,这个结果与搜索空间大小的比例吻合。
附录
jupyter notebook
完整实验代码,包括测量搜索空间大小的部分:
from itertools import *
l = 0
def cc(s, a):
global l
def p(x):
return sum([x * a for x, a in zip(x, a)])
l = len([x for x in product(*[range(s // ai + 1) for ai in a])])
return [x for x in product(*[range(s // ai + 1) for ai in a]) if p(x) == s]
x1 = cc(100, [1, 5, 10, 25, 50])
cnt = 0
def cr(s, a):
def r(s, a, x):
global cnt
if s == 0:
xs.append(tuple(x))
cnt += 1
elif s > 0 and len(a) != 0:
x = list(x)
r(s, a[:-1], x)
x = list(x)
x[len(a) - 1] += 1
r(s - a[-1], a, x)
else:
cnt += 1
xs = []
r(s, a, [0] * len(a))
return xs
x2 = cr(100, [1, 5, 10, 25, 50])
l / cnt
%timeit cc(100, [1, 5, 10, 25, 50])
%timeit cr(100, [1, 5, 10, 25, 50])
def solve(s, a):
def recursion(s, a):
if len(a) == 0:
return
elif s < 0:
return
elif s == 0:
xs.append(tuple(x))
return
else:
recursion(s, a[:-1]) # 减小变量个数
x[len(a) - 1] += 1
recursion(s - a[-1], a) # 减小S
x[len(a) - 1] -= 1
return
x = [0] * len(a)
xs = []
recursion(s, a)
return xs
x3 = solve(100, [1, 5, 10, 25, 50])
x3 == x2
set(x1) == set(x2)
%timeit solve(100, [1, 5, 10, 25, 50])
574 / 8.1
xs = []
recursion(s, a)
return xs
x3 = solve(100, [1, 5, 10, 25, 50])
x3 == x2
set(x1) == set(x2)
%timeit solve(100, [1, 5, 10, 25, 50])
574 / 8.1