Question: How to compute the matrix function permanent with Ryser's algorithm?

The permanent of a square matrix is defined similarly to the determinant, as a sum of products of sets of matrix entries that lie in distinct rows and columns. However, where the determinant weights each of these products with a ±1 sign based on the parity of the set, the permanent weights them all with a + 1 sign.While the determinant can be computed in polynomial time by Gaussian elimination, the permanent cannot. Valiant has showed that computing permanents is #P-hard, and even #P-complete for matrices in which all entries are 0 or 1 [L.G. Valiant, The Complexity of Computing the Permanent, Theoretical Computer Science 8 (1979) 189–201]. The development of both exact and approximate algorithms for computing the permanent of a matrix is an active area of research. (see Wikipedia: computing the permanent, more details on permanent)

The best known general exact algorithm is due to H. J. Ryser (1963) [Ryser, Herbert John (1963), Combinatorial Mathematics, The Carus mathematical monographs, The Mathematical Association of America. Ryser formula ]. Ryser's formula can be evaluated using O(2^{n}*n^2) arithmetic operations, whereas O(n!*n) arithmetic operations if we compute the permanent using the definition.

There have been codes of Ryser's algorithm for computing the permanent written in C language and Matlab. I cannot find a Ryser's procedure in Maple. So I try to write a Maple procedure to compute the permanent using Ryser's algorithm. The code is as follows:

permanentRyser := proc (M::Matrix) 
      local S, T, B, m, n, s, i, j, rowsum, sum, prod, perm;
      m, n := op(1, M);
      if m <> n then
               error "expecting a square Matrix, got dimensions %1, %2", m, n
      end if;
      rowsum := 0;
      sum := 0;
      prod := 1; 
      S := combinat:-subsets([seq(1 .. m)]);
      S_nextvalue();
      while not S_finished do
             T := S_nextvalue();
             s := numelems(T);
             B := LinearAlgebra:-SubMatrix(M, [1 .. m], T);
             for i to m do
                  for j to s do
                          rowsum := rowsum+B[i][j];
                 end do;
                 prod := prod*rowsum; 
                 rowsum := 0;
             end do; 
            prod := (-1)^s*prod;
            sum := sum+prod;
            prod := 1 ;
      end do;
      perm := expand((-1)^m*sum) ;
end proc:
The last second statement "perm := expand((-1)^m*sum) ;" I mean to compute the permanent of  a matrix containing a variable, e.g. the characteristic matrix. If the matrix is numeric, then "expand" should be deleted.

Now I have two questions:

1. Suppose that A is random matrix of order 20, the time(permanentRyser(A)) is about 716 seconds and time(LinearAlgebra:-Permanent(A)) is about 66 seconds. We can see that LinearAlgebra:-Permanent(A) is much faster than permanentRyser(A).  I don't know whether the code I written is accurate and efficient. Thanks to anyone who gives a more efficient procedure of computing the permanent using Ryser's formular.

2. The source code of the function permanent is as follows.

proc (A::(`~Matrix`(square)), ` $`)
     option `Copyright (c) 1999 Waterloo Maple Inc. All rights reserved.`;
     local m, n, t1;
     m, n := op(1, A);
     if m <> n then
              error "expecting a square Matrix, got dimensions %1, %2", m, n
     end if;
     if n = 1 then t1 := A[1, 1] ;
     else
            tools/ClearRememberTable(Permanent/pminor`);
            t1 := [`$`(1 .. n)];
            t1 := Permanent/pminor(A, t1, t1, n)
      end if
end proc

I cannot understand the source code. Does the source code compute the permanent by Laplacian expansion? If the LinearAlgebra:-Permanent() compute the permanent by Laplacian expansion, then it should take more time than by using Ryser's algorithm.

Thanks all.