permanent of square matrix

[Albert Chan]

(01-10-2024 08:32 AM)Werner Wrote:  Just do v/2, use dj=+/-1 again and double the result. If all vi are odd, this will not help ;-)

Numerical issue is not 2^(n-1) scaling factor (even for OP 9×9 matrix, factor is only 2^8=256).
It is really catastrophic cancellation of p's during its iterations.

For illustration purpose, I scaled lua code final result of -p/terms to simply return p
For base 2, this scaling served no purpose, except reduce tiny chance of overflow.
For base 10, this scaling may make numerical issue worse (some vi may be odd)

Code:

function permanent(M)
    local p,n,v = -2,#M,{}
    if n==1 then return M[1][1] end
    for i = 1,n do
        local t = 0
        for j = 1,n do t = t+M[j][i] end
        v[i] = t/2
        p = p*v[i]
    end
    for k = 1,2^(n-1)-1 do
        local j, t = bit.ffs(k), -2
        local d = bit.test(k,j) and -1 or 1
        for i = 1,n do v[i] = v[i]-d*M[j][i]; t = t*v[i] end
        p = t - p
    end
    return p
end

Debug with print((-1)^k*p) (★) inside for k loop, tested with OP 9×9 random matrix.

lua> permanent(m)
-5556546744184800
20429002103295824
55255400982040850
31542151598286920
...
314811585639751230
-119381930932063740
-81181015636135740
-81169499916475780
-73135689828461120
-241381971380131900
248282227143724100
...
528128789659814300
524276224665997300
-2888266392290130400
-2753810694416846000
...
-3101014455289180000
-3584102548532368000
-3432517605155301000
-3500631868051796000      -- final p = permanent(m)

This is probably a bad example, with final p very close to true result. (error = 3 ULP)
But if true result is tiny, we could see why this algorithm have cancellation issues.

(★) (-1)^k factor to account for p = t - p sign reversal.
In other words, if rest of t's contributed nothing, this is the final p we get.