permanent of square matrix

[Albert Chan]

(01-02-2024 10:41 PM)Albert Chan Wrote:  d := makelist(2,1,n);

I am beginning to understand why d cells = ±2

±aij - (±2) * aij = ∓aij

Next time we process the same row, we want to flip it back, so ±2 --> ∓2

∓aij - (∓2) * aij = ±aij

> permanent([[a11,a12,a13],[a21,a22,a23],[a31,a32,a33]])

Code:

(
+( a11+a21+a31)*( a12+a22+a32)*( a13+a23+a33)
−( a11-a21+a31)*( a12-a22+a32)*( a13-a23+a33)
+(-a11-a21+a31)*(-a12-a22+a32)*(-a13-a23+a33)
−(-a11+a21+a31)*(-a12+a22+a32)*(-a13+a23+a33)
) / 4

First term = product of (v = column sums)
Other terms, v -= d[j]*M[j], flipped only sign of M j-th row.
This make it very efficient, instead of summing them all from scratch.

sign flips caused *all* same column products to cancel (nice!)
We are left with permutations with different columns (all with plus sign!)
4 terms generated too many copies, thus shrink by 2^(3-1) = 4 correction. (*)

> expand(Ans)

Code:

+a11*a22*a33
+a11*a23*a32
+a12*a21*a33
+a12*a23*a31
+a13*a21*a32
+a13*a22*a31

(*) Because correction is always pow-of-2. We may eliminate variable s of ±1

t1 - t2 + t3 - t4 = -(t4 - (t3 - (t2 - t1)))

Also, M last row sign never flipped, d only need (n-1) cells

permanent(M), version 2

Code:

#cas
permanent(M) :=
BEGIN
 LOCAL n,d,j,f,v,p;
 if (n:=len(M))<2 THEN RETURN M[1][1] END
 d := makelist(2,2,n);
 f := range(1,n+1);
 v := map(sum,transpose(M));
 p := product(v);
 WHILE (j:=f[1]) < n DO
  v -= d[j]*M[j]
  p := product(v) - p;
  d[j] := -d[j]
  f[1] := 1;
  f[j] := f[j+1];
  f[j+1] := j+1;
 END;
 RETURN -p/2**(n-1)
END;
#end