permanent of square matrix

[John Keith]

(01-04-2024 08:06 PM)Albert Chan Wrote:  I am beginning to understand why d cells = ±2

...


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

That is fantastic! Much faster than any other method that I have seen. Even better, it turns out that we can eliminate the arrays d and f entirely, making the program shorter and a bit faster.

First of all, we can notice that the number of iterations performed by the program is 2^(n-1)-1, so we can change the WHILE loop to a FOR loop and keep track of which iteration we are in.

What f[1], and thus j are actually computing is which bit is flipped to generate the next Gray code. This value is the 2-adic valuation of k, where k is the index of the current iteration. See A001511 and A007814.

Determining which values of d[j] are negative is a bit harder but once again the OEIS comes to the rescue. The value of d[j] is negative iff the odd part of k is of the form 4k+3. See A091067.

Both of these values can be found by testing the least significant bit of k and if 0, shifting right until the lsb is a 1. If binary numbers are inconvenient, we can divide by 2 and test the result MOD 2. The number of shifts required is the valuation, and the remaining value after shifting is the odd part. If the odd part MOD 4 = 3, d[j] is negative, else it is positive.

The following program is a rough translation of your Prime program implementing the changes described above. The program converts the matrix to a nested list to take advantage of list processing commands.

Code:


\<< DUP 2 *                            @ Multiply matrix by 2
  AXL DUP SIZE 1.                      @ Convert matrix to nested list
  IF >                                 @ Size > 1?
  THEN SWAP TRN AXL                    @ List of columns of original matrix
  :: \GSLIST DOSUBS                    @ Column sums
  DUP \PILIST UNROT                    @ Product of column sums
  DUP SIZE 2. OVER 1. - ^ \-> n w      @ n = size, w = 2^(n-1)
    \<< 1. w 1. -                      @ 2^(n-1) - 1 iterations
      FOR k 1. k                       @ Row index = 1
        WHILE DUP 2. MOD NOT           @ While even
        REPEAT 2. / SWAP 1. + SWAP     @ Divide by 2 and increment index
        END 4. PICK ROT GET            @ Get matrix row
        SWAP 4. MOD 3. SAME            @ Odd part MOD 4 = 3?
        :: ADD :: - IFTE               @ Add or subtract sums
        DUP \PILIST 4. ROLL - UNROT    @ Update product
      NEXT DROP2                       @ Drop lists
      NEG w R\->I /                    @ Return negative product / w
    \>>
  ELSE DROP 1. GET                     @ Else drop list, get single element
  END
\>>