(48g/49/50) Matrix Permanent

[John Keith]

NOTE: This program is obsolete. It has been replaced by the program in post #5 which is much more efficient.

This program computes the permanent of a square matrix. The matrix may be real or complex, and on the 49 and 50, may also be exact integer or symbolic (type 29). This program is slow because computing permanents is slow. It is not practical for matrices larger than 8 x 8, and a 7 x 7 matrix may take several minutes on the physical calculator.

The program uses recursive Laplace expansion as described in this section of the Wikipedia link above. For matrix (or sub-matrix) size 2 x 2, the explicit formula is used, requiring only two multiplications and one addition. This method requires about half as many multiplications as the "brute force" method based on the definition of the permanent.

The program should be stored in a variable named 'PRMNT'. If another name is used, the name reference in the program must be changed.

Code:


\<< DUP SIZE HEAD 2 OVER
  IF SAME                                     @ Size = {2 2}?
  THEN DROP OBJ\-> DROP ROT ROT * ROT ROT * + @ 2 x 2 permanent
  ELSE \-> n
    \<< 1 COL- OBJ\-> EVAL \->LIST SWAP 1 n   @ Remove column 1 and save
      FOR k DUP k ROW- DROP PRMNT SWAP        @ Recurse over row permutations
      NEXT DROP n \->LIST * \GSLIST           @ Multiply by column 1 and sum
    \>>
  END
\>>

A shorter version for the 49 and 50:

Code:


\<< DUP SIZE HEAD 2. OVER
  IF SAME
  THEN DROP OBJ\-> DROP UNROT * UNROT * +
  ELSE \-> n
    \<< 1. COL- AXL SWAP 1. n
      FOR k DUP k ROW- DROP PRMNT SWAP
      NEXT DROP n \->LIST * \GSLIST
    \>>
  END EVAL                                    @ EVAL for symbolic matrices
\>>

Some examples:

This 7 x 7 1/0 matrix from A000794 returns 24 in about 4 minutes on my 50g.

Code:


[[1 0 0 1 1 0 0]
 [0 1 1 0 1 0 0]
 [1 0 1 0 0 1 0]
 [0 1 0 1 0 1 0]
 [0 0 1 1 0 0 1]
 [1 1 0 0 0 0 1]
 [0 0 0 0 1 1 1]]

A symbolic permanent:

Code:


[[ 'A' 'B' 'C' ]
 [ 'D' 'E' 'F' ]
 [ 'G' 'H' 'I' ]]
PRMNT

returns

'(I*E+H*F)*A+((I*D+G*F)*B+(H*D+G*E)*C)'

[Valentin Albillo]

.
Hi,

For a 5-line HP-71B recursive subprogram to compute NxN permanents and determinants using a similar procedure (expansion by minors), have a look at my PDF article which I posted about 4 years ago:

Exact Determinants and Permanents

which includes the subprogram proper and revealing examples.

Regards.
V.

[John Keith]

Very nice, thanks! You seem to be 18 years ahead of me.

I tried the idea from your article regarding the explicit computation of 3 x 3 matrices. It does result in a significant speed-up for matrices larger than 6 x 6 but at a cost of doubling the size of the program. Not to mention "stackrobatics" that would make your eyes hurt!

I did not consider giving the program the ability to compute determinants since the 49 and 50 can solve exact and symbolic determinants much faster. Might be useful for the 48g though.

[John Keith]

The new and improved(?) version, both larger and faster. Obsolete- see post #5 below!
Does direct computation of 3 x 3 matrices and also handles 2 x 2 and 1 x 1 matrices. The program uses compiled local variables to isolate the 1 x 1 and 2 x 2 cases from the main recursive program.

Real and exact matrices up to size 9 can be solved, if slowly. Complex and symbolic matrices up to size 7 or 8 can be solved, depending on the size of elements and the amount of available memory.

The program does not avoid computations with zero elements because there is no fast and simple test for zero that will work with real, complex, exact integer and symbolic matrices.

Code:


\<< DUP SIZE HEAD DUP 3
  IF \>=                                       @ Main routine for 3 * 3
  THEN DROP                                    @ and larger
    \<< DUP SIZE HEAD 3 OVER
      IF SAME
      THEN DROP OBJ\-> DROP                    @ Explode matrix onto stack
        5 PICK OVER * 3 PICK 6 PICK * +
        9 ROLLD 6 PICK * 3 PICK 5 ROLL * +
        6 ROLLD 4 ROLL * ROT ROT * + *
        ROT ROT * + ROT ROT * +                @ 3 * 3 permanent
      ELSE \-> n
        \<< 1 COL- OBJ\-> EVAL \->LIST SWAP   @ Remove column 1 and save
          1 n
          FOR k DUP k ROW- DROP \<-P EVAL SWAP @ Recurse over row permutations
          NEXT DROP n \->LIST * \GSLIST        @ Multiply by column 1 and sum
        \>>
      END
    \>> \-> \<-P
    \<< \<-P EVAL
    \>>
  ELSE 2 SAME
    { OBJ\-> DROP UNROT * UNROT * + }          @ 2 * 2 permanent
    { 1 GET } IFTE                            @ 1 * 1 permanent
  END
\>>

As before, a slightly smaller version optimized for the 49 and 50.

Code:


\<< DUP SIZE HEAD DUP 3.
  IF \>=                                       @ Main routine for 3 * 3
  THEN DROP                                    @ and larger
    \<< DUP SIZE HEAD 3. OVER
      IF SAME
      THEN DROP OBJ\-> DROP                    @ Explode matrix onto stack
        5. PICK OVER * PICK3 6. PICK * +
        9. ROLLD 6. PICK * PICK3 5. ROLL * +
        6. ROLLD 4. ROLL * UNROT * + *
        UNROT * + UNROT * +                    @ 3 * 3 permanent
      ELSE \-> n
        \<< 1. COL- AXL SWAP 1. n
          FOR k DUP k ROW- DROP \<-P EVAL SWAP @ Recurse over row permutations
          NEXT DROP n \->LIST * \GSLIST        @ Multiply by column 1 and sum
        \>>
      END EVAL                                 @ EVAL for symbolic matrices
    \>> \-> \<-P
    \<< \<-P EVAL
    \>>
  ELSE 2. SAME
    { OBJ\-> DROP UNROT * UNROT * + EVAL }     @ 2 * 2 permanent
    { 1. GET } IFTE                            @ 1 * 1 permanent
  END
\>>

[John Keith]

A new and much improved program for the HP 49 and 50 using Glynn's formula. This program is based on Albert Chan's optimization of a program by GitHub user "lesshaste". See discussion in this thread. Thanks also to Valentin Albillo, whose recent thread inspired the current flurry of activity on matrix permanents.

The program should be used in Exact mode for larger matrices to avoid overflow.

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
\>>

Updated 2024/01/08 to clean up code.
Newer version 2024/01/10, a bit smaller and simpler, now 256 bytes.
Latest version 2024/01/24, slightly larger but 33% faster.

[John Keith]

The program in post #5 above has been updated with two significant changes. Firstly, the FLASHEVAL and R~SB command have been replaced by a WHILE loop which is just as fast and only a few bytes larger. This removes the requirement for Library 256.

Second and most importantly, I realized there was a rather glaring inefficiency in the original Python program requiring a list multiplication in every iteration. That multiplication has been factored out of the loop resulting in a speedup of about 33%. The current program computes the permanent of Valentin's 12 x 12 matrix in 254 seconds (Exact mode) or 221 seconds (Approximate mode) on my 50g.