小编典典

获取线性方程的所有正整数解

python

我玩的游戏有一个谜题,涉及解决以下方程式:

x*411 + y*295 + z*161 = 3200

不想认为我只是将其拍打成sympy,到那时我还没有真正用完:

>>> from sympy import *
>>> x, y, z = symbols('x y z', integer=True, positive=True)
>>> solve(x*411 + y*295 + z*161 - 3200, [x, y, z])
[{x: -295*y/411 - 161*z/411 + 3200/411}]

嗯,这只给了我一个从属的解决方案,但是我希望将变量约束到的域中所有可能的解决方案,例如(假设没有其他解决方案)[{x: 4, y: 2, z:6}][(4, 2, 6)]

当然,我现在可以在嵌套循环中手动替换两个变量,或者手动解决它(就像我在上面的解决方案中所做的那样),但是我想知道如何让sympy(或另一个库)为我做这件事。


阅读 347

收藏
2021-01-20

共1个答案

小编典典

SymPy可以求解Diophantine方程,但没有生成正解的内置方法。使用Sage可以轻松完成此任务:这是四行代码,可生成方程式的所有
非负整数解

p = MixedIntegerLinearProgram()
w = p.new_variable(integer=True, nonnegative=True)
p.add_constraint(411*w[0] + 295*w[1] + 161*w[2] == 3200)
p.polyhedron().integral_points()

输出是 ((4, 2, 6),)

在幕后,integral_points很可能只会运行多个循环;尽管当这似乎不起作用时,它会尝试使用Smith范式。

我知道您想要积极的解决方案,但是(a)很容易从答案中排除任何包含零的元组;(b)在求解之前,也很容易用x-1等替换x;(c)坚持“否定性”使得使用
上述混合整数线性编程模块可以轻松创建多面体。

根据文档,也可以直接从不等式(“
Hrep”)构建多面体对象。这将允许人们明确地说x>
= 1,依此类推,但是我在这条路线上还没有成功。

使用SymPy

SymPy的Diophantine模块的输出是一个参数化解决方案,例如

(t_0, 2627*t_0 + 161*t_1 - 19200, -4816*t_0 - 295*t_1 + 35200)

在您的示例中。可以将其循环使用,以非常有效的方式生成解决方案。症结在于寻找参数t_0和t_1的界限。由于这只是一个示例,我查看了上面的最后一个表达式,并将限制35200/4816和35200/295直接插入下面的循环中。

from sympy import *
x, y, z = symbols('x y z')
[s] = diophantine(x*411 + y*295 + z*161 - 3200)
print(s)
t_0, t_1 = s[2].free_symbols
for t0 in range(int(35200/4816)+1):
    for t1 in range(int(35200/295)+1):
        sol = [expr.subs({t_0: t0, t_1: t1}) for expr in s]
        if min(sol) > 0:
            print(sol)

输出为[4, 2, 6]

2021-01-20