一种计算整数网格个数的高效算法
An efficient algorithm to count the number of integer grids
考虑一个由非负整数组成的 3 x 3 正方形网格。对于每一行 i，整数之和设置为 r_i。类似地，对于每一列 j，该列中的整数总和设置为 c_j。因此，问题的实例由 6 个非负整数描述。
Is there an efficient algorithm to count how many different
assignments of integers to the grid there are given the row and column
sum constraints?
显然可以枚举所有可能的非负整数矩阵，其值最大为 sum r_i 并检查每个矩阵的约束，但这会非常慢。
例子
假设行约束为 1 2 3，列约束为 3 2 1。可能的整数网格是：
┌─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┐
│0 0 1│0 0 1│0 0 1│0 1 0│0 1 0│0 1 0│0 1 0│1 0 0│1 0 0│1 0 0│1 0 0│1 0 0│
│0 2 0│1 1 0│2 0 0│0 1 1│1 0 1│1 1 0│2 0 0│0 1 1│0 2 0│1 0 1│1 1 0│2 0 0│
│3 0 0│2 1 0│1 2 0│3 0 0│2 1 0│2 0 1│1 1 1│2 1 0│2 0 1│1 2 0│1 1 1│0 2 1│
└─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┘
在实践中，我的主要兴趣是网格的总和最多为 100，但更通用的解决方案会非常有趣。
Is there an efficient algorithm to count how many different assignments of integers to the grid there are given the row and column sum constraints?
upd 当 N 固定时（即变为常量 3），我对这个特定问题的回答是错误的。在这种情况下，它是多项式的。抱歉误导信息。
TL;DR：我认为它至少是 NP-hard。没有多项式算法，但也许有一些启发式加速。
对于 N×N 网格，您有 N 行总和方程，N 列总和方程和 N^2 非负约束：
对于N > 2，这个系统通常有不止一种可能的解决方案。因为有 N^2 个未知变量 x_ij 和 2N 个方程 => N > 2: N^2 > 2N.
您可以消除 2N - 1 个变量，只剩下一个方程式，其中 K = N^2 - (2N-1) 个变量得到总和 S。然后你将不得不处理 integer partition problem 以找出 K 项的所有可能组合以获得 S。这个问题是 NP 完全的。并且组合的数量不仅取决于项的数量K，还取决于值的顺序S.
这个问题让我想起了Simplex method. My first thought was to find just one solution using something like that method and then traverse edges of the convex to find all the possible solutions. And I was hoping that there's an optimal algorithm for that. But no, integer simplex method, which is related to integer linear programming，是NP-hard :(
我希望，您可以使用一些针对相关问题的启发式方法来加速简单的暴力解决方案。
我不知道匹配算法，但我认为找出一个不会那么困难。给定任何一个解决方案，您可以通过选择网格矩形区域的四个角，将两个对角角增加某个值并将另外两个角减少相同的值来导出另一个解决方案。该值的范围将受到每个对角线对的最低值的限制。如果您确定了所有此类范围的大小，您应该能够将它们相乘以确定所有可能的解决方案。
假设您将网格描述为熟悉的电子表格，其中列按字母顺序排列，行按数字排列，您可以在以下列表中描述所有可能的区域：
A1:B2, A1:B3, A1:C2, A1:C3, B1:C2, B1:C3, A2:B3, A2:C3, B2:C3
对于每个区域，我们根据每个对角角对的最低值列出一个范围。您可以逐步减少任一对，直到一个成员达到零，因为另一对没有上限。
选择您示例的第一个解决方案，我们可以使用此技术导出所有其他可能的解决方案。
A B C
┌─────┐
1 │0 0 1│ sum=1
2 │0 2 0│ sum=2
3 │3 0 0│ sum=3
└─────┘
3 2 1 = sums
A1:B2 - 1 solution (0,0,0,2)
A1:C2 - 1 solution (0,1,0,0)
A1:B3 1 solution (0,0,3,0)
A1:C3 2 solutions (0,1,3,0), (1,0,2,1)
B1:C2 2 solutions (0,1,2,0), (1,0,1,1)
B1:C3 1 solution (0,1,0,0)
A2:B3 3 solutions (0,2,3,0), (1,1,2,1), (2,0,1,2)
A2:C3 1 solution (0,0,3,0)
B2:C3 1 solution (2,0,0,0)
将所有解数相乘得到 2*2*3=12 个解。
如果总和很小，也许一个简单的 4 嵌套循环解决方案就足够快了？
function solve(rowsum, colsum) {
var count = 0;
for (var a = 0; a <= rowsum[0] && a <= colsum[0]; a++) {
for (var b = 0; b <= rowsum[0] - a && b <= colsum[1]; b++) {
var c = rowsum[0] - a - b;
for (var d = 0; d <= rowsum[1] && d <= colsum[0] - a; d++) {
var g = colsum[0] - a - d;
for (var e = 0; e <= rowsum[1] - d && e <= colsum[1] - b; e++) {
var f = rowsum[1] - d - e;
var h = colsum[1] - b - e;
var i = rowsum[2] - g - h;
if (i >= 0 && i == colsum[2] - c - f) ++count;
}
}
}
}
return count;
}
document.write(solve([1,2,3],[3,2,1]) + "<br>");
document.write(solve([22,33,44],[30,40,29]) + "<br>");
我已经厌倦了优化慢速选项。我得到了所有组合并更改代码只是为了获得总数。这是我能得到的最快速度：
private static int count(int[] rowSums, int[] colSums)
{
int count = 0;
int[] row0 = new int[3];
int sum = rowSums[0];
for (int r0 = 0; r0 <= sum; r0++)
for (int r1 = 0, max1 = sum - r0; r1 <= max1; r1++)
{
row0[0] = r0;
row0[1] = r1;
row0[2] = sum - r0 - r1;
count += getCombinations(rowSums[1], row0, colSums);
}
return count;
}
private static int getCombinations(int sum, int[] row0, int[] colSums)
{
int count = 0;
int max1 = Math.Min(colSums[1] - row0[1], sum);
int max2 = Math.Min(colSums[2] - row0[2], sum);
for (int r0 = 0, max0 = Math.Min(colSums[0] - row0[0], sum); r0 <= max0; r0++)
for (int r1 = 0; r1 <= max1; r1++)
{
int r01 = r0 + r1;
if (r01 <= sum)
if ((r01 + max2) >= sum)
count++;
}
return count;
}
Stopwatch w2 = Stopwatch.StartNew();
int res = count(new int[] { 1, 2, 3 }, new int[] { 3, 2, 1 });//12
int res1 = count(new int[] { 22, 33, 44 }, new int[] { 30, 40, 29 });//117276
int res2 = count(new int[] { 98, 99, 100}, new int[] { 100, 99, 98});//12743775
int res3 = count(new int[] { 198, 199, 200 }, new int[] { 200, 199, 198 });//201975050
w2.Stop();
Console.WriteLine("w2:" + w2.ElapsedMilliseconds);//322 - 370 on my computer
这对解决#P-hard 问题没有帮助（如果您允许矩阵为任意大小——请参阅下面评论中的参考资料），但有一个解决方案不等于枚举所有矩阵，而是称为 semi-standard Young tableaux 的一组较小的对象。根据您的输入，它可能会运行得更快，但仍然具有指数级的复杂性。由于它是几本代数组合学书籍或 Knuth 的 AOCP 3 中的整个章节，因此我不会在这里详细介绍，仅指向相关的维基百科页面。
想法是使用Robinson–Schensted–Knuth correspondence each of these matrix is in bijection with a pair of tableaux of the same shape, where one of the tableau is filled with integers counted by the row sum, the other by the column sum. The number of tableau of shape U filled with numbers counted by V is called the Kostka Number K(U,V)。因此，您最终会得到一个公式，例如
#Mat(RowSum, ColSum) = \sum_shape K(shape, RowSum)*K(shape, ColSum)
当然如果 RowSum == ColSum == Sum:
#Mat(Sum, Sum) = \sum_shape K(shape, Sum)^2
这是您在 SageMath 系统中的示例：
sage: sum(SemistandardTableaux(p, [3,2,1]).cardinality()^2 for p in Partitions(6))
12
这里有一些更大的例子：
sage: sums = [6,5,4,3,2,1]
sage: %time sum(SemistandardTableaux(p, sums).cardinality()^2 for p in Partitions(sum(sums)))
CPU times: user 228 ms, sys: 4.77 ms, total: 233 ms
Wall time: 224 ms
8264346
sage: sums = [7,6,5,4,3,2,1]
sage: %time sum(SemistandardTableaux(p, sums).cardinality()^2 for p in Partitions(sum(sums)))
CPU times: user 1.95 s, sys: 205 µs, total: 1.95 s
Wall time: 1.94 s
13150070522
sage: sums = [5,4,4,4,4,3,2,1]
sage: %time sum(SemistandardTableaux(p, sums).cardinality()^2 for p in Partitions(sum(sums)))
CPU times: user 1.62 s, sys: 221 µs, total: 1.62 s
Wall time: 1.61 s
1769107201498
很明显，您不会得到那么快的枚举矩阵。
根据 גלעד ברקן@ 的要求，这里有一个具有不同行和列总和的解决方案：
sage: rsums = [5,4,3,2,1]; colsums = [5,4,3,3]
sage: %time sum(SemistandardTableaux(p, rsums).cardinality() * SemistandardTableaux(p, colsums).cardinality() for p in Partitions(sum(rsums)))
CPU times: user 88.3 ms, sys: 8.04 ms, total: 96.3 ms
Wall time: 92.4 ms
10233
撇开我使用 Robinson-Schensted-Knuth 双射的其他答案不谈，这里是
另一种不需要高级组合学的解决方案，但有一些技巧
编程为任意更大的矩阵解决了这个问题。第一个想法
应该用来解决这类问题的是使用递归，通过一些记忆避免重新计算的事情
或者更好的动态规划。特别是一旦您选择了候选人
对于第一行，您将第一行减去列总和，您是
剩下同样的问题只是少了一排。避免重新计算
你存储结果的东西。你可以这样做
要么基本上在一个大table(memoization)
或者以更棘手的方式存储 n 行矩阵的所有解
并推导出具有 n+1 行的矩阵的解数（动态规划）。
这是在Python中使用记忆的递归方法：
# Generator for the rows of sum s which are smaller that maxrow
def choose_one_row(s, maxrow):
if not maxrow:
if s == 0: yield []
else: return
else:
for i in range(0, maxrow[0]+1):
for res in choose_one_row(s-i, maxrow[1:]):
yield [i]+res
memo = dict()
def nmat(rsum, colsum):
# sanity check: sum by row and column must match
if sum(rsum) != sum(colsum): return 0
# base case rsum is empty
if not rsum: return 1
# convert to immutable tuple for memoization
rsum = tuple(rsum)
colsum = tuple(colsum)
# try if allready computed
try:
return memo[rsum, colsum]
except KeyError:
pass
# apply the recursive formula
res = 0
for row in choose_one_row(rsum[0], colsum):
res += nmat(rsum[1:], tuple(a - b for a, b in zip(colsum, row)))
# memoize the result
memo[(tuple(rsum), tuple(colsum))] = res
return res
然后之后：
sage: nmat([3,2,1], [3,2,1])
12
sage: %time nmat([6,5,4,3,2,1], [6,5,4,3,2,1])
CPU times: user 1.49 s, sys: 7.16 ms, total: 1.5 s
Wall time: 1.48 s
8264346
考虑一个由非负整数组成的 3 x 3 正方形网格。对于每一行 i，整数之和设置为 r_i。类似地，对于每一列 j，该列中的整数总和设置为 c_j。因此，问题的实例由 6 个非负整数描述。
Is there an efficient algorithm to count how many different assignments of integers to the grid there are given the row and column sum constraints?
显然可以枚举所有可能的非负整数矩阵，其值最大为 sum r_i 并检查每个矩阵的约束，但这会非常慢。
例子
假设行约束为 1 2 3，列约束为 3 2 1。可能的整数网格是：
┌─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┐
│0 0 1│0 0 1│0 0 1│0 1 0│0 1 0│0 1 0│0 1 0│1 0 0│1 0 0│1 0 0│1 0 0│1 0 0│
│0 2 0│1 1 0│2 0 0│0 1 1│1 0 1│1 1 0│2 0 0│0 1 1│0 2 0│1 0 1│1 1 0│2 0 0│
│3 0 0│2 1 0│1 2 0│3 0 0│2 1 0│2 0 1│1 1 1│2 1 0│2 0 1│1 2 0│1 1 1│0 2 1│
└─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┘
在实践中，我的主要兴趣是网格的总和最多为 100，但更通用的解决方案会非常有趣。
Is there an efficient algorithm to count how many different assignments of integers to the grid there are given the row and column sum constraints?
upd 当 N 固定时（即变为常量 3），我对这个特定问题的回答是错误的。在这种情况下，它是多项式的。抱歉误导信息。
TL;DR：我认为它至少是 NP-hard。没有多项式算法，但也许有一些启发式加速。
对于 N×N 网格，您有 N 行总和方程，N 列总和方程和 N^2 非负约束：
对于N > 2，这个系统通常有不止一种可能的解决方案。因为有 N^2 个未知变量 x_ij 和 2N 个方程 => N > 2: N^2 > 2N.
您可以消除 2N - 1 个变量，只剩下一个方程式，其中 K = N^2 - (2N-1) 个变量得到总和 S。然后你将不得不处理 integer partition problem 以找出 K 项的所有可能组合以获得 S。这个问题是 NP 完全的。并且组合的数量不仅取决于项的数量K，还取决于值的顺序S.
这个问题让我想起了Simplex method. My first thought was to find just one solution using something like that method and then traverse edges of the convex to find all the possible solutions. And I was hoping that there's an optimal algorithm for that. But no, integer simplex method, which is related to integer linear programming，是NP-hard :(
我希望，您可以使用一些针对相关问题的启发式方法来加速简单的暴力解决方案。
我不知道匹配算法，但我认为找出一个不会那么困难。给定任何一个解决方案，您可以通过选择网格矩形区域的四个角，将两个对角角增加某个值并将另外两个角减少相同的值来导出另一个解决方案。该值的范围将受到每个对角线对的最低值的限制。如果您确定了所有此类范围的大小，您应该能够将它们相乘以确定所有可能的解决方案。
假设您将网格描述为熟悉的电子表格，其中列按字母顺序排列，行按数字排列，您可以在以下列表中描述所有可能的区域：
A1:B2, A1:B3, A1:C2, A1:C3, B1:C2, B1:C3, A2:B3, A2:C3, B2:C3
对于每个区域，我们根据每个对角角对的最低值列出一个范围。您可以逐步减少任一对，直到一个成员达到零，因为另一对没有上限。
选择您示例的第一个解决方案，我们可以使用此技术导出所有其他可能的解决方案。
A B C
┌─────┐
1 │0 0 1│ sum=1
2 │0 2 0│ sum=2
3 │3 0 0│ sum=3
└─────┘
3 2 1 = sums
A1:B2 - 1 solution (0,0,0,2)
A1:C2 - 1 solution (0,1,0,0)
A1:B3 1 solution (0,0,3,0)
A1:C3 2 solutions (0,1,3,0), (1,0,2,1)
B1:C2 2 solutions (0,1,2,0), (1,0,1,1)
B1:C3 1 solution (0,1,0,0)
A2:B3 3 solutions (0,2,3,0), (1,1,2,1), (2,0,1,2)
A2:C3 1 solution (0,0,3,0)
B2:C3 1 solution (2,0,0,0)
将所有解数相乘得到 2*2*3=12 个解。
如果总和很小，也许一个简单的 4 嵌套循环解决方案就足够快了？
function solve(rowsum, colsum) {
var count = 0;
for (var a = 0; a <= rowsum[0] && a <= colsum[0]; a++) {
for (var b = 0; b <= rowsum[0] - a && b <= colsum[1]; b++) {
var c = rowsum[0] - a - b;
for (var d = 0; d <= rowsum[1] && d <= colsum[0] - a; d++) {
var g = colsum[0] - a - d;
for (var e = 0; e <= rowsum[1] - d && e <= colsum[1] - b; e++) {
var f = rowsum[1] - d - e;
var h = colsum[1] - b - e;
var i = rowsum[2] - g - h;
if (i >= 0 && i == colsum[2] - c - f) ++count;
}
}
}
}
return count;
}
document.write(solve([1,2,3],[3,2,1]) + "<br>");
document.write(solve([22,33,44],[30,40,29]) + "<br>");
我已经厌倦了优化慢速选项。我得到了所有组合并更改代码只是为了获得总数。这是我能得到的最快速度：
private static int count(int[] rowSums, int[] colSums)
{
int count = 0;
int[] row0 = new int[3];
int sum = rowSums[0];
for (int r0 = 0; r0 <= sum; r0++)
for (int r1 = 0, max1 = sum - r0; r1 <= max1; r1++)
{
row0[0] = r0;
row0[1] = r1;
row0[2] = sum - r0 - r1;
count += getCombinations(rowSums[1], row0, colSums);
}
return count;
}
private static int getCombinations(int sum, int[] row0, int[] colSums)
{
int count = 0;
int max1 = Math.Min(colSums[1] - row0[1], sum);
int max2 = Math.Min(colSums[2] - row0[2], sum);
for (int r0 = 0, max0 = Math.Min(colSums[0] - row0[0], sum); r0 <= max0; r0++)
for (int r1 = 0; r1 <= max1; r1++)
{
int r01 = r0 + r1;
if (r01 <= sum)
if ((r01 + max2) >= sum)
count++;
}
return count;
}
Stopwatch w2 = Stopwatch.StartNew();
int res = count(new int[] { 1, 2, 3 }, new int[] { 3, 2, 1 });//12
int res1 = count(new int[] { 22, 33, 44 }, new int[] { 30, 40, 29 });//117276
int res2 = count(new int[] { 98, 99, 100}, new int[] { 100, 99, 98});//12743775
int res3 = count(new int[] { 198, 199, 200 }, new int[] { 200, 199, 198 });//201975050
w2.Stop();
Console.WriteLine("w2:" + w2.ElapsedMilliseconds);//322 - 370 on my computer
这对解决#P-hard 问题没有帮助（如果您允许矩阵为任意大小——请参阅下面评论中的参考资料），但有一个解决方案不等于枚举所有矩阵，而是称为 semi-standard Young tableaux 的一组较小的对象。根据您的输入，它可能会运行得更快，但仍然具有指数级的复杂性。由于它是几本代数组合学书籍或 Knuth 的 AOCP 3 中的整个章节，因此我不会在这里详细介绍，仅指向相关的维基百科页面。
想法是使用Robinson–Schensted–Knuth correspondence each of these matrix is in bijection with a pair of tableaux of the same shape, where one of the tableau is filled with integers counted by the row sum, the other by the column sum. The number of tableau of shape U filled with numbers counted by V is called the Kostka Number K(U,V)。因此，您最终会得到一个公式，例如
#Mat(RowSum, ColSum) = \sum_shape K(shape, RowSum)*K(shape, ColSum)
当然如果 RowSum == ColSum == Sum:
#Mat(Sum, Sum) = \sum_shape K(shape, Sum)^2
这是您在 SageMath 系统中的示例：
sage: sum(SemistandardTableaux(p, [3,2,1]).cardinality()^2 for p in Partitions(6))
12
这里有一些更大的例子：
sage: sums = [6,5,4,3,2,1]
sage: %time sum(SemistandardTableaux(p, sums).cardinality()^2 for p in Partitions(sum(sums)))
CPU times: user 228 ms, sys: 4.77 ms, total: 233 ms
Wall time: 224 ms
8264346
sage: sums = [7,6,5,4,3,2,1]
sage: %time sum(SemistandardTableaux(p, sums).cardinality()^2 for p in Partitions(sum(sums)))
CPU times: user 1.95 s, sys: 205 µs, total: 1.95 s
Wall time: 1.94 s
13150070522
sage: sums = [5,4,4,4,4,3,2,1]
sage: %time sum(SemistandardTableaux(p, sums).cardinality()^2 for p in Partitions(sum(sums)))
CPU times: user 1.62 s, sys: 221 µs, total: 1.62 s
Wall time: 1.61 s
1769107201498
很明显，您不会得到那么快的枚举矩阵。
根据 גלעד ברקן@ 的要求，这里有一个具有不同行和列总和的解决方案：
sage: rsums = [5,4,3,2,1]; colsums = [5,4,3,3]
sage: %time sum(SemistandardTableaux(p, rsums).cardinality() * SemistandardTableaux(p, colsums).cardinality() for p in Partitions(sum(rsums)))
CPU times: user 88.3 ms, sys: 8.04 ms, total: 96.3 ms
Wall time: 92.4 ms
10233
撇开我使用 Robinson-Schensted-Knuth 双射的其他答案不谈，这里是 另一种不需要高级组合学的解决方案，但有一些技巧 编程为任意更大的矩阵解决了这个问题。第一个想法 应该用来解决这类问题的是使用递归，通过一些记忆避免重新计算的事情 或者更好的动态规划。特别是一旦您选择了候选人 对于第一行，您将第一行减去列总和，您是 剩下同样的问题只是少了一排。避免重新计算 你存储结果的东西。你可以这样做
要么基本上在一个大table(memoization)
或者以更棘手的方式存储 n 行矩阵的所有解 并推导出具有 n+1 行的矩阵的解数（动态规划）。
这是在Python中使用记忆的递归方法：
# Generator for the rows of sum s which are smaller that maxrow
def choose_one_row(s, maxrow):
if not maxrow:
if s == 0: yield []
else: return
else:
for i in range(0, maxrow[0]+1):
for res in choose_one_row(s-i, maxrow[1:]):
yield [i]+res
memo = dict()
def nmat(rsum, colsum):
# sanity check: sum by row and column must match
if sum(rsum) != sum(colsum): return 0
# base case rsum is empty
if not rsum: return 1
# convert to immutable tuple for memoization
rsum = tuple(rsum)
colsum = tuple(colsum)
# try if allready computed
try:
return memo[rsum, colsum]
except KeyError:
pass
# apply the recursive formula
res = 0
for row in choose_one_row(rsum[0], colsum):
res += nmat(rsum[1:], tuple(a - b for a, b in zip(colsum, row)))
# memoize the result
memo[(tuple(rsum), tuple(colsum))] = res
return res
然后之后：
sage: nmat([3,2,1], [3,2,1])
12
sage: %time nmat([6,5,4,3,2,1], [6,5,4,3,2,1])
CPU times: user 1.49 s, sys: 7.16 ms, total: 1.5 s
Wall time: 1.48 s
8264346