仕事
100 までの N について、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
次のステップとして、col 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 までの値に達する可能性があります。double はこの数値を保持できますが、そのような大きな数値の絶対精度は、予想される全体的な結果を下回るか、その程度になります。期待される結果 P は、c1!*c2!*c3!*c4! のオーダーです。(実際には、0 と 1 の間にあるはずの P/(c1!*c2!*c3!*c4!) に興味があります)。
中間結果の合計がほぼ 0 になるように、値 S の加算と減算を調整しようとしました。これは、N^N を超える中間結果を回避できるという意味で役立ちますが、これによって改善されるのは若干。また、中間結果に対数を使用して絶対数を抑えることも考えましたが、エンコードされた数値の相対的な精度は、浮動小数点数としてのエンコードによって制限され、同じ問題に遭遇すると思います。可能であれば、パフォーマンス上の理由から可変精度演算を実装しているデータ型の使用を避けたいです (現在、matlab を使用しています)。