任务
我想计算 NxN 矩阵的永久P，N 最多为 100。我可以利用矩阵仅具有 M=4（或稍多）不同的行和列的事实。矩阵可能看起来像
A1 ... A1 B1 ... B1 C1 ... C1 D1 ... D1 |
... | r1 identical rows
A1 ... A1 B1 ... B1 C1 ... C1 D1 ... D1 |
A2 ... A2 B2 ... B2 C2 ... C2 D2 ... D2
...
A2 ... A2 B2 ... B2 C2 ... C2 D2 ... D2
A3 ... A3 B3 ... B2 C2 ... C2 D2 ... D2
...
A3 ... A3 B3 ... B3 C3 ... C3 D3 ... D3
A4 ... A4 B4 ... B4 C4 ... C4 D4 ... D4
...
A4 ... A4 B4 ... B4 C4 ... C4 D4 ... D4
---------
c1 identical cols
c 和 r 是 cols 和 rows 的多重性。矩阵中的所有值都介于 0 和 1 之间，并被编码为双精度浮点数。
算法
我尝试使用Ryser 公式来计算永久值。对于公式，首先需要计算每一行的总和，然后将所有行的总和相乘。对于上面的矩阵，这会产生
S0 = (c1 * A1 + c2 * B1 + c3 * C1 + c4 * D1)^r1 * ...
* (c1 * A4 + c2 * B4 + c3 * C4 + c4 * D4)^r4
作为下一步，删除 col 1 也是如此
S1 = ((c1-1) * A1 + c2 * B1 + c3 * C1 + c4 * D1)^r1 * ...
* ((c1-1) * A4 + c2 * B4 + c3 * C4 + c4 * D4)^r4
从 S0 中减去这个数字。
该算法继续使用所有可能的方法来删除单个和一组列，并将剩余矩阵的行和的乘积相加（删除偶数列）和减去（删除奇数列）。如果使用相同的 cols，则可以相对有效地解决该任务（例如，结果 S1 将恰好弹出 c1 次）。
问题
即使最终结果很小，中间结果 S0、S1、...的值也可以达到 N^N。double 可以保存这个数字，但这样大数字的绝对精度低于或接近预期的整体结果。预期结果 P 的顺序为 c1!*c2!*c3!*c4! （实际上我对 P/(c1!*c2!*c3!*c4!) 感兴趣，它应该位于 0 和 1 之间）。
我试图以中间结果的总和约为 0 的方式安排值 S 的加法和减法。从某种意义上说，这有助于我可以避免超过 N ^ N 的中间结果，但这只会改善情况一点。我还考虑过对中间结果使用对数来降低绝对数字 - 但编码数字的相对精度仍将受到浮点数编码的限制，我想我会遇到同样的问题。如果可能的话，出于性能原因，我想避免使用实现可变精度算术的数据类型（目前我正在使用 matlab）。