%%HP: T(3)A(R)F(.);
DIR
@       SYMBOLIC MATRICES by Eliel Louzoun
@       translated to HP 48 by Per Stenius, Helsinki Univ. of Tech.
@       -----------------
@ This is a set of programs which handle symbolic matrices for the hp48sx.
@ The matrices are entered as a list, for example {{A B }{ C D }}. This package
@ contains programs for determinant,inverse of matrices,eigen values,
@ multiplication of matrices & multiplication by a scalar.
@ Bugs: The eigenvalue program works only if the values are all real.

det
\<< DUP SIZE \-> A B
    \<<
       IF B 1 ==                   @ det of a 1x1 matrix
       THEN A 1 GET 1 GET
       ELSE
         IF B 2 ==
         THEN A DET2               @ call basic case of 2x2 determinant
         ELSE 0
          1 B  FOR I
             A 1 GET I GET
             \-> E \<< IF E 0 SAME THEN ELSE
             A 1  I MINOR det E    @ recursive call
             * -1 1 I + ^ * +
          END \>> 
          NEXT
         END
       END
     \>>
  \>>
            @ This is an accessory program to compute A*D-B*C
DET2        @ from {{A B}{C D}}
\<<  \-> A
    \<< A 1 GET 1 GET
        A 2 GET 2 GET *
        A 1 GET 2 GET
        A 2 GET 1 GET * -
    \>>
 \>>

             @ This is also an accessory program to compute M(i,j) of
MINOR        @ an element in the matrix
\<< \-> A B C
  \<< A SIZE \-> S
     \<< 1 S FOR I
          IF B I \=/ THEN A I GET
       \-> D
     \<< 1 S FOR J
           IF J C \=/ THEN D J GET
           END NEXT
           S 1 - \->LIST
     \>>
      END NEXT S 1 - \->LIST
    \>>
  \>>
\>>

           @ Compute the inverse of the matrix in level 1.
           @ The  output is 2: B (a matrix)
           @                1: C (an algebraic expression)
           @ where inv(A) = B / C .
           @ This program use the cramer rule to compute the inverse.
inv
\<<  DUP SIZE \-> A B
  \<< 1 B FOR I
       1 B FOR J
         A J I MINOR det
         -1 I J + ^ *
       NEXT B \->LIST
     NEXT B \->LIST
     A det
  \>>
\>>

            @ matrices multiplication
            @ input 2: A
            @       1: B
            @ output 1: A*B
MMUL
\<< \-> A B
  \<< A SIZE B SIZE B 1 GET SIZE
     \-> L C R
    \<< 1 L FOR I A I GET
       \-> LI
       \<< 1 R FOR J 0
            1 C FOR K LI K GET
                 B K GET J GET
                 * +
            NEXT
          NEXT R \->LIST
        \>>
        NEXT L \->LIST
     \>>
   \>>
\>>

             @ scalar multiplication program
             @ input  2: scalar
             @        1: matrix
             @ output 1: matrix = scalar * matrix
SCMUL 
\<< DUP SIZE \-> M A S
  \<< 1 S FOR I A I GET \-> C
    \<< 1 S FOR J C J GET M *
       NEXT
    \>> S \->LIST
    NEXT S \->LIST
  \>>
\>>

            @ eliminate one real root from polynom
            @ input  3: a polynom in S [P(S)]
            @        2: the degree of this polynom
            @        1: the root we want to eliminate [a]
            @ output 1: P(S)/(S-a)
DIVP
\<< \-> P K R
  \<< P 'S' R - / 'S' K 1 - TAYLR \>>
\>>

eigf        @ This program returns the eigen function of a given matrix
\<< DUP SIZE \-> A B
  \<< B '-S' SIDN A
     add det
  \>>
\>>

           @ sum of two matrices
           @ input  2: A (matrix)
           @        1: B (matrix)
           @ output 1: A + B
add
\<<  DUP SIZE \-> B A S
  \<< 1 S FOR I A I GET B I GET
     \-> C D
      \<< 1 S FOR J C J GET
             D J GET +
         NEXT
      \>> S \->LIST
     NEXT S \->LIST
  \>>
\>>

            @ This program return an identity matrix multiplied by a scalar
            @ input  1: size (scalar)
            @        2: const (algbraic object)
            @ output 1: const * I
SIDN
\<< \-> K S
  \<< 1 K FOR I
       1 K FOR J
          I J == S *
       NEXT K \->LIST
     NEXT K \->LIST
  \>>
\>>

EIGV        @ returns the eigenvalues of a matrix
\<< \-> A
  \<< A eigf STEQ A SIZE
     1 FOR I RCEQ 'S' 0 ROOT
         RCEQ I 3 PICK DIVP
         STEQ
     -1 STEP { S EQ } PURGE
  \>>
\>>
END