找到两个线性等式成立的整数集
find the set of integers for which two linear equalities holds true
我可以使用什么算法来找到 n1, n2, ... ,n7 的所有正整数值的集合，满足以下不等式。
97n1 + 89n2 + 42n3 + 20n4 + 16n5 + 11n6 + 2n7 - 185 > 0
-98n1 - 90n2 - 43n3 - 21n4 - 17n5 - 12n6 - 3n7 + 205 > 0
n1 >= 0, n2 >= 0, n3 >=0. n4 >=0, n5 >=0, n6 >=0, n7 >= 0
例如一组 n1= 2, n2 = n3 = ... = n7 =0 使不等式成立。我如何找出所有其他值集？类似问题已发布在 M.SE.
ADDED:: 我需要概括 n 个变量（可能很大）的解决方案。我可以申请什么程序？对于 另一个 特殊情况 n=8
97n1 + 89n2 + 42n3 + 20n4 + 16n5 + 11n6 + 6n7 + 2n8 - 185 > 0
-98n1 - 90n2 - 43n3 - 21n4 - 17n5 - 12n6 - 7 - 3n8 + 205 > 0
n1 >= 0, n2 >= 0, n3 >=0. n4 >=0, n5 >=0, n6 >=0, n7 >= 0, n8 >= 0
Python 需要很长时间。 Wolfram Mathematica 显示在不到一分钟内有 4015 个解决方案。
Length[Solve[{97 n1 + 89 n2 + 42 n3 + 20 n4 + 16 n5 + 11 n6 + 6 n7 +
2 n8 - 185 > 0,
-98 n1 - 90 n2 - 43 n3 - 21 n4 - 17 n5 - 12 n6 - 7 n7 - 3 n8 +
205 > 0,
n1 >= 0, n2 >= 0, n3 >= 0, n4 >= 0, n5 >= 0, n6 >= 0, n7 >= 0,
n8 >= 0}, {n1, n2, n3, n4, n5, n6, n7, n8}, Integers]]
有 1013 个解决方案，但我不知道最有效的解决方法。
看第二个不等式，17 * n5不能大于205（否则整个左边不能为正）。这导致 n5 <= 12。通过为每个其他变量计算类似的界限，您可以将问题简化为可以使用嵌套循环快速解决的问题。
此 Java 代码打印出所有解决方案。
for (int n1 = 0; n1 <= 2; n1++) {
for (int n2 = 0; n2 <= 2; n2++) {
for (int n3 = 0; n3 <= 4; n3++) {
for (int n4 = 0; n4 <= 9; n4++) {
for (int n5 = 0; n5 <= 12; n5++) {
for (int n6 = 0; n6 <= 17; n6++) {
for (int n7 = 0; n7 <= 68; n7++) {
if (97 * n1 + 89 * n2 + 42 * n3 + 20 * n4 + 16 * n5 + 11 * n6 + 2 * n7 - 185 > 0
&& -98 * n1 - 90 * n2 - 43 * n3 - 21 * n4 - 17 * n5 - 12 * n6 - 3 * n7 + 205 > 0) {
System.out.println(Arrays.asList(n1, n2, n3, n4, n5, n6, n7));
}
}
}
}
}
}
}
}
只要
的前几项之和停止每个循环，就可以提高效率
98n1 + 90n2 + 43n3 + 21n4 + 17n5 + 12n6 + 3n7
达到 205。
例如，n4循环可以替换为
for (int n4 = 0; 98 * n1 + 90 * n2 + 43 * n3 + 21 * n4 < 205; n4++)
如果你把7个循环都改成这样，那么1013个解都可以非常快的找到。
好像是整数解的线性规划。我认为已经实施了一些算法。
尝试考虑 Linear Programming Tool/Libraries for Java.
在你给出的两个例子中，我注意到了相同的模式。例如，对于第一种情况，如果将两个方程相加，您将得到
-n1 - n2 - n3 - n4 - n5 - n6 - n7 + 20 > 0
可以重新排列为
n1 + n2 + n3 + n4 + n5 + n6 + n7 < 20
这是一个很好的有界方程，您可以对其进行暴力破解。具体来说，您可以从 0 到 19 迭代 n1，从 0 到 19-n1 迭代 n2，等等。一个可能的解决方案是 (0, 0, 0, 0, 0, 0 , 0), 但我们注意到这不满足我们原来的方程。因此，只需为 (n1, n2, ..., n7) 生成所有可能的值，并只保留满足您的等式的值。将所有这些结果硬编码为
def find_solutions(N):
sols = []
for n1 in xrange(N):
for n2 in xrange(N-n1):
for n3 in xrange(N-n1-n2):
for n4 in xrange(N-n1-n2-n3):
for n5 in xrange(N-n1-n2-n3-n4):
for n6 in xrange(N-n1-n2-n3-n4-n5):
for n7 in xrange(N-n1-n2-n3-n4-n5-n6):
if (97*n1 + 89*n2 + 42*n3 + 20*n4 + 16*n5 + 11*n6 + 2*n7 - 185 > 0 and
-98*n1 - 90*n2 - 43*n3 - 21*n4 - 17*n5 - 12*n6 - 3*n7 + 205 > 0):
sols.append((n1, n2, n3, n4, n5, n6, n7))
return sols
find_solutions(20) 在 0.6 秒内找到所有 1013 个解决方案。同样，对于第二种情况，它在 2.3 秒内找到所有 4015 个解决方案。现在，这根本不容易概括，但它表明，使用聪明的方法，Python 或任何其他语言，不必很慢。
另一方面，递归允许我们将其推广到任意数量的嵌套循环，但代价是 运行ning 有点慢。
def find_solutions(N, coeffs, depth=0, variables=None, subtotal=None, solutions=None):
if variables is None:
solutions = []
subtotal = [0 for _ in xrange(len(coeffs[0]))]
variables = [0 for _ in xrange(len(coeffs[0])-1)]
if depth == len(coeffs[0])-2:
for v in xrange(N-sum(variables[:depth])):
conditions = all(
subtotal[i]+coeffs[i][depth]*v > coeffs[i][-1]
for i in xrange(len(coeffs))
)
if conditions:
variables[depth] = v
solutions.append(tuple(variables))
else:
for v in xrange(N-sum(variables[:depth])):
variables[depth] = v
total = [subtotal[i]+coeffs[i][depth]*v for i in xrange(len(coeffs))]
find_solutions(N, coeffs, depth+1, variables, total, solutions)
if depth == 0:
return solutions
为了运行这个，生成每个方程的系数并将它们传递给函数。请记住常量的符号是倒转的！
coeffs = [
[97, 89, 42, 20, 16, 11, 2, 185],
[-98, -90, -43, -21, -17, -12, -3, -205]
]
solutions = find_solutions(20, coeffs)
print(len(solutions))
这一个在 1.6 秒内完成了 n=7 个案例，在 5.8 秒内完成了 n=8 个案例。如果您希望 n 变得非常大，我会研究任何可能的优化，但目前看起来令人满意。
剩下的问题是你的方程式之和是否总是会简化为 n1 + n2 + ... nn < N 之类的东西。如果不是这种情况，有一个简单的解决方案，但我选择不过早地过度概括代码，超出您提供给我们的示例。
最后，我想可以在 Java 或 C# 中实现相同的方法，而且速度可能会更快。如果您的一般案例预计需要更长的时间才能解决，我不介意这样做。
Reti43 的想法是正确的，但是有一个快速的递归解决方案，可以在对不等式的限制较少的假设下工作。
def solve(smin, smax, coef1, coef2):
"""
Return a list of lists of non-negative integers `n` that satisfy
the inequalities,
sum([coef1[i] * n[i] for i in range(len(coef1)]) > smin
sum([coef2[i] * n[i] for i in range(len(coef1)]) < smax
where coef1 and coef2 are equal-length lists of positive integers.
"""
if smax < 0:
return []
n_max = ((smax-1) // coef2[0])
solutions = []
if len(coef1) > 1:
for n0 in range(n_max + 1):
for solution in solve(smin - n0 * coef1[0],
smax - n0 * coef2[0],
coef1[1:], coef2[1:]):
solutions.append([n0] + solution)
else:
n_min = max(0, (smin // coef1[0]) + 1)
for n0 in range(n_min, n_max + 1):
if n0 * coef1[0] > smin and n0 * coef2[0] < smax:
solutions.append([n0])
return solutions
您可以像这样将其应用于您的原始问题，
smin, coef1 = 185, (97, 89, 42, 20, 16, 11, 2)
smax, coef2 = 205, (98, 90, 43, 21, 17, 12, 3)
solns7 = solve(smin, smax, coef1, coef2)
len(solns7)
1013
对于这样更长的问题，
smin, coef1 = 185, (97, 89, 42, 20, 16, 11, 6, 2)
smax, coef2 = 205, (98, 90, 43, 21, 17, 12, 7, 3)
solns8 = solve(smin, smax, coef1, coef2)
len(solns8)
4015
在我的 Macbook 上，这两种情况都在几毫秒内解决了。这应该可以很好地扩展到稍大的问题，但从根本上说，它在系数 N 的数量上是 O(2^N)。它实际扩展的好坏取决于附加系数的大小——更大的系数（与 smax- 相比） smin），可能的解决方案越少，它就会越快运行。
更新：从链接M.SE post上的讨论，我看出这里两个不等式之间的关系是问题结构的一部分。鉴于此，可以给出一个稍微简单的解决方案。下面的代码还包括一些额外的优化，这些优化在我的笔记本电脑上将 8 变量情况的解决方案从 88 毫秒加速到 34 毫秒。我已经在包含多达 22 个变量的示例中尝试过此操作，并在不到一分钟的时间内得到了结果，但它永远不会适用于数百个变量。
def solve(smin, smax, coef):
"""
Return a list of lists of non-negative integers `n` that satisfy
the inequalities,
sum([coef[i] * n[i] for i in range(len(coef)]) > smin
sum([(coef[i]+1) * n[i] for i in range(len(coef)]) < smax
where coef is a list of positive integer coefficients, ordered
from highest to lowest.
"""
if smax <= smin:
return []
if smin < 0 and smax <= coef[-1]+1:
return [[0] * len(coef)]
c0 = coef[0]
c1 = c0 + 1
n_max = ((smax-1) // c1)
solutions = []
if len(coef) > 1:
for n0 in range(n_max + 1):
for solution in solve(smin - n0 * c0,
smax - n0 * c1,
coef[1:]):
solutions.append([n0] + solution)
else:
n_min = max(0, (smin // c0) + 1)
for n0 in range(n_min, n_max + 1):
solutions.append([n0])
return solutions
您可以像这样将其应用于 8 变量示例，
solutions = solve(185, 205, (97, 89, 42, 20, 16, 11, 6, 2))
len(solutions)
4015
这个方案直接枚举了有界区域内的格点。由于您需要所有这些解决方案，因此获得它们所需的时间将（最多）与绑定格点的数量成正比，绑定格点的数量随维度（变量）的数量呈指数增长。
（最多 13 个整数）
这是在 1600 个内核上使用 gpgpu(opencl) 的丑陋蛮力，在 9.3 毫秒内找到 1013(7 ints) 个解决方案，包括从 gpu 到 cpu 内存的数组下载时间：
编辑： 更正了 n1、n2、n3，因为它们是 1,20,400，而不是 20,20,20。
__kernel void solver1(__global __write_only float * subCount){
int threadId=get_global_id(0);
int ctr=0;
int n1=threadId/400;
int n2=(threadId/20)%20;
int n3=threadId%20;
for(int n4=0;n4<=20;n4++)
for(int n5=0;n5<=20;n5++)
for(int n6=0;n6<=20;n6++)
for(int n7=0;n7<=20;n7++)
if (
(97*n1 + 89*n2 + 42*n3 + 20*n4 + 16*n5 + 11*n6 + 2*n7 - 185 > 0) &&
(-98*n1 - 90*n2 - 43*n3 - 21*n4 - 17*n5 - 12*n6 - 3*n7 + 205 > 0))
{ctr++;}
subCount[threadId]=ctr;
}
然后 subCount.Sum() 给出解决方案的数量。（以微秒为单位）
global work size = 8000(threads)（将 n4 加入其中将使它达到 160000 并提高性能）
本地工作大小= 160（在我的机器上效率低下，使其成为 2 的幂是最优的）
这只是一个要发送到gpu 进行编译的字符串。您只需在字符串中添加额外的循环（或只是 threadId 关系），如 n8、n9 并更改“if”正文以对它们求和。
编辑： 将 1 个额外的整数添加到方程中，即使进行更多优化（找到 4015 个解），求解时间也会增加到 101 毫秒。
__kernel void solver2(__global __write_only float * subCount){
int threadId=get_global_id(0);
int ctr=0;
int n1=threadId/160000; int c1n1=97*n1; int c2n1=-98*n1;
int n2=(threadId/8000)%20; int c1n2=89*n2; int c2n2=- 90*n2 ;
int n3=(threadId/400)%20; int c1n3=42*n3; int c2n3=- 43*n3 ;
int n4=(threadId/20)%20; int c1n4=20*n4 ;int c2n4=- 21*n4 ;
int n5=threadId%20; int c1n5=16*n5 ;int c2n5=- 17*n5 ;
int t1=c1n1+c1n2+c1n3+c1n4+c1n5;
int t2=c2n1+c2n2+c2n3+c2n4+c2n5;
for(int n6=0;n6<=20;n6++)
for(int n7=0;n7<=20;n7++)
for(int n8=0;n8<=20;n8++)
if(t1+ 11*n6 + 2*n7+6*n8 > 185 && t2 - 12*n6 - 3*n7-7*n8 > -205)
ctr++;
subCount[threadId]=ctr;
}
工作组规模 = 256
全球工作规模 = 3200000 个线程
编辑： 仅含另一个 for 循环的 9 整数版本在 1581 毫秒内找到了 c1=3 和 c2=-4 的 14k 解。将 n6 到 n9 的上限从 20 更改为 19 在 1310 毫秒内给出了相同的结果。
编辑： 添加一些来自 Reti43 和 Molecule 的答案，9整数版本使用 3.2M 线程花费 14 毫秒（但对于前 5 个整数仍然效率低下）。
__kernel void solver3(__global __write_only float * subCount){
int threadId=get_global_id(0);
int ctr=0;
int n1=threadId/160000; int c1n1=97*n1; int c2n1=-98*n1;
int n2=(threadId/8000)%20; int c1n2=89*n2; int c2n2=- 90*n2 ;
int n3=(threadId/400)%20; int c1n3=42*n3; int c2n3=- 43*n3 ;
int n4=(threadId/20)%20; int c1n4=20*n4 ;int c2n4=- 21*n4 ;
int n5=threadId%20; int c1n5=16*n5 ;int c2n5=- 17*n5 ;
int t1=c1n1+c1n2+c1n3+c1n4+c1n5;
int t2=c2n1+c2n2+c2n3+c2n4+c2n5;
int m=max( max( max( max(n1,n2),n3),n4),n5);
for(int n6=0;n6<=20-m;n6++)
for(int n7=0;n7<=20-m-n6;n7++)
for(int n8=0;n8<=20-m-n6-n7;n8++)
for(int n9=0;n9<=20-m-n6-n7-n8;n9++)
if(t1+ 11*n6 + 2*n7+6*n8 +3*n9> 185 && t2 - 12*n6 - 3*n7-7*n8-4*n9 > -205)
ctr++;
subCount[threadId]=ctr;
}
10 个整数：46k 解，35ms (c1=7,c2=-8)
11 个整数：140k 解，103ms (c1=1,c2=-2)
12 个整数：383k 解，274ms (c1=5,c2=-6)
由于硬件限制，15 个整数可能太慢了。甚至不会编译 20 个整数。
注意：
m=max( max( max( max(n1,n2),n3),n4),n5) 不是最优的，但 m=n1+n2+n3+n4+n5 是。
我现在会尝试Monte Carlo方法来克服硬件资源瓶颈。也许以某种方式收敛到解决方案编号。
伪代码：
for each inequation:
find all real roots of the equivalent equation, i.e. the zero-crossings
for each interval between two adjacent roots:
pick any number strictly inside the interval
evaluate the polynomial in that point
if the evaluated polimonial is positive:
add every integer in the interval to the list of solutions to that inequation
(treat the open-ended intervals outside the extreme roots as special cases, they may contain infinite solutions)
find the integers that are in all the lists of solutions to the individual equations
我可以使用什么算法来找到 n1, n2, ... ,n7 的所有正整数值的集合，满足以下不等式。
97n1 + 89n2 + 42n3 + 20n4 + 16n5 + 11n6 + 2n7 - 185 > 0
-98n1 - 90n2 - 43n3 - 21n4 - 17n5 - 12n6 - 3n7 + 205 > 0
n1 >= 0, n2 >= 0, n3 >=0. n4 >=0, n5 >=0, n6 >=0, n7 >= 0
例如一组 n1= 2, n2 = n3 = ... = n7 =0 使不等式成立。我如何找出所有其他值集？类似问题已发布在 M.SE.
ADDED:: 我需要概括 n 个变量（可能很大）的解决方案。我可以申请什么程序？对于 另一个 特殊情况 n=8
97n1 + 89n2 + 42n3 + 20n4 + 16n5 + 11n6 + 6n7 + 2n8 - 185 > 0
-98n1 - 90n2 - 43n3 - 21n4 - 17n5 - 12n6 - 7 - 3n8 + 205 > 0
n1 >= 0, n2 >= 0, n3 >=0. n4 >=0, n5 >=0, n6 >=0, n7 >= 0, n8 >= 0
Python 需要很长时间。 Wolfram Mathematica 显示在不到一分钟内有 4015 个解决方案。
Length[Solve[{97 n1 + 89 n2 + 42 n3 + 20 n4 + 16 n5 + 11 n6 + 6 n7 +
2 n8 - 185 > 0,
-98 n1 - 90 n2 - 43 n3 - 21 n4 - 17 n5 - 12 n6 - 7 n7 - 3 n8 +
205 > 0,
n1 >= 0, n2 >= 0, n3 >= 0, n4 >= 0, n5 >= 0, n6 >= 0, n7 >= 0,
n8 >= 0}, {n1, n2, n3, n4, n5, n6, n7, n8}, Integers]]
有 1013 个解决方案，但我不知道最有效的解决方法。
看第二个不等式，17 * n5不能大于205（否则整个左边不能为正）。这导致 n5 <= 12。通过为每个其他变量计算类似的界限，您可以将问题简化为可以使用嵌套循环快速解决的问题。
此 Java 代码打印出所有解决方案。
for (int n1 = 0; n1 <= 2; n1++) {
for (int n2 = 0; n2 <= 2; n2++) {
for (int n3 = 0; n3 <= 4; n3++) {
for (int n4 = 0; n4 <= 9; n4++) {
for (int n5 = 0; n5 <= 12; n5++) {
for (int n6 = 0; n6 <= 17; n6++) {
for (int n7 = 0; n7 <= 68; n7++) {
if (97 * n1 + 89 * n2 + 42 * n3 + 20 * n4 + 16 * n5 + 11 * n6 + 2 * n7 - 185 > 0
&& -98 * n1 - 90 * n2 - 43 * n3 - 21 * n4 - 17 * n5 - 12 * n6 - 3 * n7 + 205 > 0) {
System.out.println(Arrays.asList(n1, n2, n3, n4, n5, n6, n7));
}
}
}
}
}
}
}
}
只要
的前几项之和停止每个循环，就可以提高效率98n1 + 90n2 + 43n3 + 21n4 + 17n5 + 12n6 + 3n7
达到 205。
例如，n4循环可以替换为
for (int n4 = 0; 98 * n1 + 90 * n2 + 43 * n3 + 21 * n4 < 205; n4++)
如果你把7个循环都改成这样，那么1013个解都可以非常快的找到。
好像是整数解的线性规划。我认为已经实施了一些算法。 尝试考虑 Linear Programming Tool/Libraries for Java.
在你给出的两个例子中，我注意到了相同的模式。例如，对于第一种情况，如果将两个方程相加，您将得到
-n1 - n2 - n3 - n4 - n5 - n6 - n7 + 20 > 0
可以重新排列为
n1 + n2 + n3 + n4 + n5 + n6 + n7 < 20
这是一个很好的有界方程，您可以对其进行暴力破解。具体来说，您可以从 0 到 19 迭代 n1，从 0 到 19-n1 迭代 n2，等等。一个可能的解决方案是 (0, 0, 0, 0, 0, 0 , 0), 但我们注意到这不满足我们原来的方程。因此，只需为 (n1, n2, ..., n7) 生成所有可能的值，并只保留满足您的等式的值。将所有这些结果硬编码为
def find_solutions(N):
sols = []
for n1 in xrange(N):
for n2 in xrange(N-n1):
for n3 in xrange(N-n1-n2):
for n4 in xrange(N-n1-n2-n3):
for n5 in xrange(N-n1-n2-n3-n4):
for n6 in xrange(N-n1-n2-n3-n4-n5):
for n7 in xrange(N-n1-n2-n3-n4-n5-n6):
if (97*n1 + 89*n2 + 42*n3 + 20*n4 + 16*n5 + 11*n6 + 2*n7 - 185 > 0 and
-98*n1 - 90*n2 - 43*n3 - 21*n4 - 17*n5 - 12*n6 - 3*n7 + 205 > 0):
sols.append((n1, n2, n3, n4, n5, n6, n7))
return sols
find_solutions(20) 在 0.6 秒内找到所有 1013 个解决方案。同样，对于第二种情况，它在 2.3 秒内找到所有 4015 个解决方案。现在，这根本不容易概括，但它表明，使用聪明的方法，Python 或任何其他语言，不必很慢。
另一方面，递归允许我们将其推广到任意数量的嵌套循环，但代价是 运行ning 有点慢。
def find_solutions(N, coeffs, depth=0, variables=None, subtotal=None, solutions=None):
if variables is None:
solutions = []
subtotal = [0 for _ in xrange(len(coeffs[0]))]
variables = [0 for _ in xrange(len(coeffs[0])-1)]
if depth == len(coeffs[0])-2:
for v in xrange(N-sum(variables[:depth])):
conditions = all(
subtotal[i]+coeffs[i][depth]*v > coeffs[i][-1]
for i in xrange(len(coeffs))
)
if conditions:
variables[depth] = v
solutions.append(tuple(variables))
else:
for v in xrange(N-sum(variables[:depth])):
variables[depth] = v
total = [subtotal[i]+coeffs[i][depth]*v for i in xrange(len(coeffs))]
find_solutions(N, coeffs, depth+1, variables, total, solutions)
if depth == 0:
return solutions
为了运行这个，生成每个方程的系数并将它们传递给函数。请记住常量的符号是倒转的！
coeffs = [
[97, 89, 42, 20, 16, 11, 2, 185],
[-98, -90, -43, -21, -17, -12, -3, -205]
]
solutions = find_solutions(20, coeffs)
print(len(solutions))
这一个在 1.6 秒内完成了 n=7 个案例，在 5.8 秒内完成了 n=8 个案例。如果您希望 n 变得非常大，我会研究任何可能的优化，但目前看起来令人满意。
剩下的问题是你的方程式之和是否总是会简化为 n1 + n2 + ... nn < N 之类的东西。如果不是这种情况，有一个简单的解决方案，但我选择不过早地过度概括代码，超出您提供给我们的示例。
最后，我想可以在 Java 或 C# 中实现相同的方法，而且速度可能会更快。如果您的一般案例预计需要更长的时间才能解决，我不介意这样做。
Reti43 的想法是正确的，但是有一个快速的递归解决方案，可以在对不等式的限制较少的假设下工作。
def solve(smin, smax, coef1, coef2):
"""
Return a list of lists of non-negative integers `n` that satisfy
the inequalities,
sum([coef1[i] * n[i] for i in range(len(coef1)]) > smin
sum([coef2[i] * n[i] for i in range(len(coef1)]) < smax
where coef1 and coef2 are equal-length lists of positive integers.
"""
if smax < 0:
return []
n_max = ((smax-1) // coef2[0])
solutions = []
if len(coef1) > 1:
for n0 in range(n_max + 1):
for solution in solve(smin - n0 * coef1[0],
smax - n0 * coef2[0],
coef1[1:], coef2[1:]):
solutions.append([n0] + solution)
else:
n_min = max(0, (smin // coef1[0]) + 1)
for n0 in range(n_min, n_max + 1):
if n0 * coef1[0] > smin and n0 * coef2[0] < smax:
solutions.append([n0])
return solutions
您可以像这样将其应用于您的原始问题，
smin, coef1 = 185, (97, 89, 42, 20, 16, 11, 2)
smax, coef2 = 205, (98, 90, 43, 21, 17, 12, 3)
solns7 = solve(smin, smax, coef1, coef2)
len(solns7)
1013
对于这样更长的问题，
smin, coef1 = 185, (97, 89, 42, 20, 16, 11, 6, 2)
smax, coef2 = 205, (98, 90, 43, 21, 17, 12, 7, 3)
solns8 = solve(smin, smax, coef1, coef2)
len(solns8)
4015
在我的 Macbook 上，这两种情况都在几毫秒内解决了。这应该可以很好地扩展到稍大的问题，但从根本上说，它在系数 N 的数量上是 O(2^N)。它实际扩展的好坏取决于附加系数的大小——更大的系数（与 smax- 相比） smin），可能的解决方案越少，它就会越快运行。
更新：从链接M.SE post上的讨论，我看出这里两个不等式之间的关系是问题结构的一部分。鉴于此，可以给出一个稍微简单的解决方案。下面的代码还包括一些额外的优化，这些优化在我的笔记本电脑上将 8 变量情况的解决方案从 88 毫秒加速到 34 毫秒。我已经在包含多达 22 个变量的示例中尝试过此操作，并在不到一分钟的时间内得到了结果，但它永远不会适用于数百个变量。
def solve(smin, smax, coef):
"""
Return a list of lists of non-negative integers `n` that satisfy
the inequalities,
sum([coef[i] * n[i] for i in range(len(coef)]) > smin
sum([(coef[i]+1) * n[i] for i in range(len(coef)]) < smax
where coef is a list of positive integer coefficients, ordered
from highest to lowest.
"""
if smax <= smin:
return []
if smin < 0 and smax <= coef[-1]+1:
return [[0] * len(coef)]
c0 = coef[0]
c1 = c0 + 1
n_max = ((smax-1) // c1)
solutions = []
if len(coef) > 1:
for n0 in range(n_max + 1):
for solution in solve(smin - n0 * c0,
smax - n0 * c1,
coef[1:]):
solutions.append([n0] + solution)
else:
n_min = max(0, (smin // c0) + 1)
for n0 in range(n_min, n_max + 1):
solutions.append([n0])
return solutions
您可以像这样将其应用于 8 变量示例，
solutions = solve(185, 205, (97, 89, 42, 20, 16, 11, 6, 2))
len(solutions)
4015
这个方案直接枚举了有界区域内的格点。由于您需要所有这些解决方案，因此获得它们所需的时间将（最多）与绑定格点的数量成正比，绑定格点的数量随维度（变量）的数量呈指数增长。
（最多 13 个整数） 这是在 1600 个内核上使用 gpgpu(opencl) 的丑陋蛮力，在 9.3 毫秒内找到 1013(7 ints) 个解决方案，包括从 gpu 到 cpu 内存的数组下载时间：
编辑： 更正了 n1、n2、n3，因为它们是 1,20,400，而不是 20,20,20。
__kernel void solver1(__global __write_only float * subCount){
int threadId=get_global_id(0);
int ctr=0;
int n1=threadId/400;
int n2=(threadId/20)%20;
int n3=threadId%20;
for(int n4=0;n4<=20;n4++)
for(int n5=0;n5<=20;n5++)
for(int n6=0;n6<=20;n6++)
for(int n7=0;n7<=20;n7++)
if (
(97*n1 + 89*n2 + 42*n3 + 20*n4 + 16*n5 + 11*n6 + 2*n7 - 185 > 0) &&
(-98*n1 - 90*n2 - 43*n3 - 21*n4 - 17*n5 - 12*n6 - 3*n7 + 205 > 0))
{ctr++;}
subCount[threadId]=ctr;
}
然后 subCount.Sum() 给出解决方案的数量。（以微秒为单位）
global work size = 8000(threads)（将 n4 加入其中将使它达到 160000 并提高性能）
本地工作大小= 160（在我的机器上效率低下，使其成为 2 的幂是最优的）
这只是一个要发送到gpu 进行编译的字符串。您只需在字符串中添加额外的循环（或只是 threadId 关系），如 n8、n9 并更改“if”正文以对它们求和。
编辑： 将 1 个额外的整数添加到方程中，即使进行更多优化（找到 4015 个解），求解时间也会增加到 101 毫秒。
__kernel void solver2(__global __write_only float * subCount){
int threadId=get_global_id(0);
int ctr=0;
int n1=threadId/160000; int c1n1=97*n1; int c2n1=-98*n1;
int n2=(threadId/8000)%20; int c1n2=89*n2; int c2n2=- 90*n2 ;
int n3=(threadId/400)%20; int c1n3=42*n3; int c2n3=- 43*n3 ;
int n4=(threadId/20)%20; int c1n4=20*n4 ;int c2n4=- 21*n4 ;
int n5=threadId%20; int c1n5=16*n5 ;int c2n5=- 17*n5 ;
int t1=c1n1+c1n2+c1n3+c1n4+c1n5;
int t2=c2n1+c2n2+c2n3+c2n4+c2n5;
for(int n6=0;n6<=20;n6++)
for(int n7=0;n7<=20;n7++)
for(int n8=0;n8<=20;n8++)
if(t1+ 11*n6 + 2*n7+6*n8 > 185 && t2 - 12*n6 - 3*n7-7*n8 > -205)
ctr++;
subCount[threadId]=ctr;
}
工作组规模 = 256 全球工作规模 = 3200000 个线程
编辑： 仅含另一个 for 循环的 9 整数版本在 1581 毫秒内找到了 c1=3 和 c2=-4 的 14k 解。将 n6 到 n9 的上限从 20 更改为 19 在 1310 毫秒内给出了相同的结果。
编辑： 添加一些来自 Reti43 和 Molecule 的答案，9整数版本使用 3.2M 线程花费 14 毫秒（但对于前 5 个整数仍然效率低下）。
__kernel void solver3(__global __write_only float * subCount){
int threadId=get_global_id(0);
int ctr=0;
int n1=threadId/160000; int c1n1=97*n1; int c2n1=-98*n1;
int n2=(threadId/8000)%20; int c1n2=89*n2; int c2n2=- 90*n2 ;
int n3=(threadId/400)%20; int c1n3=42*n3; int c2n3=- 43*n3 ;
int n4=(threadId/20)%20; int c1n4=20*n4 ;int c2n4=- 21*n4 ;
int n5=threadId%20; int c1n5=16*n5 ;int c2n5=- 17*n5 ;
int t1=c1n1+c1n2+c1n3+c1n4+c1n5;
int t2=c2n1+c2n2+c2n3+c2n4+c2n5;
int m=max( max( max( max(n1,n2),n3),n4),n5);
for(int n6=0;n6<=20-m;n6++)
for(int n7=0;n7<=20-m-n6;n7++)
for(int n8=0;n8<=20-m-n6-n7;n8++)
for(int n9=0;n9<=20-m-n6-n7-n8;n9++)
if(t1+ 11*n6 + 2*n7+6*n8 +3*n9> 185 && t2 - 12*n6 - 3*n7-7*n8-4*n9 > -205)
ctr++;
subCount[threadId]=ctr;
}
10 个整数：46k 解，35ms (c1=7,c2=-8)
11 个整数：140k 解，103ms (c1=1,c2=-2)
12 个整数：383k 解，274ms (c1=5,c2=-6)
由于硬件限制，15 个整数可能太慢了。甚至不会编译 20 个整数。
注意： m=max( max( max( max(n1,n2),n3),n4),n5) 不是最优的，但 m=n1+n2+n3+n4+n5 是。
我现在会尝试Monte Carlo方法来克服硬件资源瓶颈。也许以某种方式收敛到解决方案编号。
伪代码：
for each inequation:
find all real roots of the equivalent equation, i.e. the zero-crossings
for each interval between two adjacent roots:
pick any number strictly inside the interval
evaluate the polynomial in that point
if the evaluated polimonial is positive:
add every integer in the interval to the list of solutions to that inequation
(treat the open-ended intervals outside the extreme roots as special cases, they may contain infinite solutions)
find the integers that are in all the lists of solutions to the individual equations