问：

任务

我想计算 N 到 100 的 NxN 矩阵的永久 P。我可以利用矩阵仅具有 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 是列和行的多重性。矩阵中的所有值都介于 0 和 1 之间，并编码为双精度浮点数。

算法

我尝试使用 Ryser 公式来计算永久。对于公式，需要首先计算每行的总和，然后将所有行的总和相乘。对于上面的矩阵，这会产生

 S0 = (c1 * A1 + c2 * B1 + c3 * C1 + c4 * D1)^r1 * ... 
    * (c1 * A4 + c2 * B4 + c3 * C4 + c4 * D4)^r4

下一步，删除列 1 也会完成相同的操作

 S1 = ((c1-1) * A1 + c2 * B1 + c3 * C1 + c4 * D1)^r1 * ... 
    * ((c1-1) * A4 + c2 * B4 + c3 * C4 + c4 * D4)^r4

并从 S0 中减去这个数字。

该算法继续采用所有可能的方法来删除单个列和组列，并将剩余矩阵的行和的乘积相加（删除的列数为偶数）和减去（删除的列数为奇数）。 如果使用相同的列，则可以相对有效地解决任务（例如，结果 S1 将恰好弹出 c1 次）。

问题

即使最终结果很小，中间结果 S0、S1、...可以达到 N^N 的值。双精度可以保持这个数字，但如此大的数字的绝对精度低于或等于预期的总体结果。预期结果 P 的量级为 c1！*c2！*c3！*c4！（实际上我对 P/（c1！*c2！*c3！*c4！） 感兴趣，它应该介于 0 和 1 之间）。

我试图以中间结果的总和约为 0 的方式排列值 S 的加法和减法。从某种意义上说，这有助于避免超过 N^N 的中间结果，但这只会稍微改善一点。我还考虑过对中间结果使用对数来降低绝对数字 - 但编码数字的相对准确性仍将受到编码为浮点数的限制，我想我会遇到同样的问题。如果可能的话，出于性能原因，我想避免使用实现可变精度算术的数据类型（目前我正在使用 matlab）。