线性方程的整数解

最新推荐文章于 2021-08-30 16:39:38 发布

__Lysias__

最新推荐文章于 2021-08-30 16:39:38 发布

阅读量1.2k

点赞数

分类专栏：离散数学文章标签： sicp python

本文链接：https://blog.csdn.net/weixin_42156883/article/details/104137309

版权

离散数学专栏收录该内容

9 篇文章 0 订阅

订阅专栏

多元线性方程整数解的计数

本文基于 SICP 1.2 换零钱方式的统计。

文章目录

多元线性方程整数解的计数

问题描述

统计
$S=a_1x_1+a_2x_2+...+a_nx_n(a>0)$
的非负整数解的个数，列出这些解。

枚举

静态地看，单个 $x_i$ 的范围容易确定，即 $x_i\in[0,\frac{S}{a_i}]$ ，记
$X=[0,\frac{S}{a_1}]\times[0,\frac{S}{a_2}]\times ...\times[0,\frac{S}{a_n}]$
这就是解空间的可能范围，枚举法代码如下：

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$ 的范围有影响。例如，在已经决定 $x_1$ 的情况下， $x_2$ 的上限从 $\frac{S}{a_2}$ 变为 $\frac{S-a_1x_1}{a_2}$ 。

问题的规模分布在两个维度上：

变量的个数
总和 $S$

一个成功的递归算法势必要以某种方式缩减需要考虑的变量个数以及总和 $S$ 。SICP 提供了一种分解思路，可以将方程的解分为两组：

$x_n=0$ ，此时相当于解 $S=\sum\limits_{i=1}^{n-1}a_ix_i$
$x_n\geq1$ ，此时相当于解 $S-a_n=\sum\limits_{i=1}^{n}a_ix_i$

这种划分是不重复、无遗漏的，且产生的子问题分别在变量个数与总和上比原问题小。分解不断进行下去就会触及边界，边界条件是：

$n = 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$ 中元素个数做对比，发现递归空间的大小只有枚举空间大小的 $\frac{1}{40}$ 。测试运行时间：

%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