一种计算整数网格数量的有效算法

Question

考虑一个由非负整数组成的 3 x 3 方格。对于每一行i，整数之和设置为r_i。类似地，对于每一列j，该列中的整数总和设置为c_j。因此，问题的一个实例由6 非负整数描述。

是否有一种有效的算法来计算有多少不同 将整数分配给网格有给定行和列 总和约束？

显然，可以枚举所有可能的非负整数矩阵，其值不超过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，但更通用的解决方案会非常有趣。

Answer 1

在给定行和列总和约束的情况下，是否有一种有效的算法来计算有多少不同的整数分配给网格？

upd 当N 固定时（即变为常量3），我对这个特定问题的回答是错误的。在这种情况下，它是多项式的。对于误导性信息，我们深表歉意。

TL;DR：我认为它至少是 NP 难的。没有多项式算法，但也许有一些启发式加速。

---

对于 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的term数量，还取决于S的值的顺序。

这个问题让我想起了Simplex method。我的第一个想法是使用类似方法找到一个解决方案，然后遍历凸面的边缘以找到所有可能的解决方案。我希望有一个最佳算法。但是不，与integer linear programming 相关的整数单纯形法是NP-hard :(

我希望，对于相关问题，您可以使用某种启发式方法来加速天真的蛮力解决方案。

Answer 2

我不知道匹配算法，但我认为解决一个问题并不难。给定任何一种解决方案，您可以通过选择网格矩形区域的四个角，将两个对角角增加某个值并将其他两个角减小相同的值来得出另一个解决方案。该值的范围将受到每个对角线对的最小值的限制。如果您确定所有这些范围的大小，您应该能够将它们相乘以确定可能的总解。

假设您将网格描述为熟悉的电子表格，按字母顺序排列列，按数字排列行，您可以在以下列表中描述所有可能的区域：

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 个解决方案。

Answer 3

如果总和很小，也许一个简单的 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>");

Answer 4

#P-hard 的问题无济于事（如果您允许矩阵为任何大小 - 请参阅下面评论中的参考），但有一个解决方案并不等于枚举所有矩阵，而是一组较小的对象，称为semi-standard Young tableaux。根据您的输入，它可能会更快，但仍然具有指数级复杂性。由于它是几本代数组合学书籍或 Knuth 的 AOCP 3 中的一整章，因此我不会在这里详细介绍，仅指向相关的 wikipedia 页面。

这个想法是，使用Robinson–Schensted–Knuth correspondence，这些矩阵中的每一个都与一对相同形状的表格成双射，其中一个表格用行和计算的整数填充，另一个用列和计算.用 V 计数的数字填充的 U 形图形的数量称为Kostka NumberK(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

Answer 5

我已经厌倦了优化慢速选项。我得到所有组合并仅更改代码以获得总计数。这是我能得到的最快速度：

    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

Answer 6

除了我使用 Robinson-Schensted-Knuth 双射的其他答案，这里是 另一个不需要高级组合的解决方案，但有一些技巧 编程解决任意较大矩阵的这个问题。第一个想法 应该用来解决这类问题的方法是使用递归，通过一些记忆来避免重新计算 或者更好的动态规划。具体来说，一旦您选择了候选人 对于第一行，您将第一行减去列总和，您就是 留下同样的问题，只有少一排。避免重新计算 你存储结果的东西。你可以这样做

• 基本上都在一张大表中（记忆）

• 或者以一种更棘手的方式存储具有 n 行的矩阵的所有解 并推导出具有 n+1 行的矩阵的解数（动态规划）。

这是在 Python 中使用 memoization 的递归方法：

 # 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