-JMB_MATRIX   Manual


 Overview
 

-The elements aij of a matrix are to be stored in column order into successive data registers
  and each matrix is identified by a control number of the form   bbb.eeerr

  where  Rbb  is the first register  ,   Ree   is the last register  and     rr     is the number of rows of the matrix.

                                 3   1   4   1                                 R11   R14   R17   R20
-For example   A =   5   9   2   6    may be stored in    R12   R15   R18   R21   respectively   and its control number =  11.02203
                                5   3   5   8                                  R13   R16   R19   R22

-"MSTO" may help you to store matrices in this way.

-Likewise, the control number of a polynomial is  bbb.eee
-If you have used "MSTO" to store its coefficients, subtract  rr E-5.
 
 
 
XROM  Function  Desciption
 08,00
 08,01
 08,02
 08,03
 08,04
 08,05
 08,06
 08,07
 08,08
 08,09
 08,10
 08,11
 08,12
 08,13
 08,14
 08,15
 08,16
 08,17
 08,18
 08,19
 08,20
 08,21
 08,22
 08,23
 08,24
 08,25
 08,26
 08,27
 08,28
 08,29
 08,30
 08,31
 08,32
 08,33
 08,34
 08,35
 08,36
 08,37
 08,38
 08,39
 08,40
 08,41
 08,42
 08,43
 08,44
 08,45
 08,46
 08,47
 08;48
 08,49
 08,50
 08,51
 08,52
 08,53
 08,54
 08,55
 08,56
 -JMB MATRX
 "CRMAT"
 "MCO"
 "MSTO"
 "RANM"
 "RANSYM"
 "TRN"
 "TRN2"
 -MTR FNS
 "EXPM"
 "GRVL"
 "LNM"
 "M-"
 "M+"
 "M*"
 "MINV"
 "MNORM"
 "MPOW"
 "MROOT"
 "PMAT"
 "SQRTM"
 "TM*"
 "TRACE"
 -EIGENVAL
 "DFL"
 "JCB"
 "JCBZ"
 "PCAR"
 "PMIN"
 -LNR SYS
 "D0"
 "D3"
 "D4"
 "D5"
 "DET"
 "LS"
 "LS3"
 "LS3-"
 "LS3+"
 -MAGIC MTR
 "MG4DHC"
 "MGCB"
 "MGSQ"
 "PMGCB"
 "PMGSQ"
-POLYNOM
 "BRST"
 "COHN"
 "DIV"
 "FCTR"
 "IRR?"
 "P2"
 "P3"
 "P4"
 "PL"
 "PPRD"
 "PR?"
 Section Header
 Creating a Matrix
 Copying a Matrix
 Storing the Elements of a Matrix & Control Number
 Random Matrices
 Random Symmetric Matrices
 Transpose of a Matrix
 Transpose of a Symmetric Matrix
 Section Header
 Exponential of a Matrix
 Pseudo-Inverse of a Matrix ( Gréville's Method )
 Logarithm of a Matrix
 Subtraction of 2 Matrices
 Addition of 2 Matrices
 Multiplication of 2 Matrices
 Inverse of a Matrix
 Euclidean Norm of a Matrix
 Power of a Matrix
 P-th Root of a Matrix
 Matrix Polynomial
 Square Root of a Matrix
 Product of a Matrix Transpose by another Matrix
 Trace of a Square Matrix
 Section Header
 Deflation Method
 Jacobi's Method
 Jacobi's Method for Complex Matrices
 Characteristic Polynomial
 Minimal Polynomial
 Section Header
 A Subroutine used by "D5"
 Determinant of order 3
 Determinant of order 4
 Determinant of order 5
 Determinant of order n
 Linear Systems
 Linear Systems ( variant )
 Underdetermined Linear Systems
 Overdetermined Linear Systems
 Section Header
 4-Dimensional Magic Hypercubes of odd orders
 Magic Cubes
 Magic Squares
 Pandiagonal & Panmagic Cubes
 Panmagic Squares
 Section Header
 Bairstow's Method
 Cohn's Irreducibility Criterion
 Polynomial Euclidean Division
 Factorization over Z[X]
 A Better Irreducibility Criterion
 Quadratic Equations
 Cubic Equations
 Quartic Equations
 Evaluation of a Polynomial P(x)
 Product of 2 Polynomials
 Primality Test ( Non-negative integers )

 

Creating a Matrix
 

Data Registers:        R00 = function name                   ( Register R00 is to be initialized before executing "CRMAT" )

                                  and the registers containing the matrix elements ( no other register is used )
Flags: /
Subroutine:      The function f  ( assuming i is in X-register and j is in Y-register upon entry )

        If needed, you have also 1 in L-register,
        n = the number of rows in N-register and
        k = the pointer value of the current element in M-register ( from 1 to n.m , counted in column order )
 
 
      STACK        INPUTS      OUTPUTS
           X      bbb.eeerr             /

   where bbb.eeerr is the control number of the array.  rr = number of rows.

Example:   Store the elements of the 3x4 matrix  [ai,j]  defined by  ai,j = i2 + j3  into registers R01 thru R12.

  01  LBL "T"      defines the function  f(i,j) = i2 + j3
  02  X^2
  03  X<>Y
  04  3
  05  Y^X
  06  +
  07  RTN

      "T"    ASTO 00
1.01203  XEQ "CRMAT"   and it yields
 

      2   9   28  65                  R01  R04  R07  R10
      5  12  31  68        =        R02  R05  R08  R11       respectively
     10 17  36  73                  R03  R06  R09  R12
 

Note:

-If your matrix is actually a vector ( only one row )  you can key in  bbb.eee  instead of  bbb.eee01
 

Copying a Matrix
 

-You want to copy the elements of a matrix ( control number = bbb.eeerr1 )  into a block of registers that starts with Rbb2
 
 
      STACK        INPUTS      OUTPUTS
           Y     bbb.eeerr1      bbb.eeerr1
           X         bbb2      bbb.eeerr2

 
Example:

           2  7  1  3          R01  R04  R07  R10
   A =  1  9  4  2    =    R02  R05  R08  R11    and you want to copy this matrix in registers  R21 ...etc...
           4  6  2  1          R03  R06  R09  R12
 

   1.01203  ENTER^                                                                  R21  R24  R27  R30
        21      XEQ "MCO"  >>>>  21.03203  and  A is also in     R22  R25  R28  R31
                                                                                                  R23  R26  R29  R32

Warning:

-The 2 blocks of registers cannot overlap unless  bbb2 <  bbb1
 

Storing the Elements of a Matrix & Control Number
 
 
 
      STACK        INPUTS      OUTPUTS
           X           bbb      bbb.eeerr

 
Example:    You want to store

              3   1   4   1
     A =   5   9   2   6    starting with register  R11
              5   3   5   8

-You enter the different elements in column order.
-When the first column is stored, simply press R/S without any digit entry,
  the next column indexes will be automatically incremented.
-When all the elements are stored, press R/S and you'll get the control number of the matrix!

-More explicitly:

  11  XEQ "MSTO"  the HP-41 displays   A1,1=?

   3   R/S                   -------------------   A2,1=?
   5   R/S                   -------------------   A3,1=?
   5   R/S                   -------------------   A4,1=?        the first column is stored so:
Simply  R/S               -------------------   A1,2=?        now, the HP-41 knows the matrix has 3 rows

   1   R/S                   -------------------   A2,2=?
   9   R/S                   -------------------   A3,2=?
   3   R/S                   -------------------   A1,3=?

   4   R/S                   -------------------   A2,3=?
   2   R/S                   -------------------   A3,3=?
   5   R/S                   -------------------   A1,4=?

   1   R/S                   -------------------   A2,4=?
   6   R/S                   -------------------   A3,4=?
   8   R/S                   -------------------   A1,5=?       all the elements are stored so:

Simply R/S   >>>>    control number  =  11.02203      the first register is R11 , the last one is R22 and there are 3 rows.

-The number of rows must be < 100
-If you put a number like PI in X-register , set flag F22 before pressing R/S   ( otherwise, the HP-41 would think there is no digit entry... )

-You can use registers X and Y to key in an element like  2+LN(3):      3  LN  2  +  R/S
 

Random Matrices
 

-"RANM" stores pseudo-random integers ( between 1 and N ) into registers Rbb , ........ , Ree
-No other register is used.
 
 
      STACK        INPUTS      OUTPUTS
           Y       bbb.eeerr             /
           X            N             /

   where bbb.eeerr is the control number of the array

Example:    Store random integers between 1 and 12 into registers   R01 , R02 , ........ , R07

  1.007  ENTER^
    12     XEQ "RANM"  gave   R01 = 10   R02 = 5   R03 = 8   R04 = 6   R05 = 6   R06 = 8   R07 = 10  ...  when I did it !

-Since the date & time are used to initialize the random number generator, you will probably get different values...
 

Random Symmetric Matrices
 

-"RANSYM" directly stores pseudo-random integers ( between 1 and N ) in registers  R01 thru  Rn2
-No other register is used.
 
 
      STACK        INPUTS      OUTPUTS
           Y             n             /
           X            N             /

Example:   Store a random symmetric matrix of order 5 - elements between 1 and 12 - in registers R01 thru R25

   5    ENTER^
  12   XEQ "RANSYM"  gave ( when I executed "RANSYM" )

   R01   R06   R11   R16   R21                8     2    12    9    11
   R02   R07   R12   R17   R22                2    12   11    8    11
   R03   R08   R13   R18   R23      =       12   11   12   10    6
   R04   R09   R14   R19   R24                9     8    10   12    6
   R05   R10   R15   R20   R25               11   11    6     6     4

-Like "RANM", "RANSYM" uses the date & time to initialize the random number generator, so you will probably get different values...
 

Transpose of a Matrix
 

-The transpose of a nxm matrix A = [ ai j ] is the mxn matrix tA = [ bi j ] defined by  bi j = aj i
-"TRN" stores the transpose of a matrix  in a different block of registers.
-The 2 blocks cannot overlap.
 
 
      STACK        INPUTS      OUTPUTS
           Y     bbb.eeerr1   rr1+bbb.eeerr1
           X         bbb2      bbb.eeerr2

 
Example:

            2  7  1  3         R01  R04  R07  R10
    A =  1  9  4  2    =   R02  R05  R08  R11    and you want to store tA in registers  R21 ...etc...
            4  6  2  1         R03  R06  R09  R12

   1.01203  ENTER^
        21      XEQ "TRN"  >>>>   21.03204  and

             2  1  4          R21  R25  R29
    tA =  7  9  6    =    R22  R26  R30
             1  4  2          R23  R27  R31
             3  2  1          R24  R28  R32
 

Transpose of a Symmetric Matrix
 

-Unlike "TRN" , "TRN2"  replaces a square matrix by its transpose.
 
 
      STACK         INPUT      OUTPUT
           X       bbb.eeerr      bbb.eeerr

 
Example:

              3  1  4          R01  R04  R07
     A =   1  5  9    =    R02  R05  R08
              2  6  5          R03  R06  R09

  1.00903  XEQ "TRN2"  >>>>   1.00903   and

    R01  R04  R07             3  1  2
    R02  R05  R08     =      1  5  6    =  tA
    R03  R06  R09              4  9  5
 

Exponential of a Matrix
 

-The exponential of a nxn matrix A is defined as

   exp(A) = I + A + A2/2! + A3/3! + ..... + Ak/k! + ....                ( I = the Unit matrix = the Identity matrix:  1 in the diagonal, 0 elsewhere )
 
 
      STACK        INPUTS      OUTPUTS
           X     bbb.eeerr1      bbb.eeerr2

      with  bbb > 005

Example:

             1   2   3           R11   R14   R17
    A =   0   1   2     =    R12   R15   R18   respectively   ( control number = 11.01903 )
             1   3   2           R13   R16   R19

   11.01903  XEQ "EXPM"  >>>>   20.02803  = control number of  exp(A)       ( in 9mn17s )

                    19.45828375     63.15030507   66.98787675             R20   R23   R26
    Exp(A) =  8.534640269    32.26024414   33.27906416      =     R21   R24   R27     respectively
                    16.63953207     58.45323648   61.70173665             R22   R25   R28

Notes:

-The elements of exp(A) will be accurately computed if all the elements of A are non-negative.
-Otherwise - especially if some elements are large negative numbers - the results may even be meaningless because of cancellation of leading digits !

-If it's possible, try to diagonalize A ( i-e find a regular matrix B so that B-1A.B is diagonal ) and apply the formula  exp ( B-1AB ) = B-1exp(A).B
-The exponential of a diagonal matrix A = [ aij ] is a diagonal matrix too where the diagonal elements are exp(aii)
 

Pseudo-Inverse of a Matrix ( Gréville's Method )
 

-If A is a nxm matrix ( n rows, m columns ), "GRVL" computes the pseudoinverse A+ , also called Moore-Penrose inverse, of the matrix A
  A+ is the unique mxn matrix ( m rows, n columns ) such that:  A A+ A = A ,  A+ A A+  = A+ ,  A A+ and  A+ A  symmetric

-The method of Greville gradually builds the pseudo-inverse row by row:

-Let

        ak = the k-th column of A
        Ak = the matrix built with the k first columns of A
       A+k  = the pseudoinverse at step k

-We start with  A+1 = ta1 / ( ta1 a1 )   if   a1 # 0
                  or  A+1 = ta1                   if   a1 = 0     and then,  for k = 2 , 3 , ..... , m

     Let    dk = A+k-1  ak
              ck = ak - Ak-1 dk
              bk = tck / ( tck ck )   if   ck # 0      or      bk = tdk A+k-1 / ( 1 +  tdk dk )   if  ck = 0

   >>>   A+k  = [ A+k-1 - dk bk ]           In other words,  A+k  is obtained by placing  the row   bk  under the matrix  A+k-1 - dk bk
                          [      bk      ]

-Finally,   A+ = A+m
 

Data Registers:           R00 thru R09: temp                        ( The " Registers" are to be initialized before executing "GRVL" )

                                    Rbb thru    Ree  the coefficients of the matrix A, in column order   bb > 09

                                      Registers RBB to RBB+N-1  are also used with N = m+n-1+2(m-1)n
Flags: /
Subroutines:  "TRN"  "M*"  "TM*"
 
 
      STACK        INPUTS      OUTPUTS
           Y       bbb.eeerr            /
           X           BBB      bbb'.eee'rr'

  where  bbb.eeerr = control number of the matrix A        bbb > 009
             BBB is chosen so that the blocks of registers do not overlap, for instance:  BBB = eee + 1
    and   bbb'.eee'rr' = control number of the matrix A+ = R06

Example:

                 1  1  4  2                    R10  R13  R16  R19
      A =     0  1  2  3         =         R11  R14  R17  R20      ( respectively )      control number =  10.02103
                 3  2  6  7                    R12  R15  R18  R21

  and if you choose  BBB = 22

   10.02103   ENTER^
        22         XEQ "GRVL"  >>>>   28.03904

                   R28   R32   R36             -21/112   -85/112   43/112
-We get      R29   R33   R37     =         7/112      23/112   -9/112         =   A+   ( with some roundoff errors of course ).
                  R30   R34   R38                49/112      1/112   -15/112
                  R31   R35   R39               -35/112    29/112    13/112

-In this example,  A A+ = Id3
 

Logarithm of a Matrix
 

-If  || A - I || < 1 , the logarithm of a nxn matrix A is defined by the series expansion:

       Ln(A) = ( A - I ) - ( A - I )2/2 + ( A - I )3/3 -  ( A - I )4/4 + ......               ( I = the Unit matrix = the Identity matrix:  1 in the diagonal, 0 elsewhere )
 
 
      STACK        INPUTS      OUTPUTS
           X       bbb.eeerr1      bbb.eeerr2

      with  bbb > 005

Example:

               1.2   0.1   0.3          R11   R14   R17
      A =   0.1   0.8   0.1    =    R12   R15   R18    respectively   ( control number = 11.01903 )
               0.1   0.2   0.9          R13   R16   R19

   11.01903  XEQ "LNM"  >>>>   20.02803  = control number of  Ln(A)

                   0.167083396    0.069577923   0.287707999          R20   R23   R26
    Ln(A) =  0.097783005   -0.240971674   0.103424021    =    R21   R24   R27     respectively
                   0.086500972    0.235053124  -0.131906636          R22   R25   R28

-In this example,   || A - I || = 0.5099...  < 1

Notes:

-If  || A - I || < 1   Exp(Ln(A)) = A
-If  || A || < Ln 2   Ln(Exp(A)) = A
 

Subtraction of 2 Matrices
 

 "M-" calculates the difference between 2 matrices A & B and returns the control number of A-B
 
 
      STACK        INPUTS      OUTPUTS
           Z     bbb.eeerr1             /
           Y     bbb.eeerr2             /
           X         bbb3     bbb.eeerr3

                            3  1  4       R01  R03  R05                   2  7  1        R07  R09  R11
Example:    A =  1  5  9   =  R02  R04  R06   and  B =   8  2  8   =  R08  R10  R12  respectively

-The control number of A is  1.00602
-The control number of B is   7.01202     if, for instance, you choose to store A+B in registers R15 .....

      1.00602  ENTER^
      7.01202  ENTER^
          15       XEQ "M-"  >>>>  15.02002  and

                 1  -6   3       R15  R17  R19
    A-B =  -7   3   1  =   R16  R18  R20   respectively

-The destination block of registers may be the same as the one of matrix A or B.
 

Addition of 2 Matrices
 

 "M+" calculates the sum of 2 matrices A & B and returns the control number of A+B
 
 
      STACK        INPUTS      OUTPUTS
           Z     bbb.eeerr1             /
           Y     bbb.eeerr2             /
           X         bbb3     bbb.eeerr3

                            3  1  4       R01  R03  R05                   2  7  1        R07  R09  R11
Example:    A =  1  5  9   =  R02  R04  R06   and  B =   8  2  8   =  R08  R10  R12  respectively

-The control number of A is  1.00602
-The control number of B is   7.01202     if, for instance, you choose to store A+B in registers R15 .....

      1.00602  ENTER^
      7.01202  ENTER^
          15       XEQ "M+"  >>>>  15.02002  and

                  5  8  5        R15  R17  R19
    A+B =   9  7 17  =   R16  R18  R20   respectively

-The destination block of registers may be the same as the one of matrix A or B.
 

Multiplication of 2 Matrices
 

 "M*" calculates the product of 2 matrices A & B and returns the control number of A*B
 
 
 
      STACK        INPUTS      OUTPUTS
           Z     bbb.eeerr1        rr1 = rr3
           Y     bbb.eeerr2        rr1 = rr3
           X         bbb3      bbb.eeerr3

 
Example:  Calculate  C = A.B where

                                                        3  1
                 2  7  1  3                         4  2
         A =  1  9  4  2                B =   7  5        assuming   A is stored in registers R01 thru R12         and choosing R26
                 4  6  2  1                         2  6                         B is stored in registers R15 thru R22         as the first register of C

-In other words,

     R01  R04  R07  R10                2  7  1  3                                   R15  R19         3  1
     R02  R05  R08  R11       =       1  9  4  2               and              R16  R20   =    4  2      respectively.
     R03  R06  R09  R12                4  6  2  1                                   R17  R21         7  5
                                                                                                      R18  R22         2  6

-Key in:           1.01203  ENTER^
                     15.02204  ENTER^
                          26        XEQ "M*"   >>>>   26.03103   =   the control number of the matrix C and the result is:

                47  39       R26  R29
       C =   71  51  =   R27  R30
                52  32       R28  R31

-The number of columns of the first matrix must equal the number of rows of the second matrix.
 

Inverse of a Matrix
 

 "MINV" can invert up to a 17x17 matrix.
  Gauss-Jordan elimination - also called the "exchange method" - is used.

 Here, the first element of the matrix must be stored into R06. ( R01 thru R05 are used for control numbers of different rows and columns.)
 You put the order of the matrix into X-register and XEQ "MINV"
 the determinant is in X-register and in R00 and the inverse matrix has replaced the original one ( in registers R06 ... etc ... )

 If flag F01 is clear, Gaussian elimination with partial pivoting is used.
 If flag F01 is set, the pivots are the successive elements of the diagonal.
 
 
      STACK        INPUTS      OUTPUTS
           X             n          det A

  where  n is the order of the matrix A

Example:   Let's take the 5x5 Pascal's matrix

          1  1  1   1  1                 R06  R11  R16  R21  R26
          1  2  3   4   5                R07  R12  R17  R22  R27
          1  3  6  10 15      =       R08  R13  R18  R23  R28        respectively
          1  4 10 20 35               R09  R14  R19  R24  R29
          1  5 15 35 70               R10  R15  R20  R25  R30

    5  XEQ "MINV"     the determinant ( 1 in this example ) will be in X-register and in R00 and the inverse matrix:

          5  -10  10  -5   1                    R06  R11  R16  R21  R26
       -10  30 -35  19 -4                     R07  R12  R17  R22  R27
        10 -35  46 -27  6         =          R08  R13  R18  R23  R28         respectively
        -5   19 -27  17 -4                     R09  R14  R19  R24  R29
         1   -4    6   -4   1                     R10  R15  R20  R25  R30

-If you try to invert a Pascal's matrix of order 16 for instance with flag F01 cleared, you'll obtain great roundoff errors ( even with an HP48 )
 because these matrices are very ill-conditioned. But if you set F01, all the coefficients of the inverse will be computed exactly !
 

Euclidean Norm of a Matrix
 

- "MNORM" computes  || A || = ( SUMi,j a2i,j )1/2
 
 
      STACK        INPUT      OUTPUT
           X       bbb.eeerr         || A ||

 
Example:

            2  7  1  3         R01  R04  R07  R10
    A =  1  9  4  2   =    R02  R05  R08  R11
            4  6  2  1         R03  R06  R09  R12
 

   1.01203  XEQ "MORM"  >>>>  || A || = 14.8996643
 

Power of a Matrix
 

-"MPOW" calculates Ap  where A is a square matrix and p is a positive integer   ( A0 = I = Identity Matrix )
-If p is a negative integer, compute the inverse A-1 first, and then (A-1) -p
-It uses "MCO" & "M*" as subroutines
 
 
      STACK        INPUTS      OUTPUTS
           Y      bbb.eeerr1            /
           X          p > 0      bbb.eeerr2

    ( bbb > 003 )

Example:

          1  4  9          R11  R14  R17
  A =  3  5  7    =    R12  R15  R18    and you want to compute  A^7
          2  1  8          R13  R16  R19

  11.01903  ENTER^
       7          XEQ "MPOW"  >>>>  20.02803 = the control number of A7

             7851276   8652584   31076204            R20  R23  R26
   A7 =   8911228  9823060   35267932     =     R21  R24  R27
             5829472  6422156   23076808             R22  R25  R28
 

P-th Root of a Matrix
 

-The principal p-th root of a non-singular matrix A ( det A # 0 )  may be computed by the algorithm:

    M0 = A      Mk+1 = Mk [ ( 2.I + ( p - 2 ) Mk ) ( I + ( p - 1) Mk ) -1 ] p      where  I  is the Identity matrix
    X0 =  I       Xk+1  = Xk  ( 2.I + ( p - 2 ) Mk ) -1 ( I + ( p - 1) Mk )

    Mk  tends to  I        as k tends to infinity
    Xk   tends to A 1/p  as k tends to infinity

-The convergence is quadratic if A has no negative real eigenvalue.
 

Data Registers:         R00 = n > 1                                   ( Registers R00 & R01 thru Rn2 are to be initialized before executing "MROOT" )

                       R01   thru     Rn2 =  the coefficients of the matrix A, in column order.                Rn2+1 to R5n2+1 temp  ( R5n2+1 = p )

          >>>>   When the program stops,  A1/p  is in registers  R01 thru Rn2

Flag:  F10
Subroutines:    "M*"   &   "LS3"
 
 
      STACK        INPUTS      OUTPUTS
           X             p             /

 
Example:      Calculate  A1/3  for the following matrix:

              4  2  3          R01  R04  R07
    A  =   3  2  5    =    R02  R05  R08    respectively.
              2  1  4          R03  R06  R09

  n = 3   STO 00

  p = 3   XEQ "MROOT"   the HP-41 displays the successive || Mk - I ||  and eventually,
 

                           R01  R04  R07              1.437771414   0.396760708   0.231139046
    A^(1/3)   =     R02  R05  R08     =       0.421053139   0.977706016   0.944423239
                           R03  R06  R09              0.270151310   0.119249482   1.469423763

Notes:

-If A has a negative eigenvalue and if p is odd,  the program works in some cases, but the convergence is slower. It may even seem to diverge in the first iterations.
-Newton's method  Xk+1 = (1/p) [ ( p - 1 ) Xk + A/Xkp-1 ]  works sometimes if A is very well conditioned. Otherwise, there is a numerical instability.
 

Another Checking:
 

                    1  2  4  7                             0.624151409   -0.207980111   0.372153250    0.846900264
                    2  4  1  9                             0.678031718    1.367346542  -0.226655590   -0.030942559
       A   =     4  1  6  3         A^(1/5)  =   0.533389477     0.168745809   1.273404725   -0.262412182
                    1  4  2  9                           -0.224295471     0.139205813   0.171840655    1.642919190
 

-The turbo mode is recommended...
 

Matrix Polynomial
 

-Given a polynomial  p(x) = am xm + am-1 xm-1 + ............... + a0     and a  nxn matrix A,
 "PMAT" calculates  p(A) = am Am + am-1 Am-1 + ............... + a0 I   ( I = Identity matrix )
 

Data Registers:       R00 to R03: temp               ( Registers Rbb thru Ree & Rbb' thru Ree' are to be initialized before executing "PMAT" )

                                Rbb = am  .......     Ree = a0
                                Rbb'   thru     Ree' =  the coefficients of the matrix A, in column order.

                                  the 2n2  registers from  RBB: temp

Flag:  F07
Subroutines:  "MCO"  "M*"
 
 
      STACK        INPUTS      OUTPUTS
           Z      bbb.eee             /
           Y     bbb'.eee'nn             /
           X         BBB     BBB.EEEnn

     where     bbb.eee    = control number of the polynomial p
       and     bbb'.eee'nn = control number of the  matrix A

Example:     p(x) = 2 x4 - x3 + 3 x2 - 4 x + 5

                 4  2  3
       A  =   3  2  5
                 2  1  4

-If you choose  bbb = 010  &   bbb' = 015    store  2  -1  3  -4   5   into    R10  R11  R12  R13 R14    ( control number = 10.014 )

              4  2  3            R15  R18  R21
    and    3  2  5    into   R16  R19  R22    respectively  ( control number = 15.02303 )
              2  1  4            R17  R20  R23

-With, say  BBB = 024

     10.014    ENTER^
   15.02303  ENTER^
       24         XEQ "PMAT"  >>>>   24.03203

                     3548  1887  4705            R24  R27  R30
     p(A)  =    3727  1987  4962     =     R25  R28  R31
                     2539  1351  3385            R26  R29  R32
 

Square Root of a Matrix
 

-Given a non-singular nxn matrix A  ( det A # 0 ), "SQRTM" uses the iterations:

    M0 = A      Mk+1 = (1/2) [ I + ( Mk + Mk-1 ) / 2 ]      where  I  is the Identity matrix
    X0 = A      Xk+1 = Xk ( I + Mk-1 ) / 2

-If A has no negative real eigenvalue,  Xk quadratically tends to  A1/2  and  Mk tends to I  as k tends to infinity.
-One could think to use Newton's method, but there is usually a numerical instability.

-The coefficients of A are to be stored column/column from register R01.
-When the program stops, A1/2  has replaced A in registers R01 to Rn2
 

Data Registers:         R00 = n                                    ( Registers R00 & R01 thru Rn2 are to be initialized before executing "SQRTM" )

                                    R01   thru     Rn2 =  the coefficients of the matrix A, in column order.                Rn2+1 to R4n2  temp

Flag:  F07
Subroutines:  "MCO"  "M*"  "LS3"
 
 
      STACK        INPUTS      OUTPUTS
           X             /         1.eee,nn

   with  eee = n2

Example:

              4  2  3            R01  R04  R07
    A  =   3  2  5     =     R02  R05  R08    respectively.
              2  1  4            R03  R06  R09

  n = 3   STO 00

  XEQ "SQRTM"   >>>>   X =  1.00903  = control number and we get:

                      1.794981016   0.656367529   0.540425261
    sqrt(A)  =   0.772143191   1.061374174   1.582989342
                      0.501888909   0.231634627   1.833600670
 

Product of a Matrix Transpose by another Matrix
 

-It is sometimes useful to multiply the transpose of a matrix A by another matrix B directly. ( product tA.B )
 
 
      STACK       INPUTS      OUTPUTS
           Z      bbb.eeerr1            rr3
           Y      bbb.eeerr2            rr3
           X         bbb3      bbb.eeerr3

-We must have  rr1  =  rr2  , the 2 matrices must have the same number of rows.

Example:  Calculate  C = tA.B where

                2  1  4                           3  1
                7  9  6                           4  2
        A =  1  4  2                  B =   7  5        assuming   A is stored in registers R01 thru R12         and choosing R26
                3  2  1                           2  6                         B is stored in registers R15 thru R22         as the first register of C

        R01  R05  R09                2  1  4                                   R15  R19         3  1
        R02  R06  R10       =       7  9  6               and              R16  R20   =    4  2      respectively.
        R03  R07  R11                1  4  2                                   R17  R21         7  5
        R04  R08  R12                3  2  1                                   R18  R22         2  6

-Key in:           1.01204  ENTER^
                     15.02204  ENTER^
                          26        XEQ "TM*"   >>>>   26.03103   =   the control number of the matrix C and the result is:

                47  39       R26  R29
       C =   71  51  =   R27  R30
                52  32       R28  R31
 

Trace of a Square Matrix
 

-The trace of a square matrix equals the sum of its diagonal elements.
 
 
      STACK        INPUT      OUTPUT
           X       bbb.eeerr         Tr(A)

 
Example:

             1  2  4        R03  R06  R09
    A =   3  5  7   =   R04  R07  R10
             7  9  8        R05  R08  R11

   3.01103  XEQ "TRACE"  >>>>  Tr(A) = 14
 

Deflation Method
 

-If k1 is an eigenvalue , U1 the corresponding unit eigenvector of A and U1' the corresponding unit eigenvector of  tA
  then A and the new matrix  A' =  A - k1 ( U1 tU1' ) / ( U1.U1' )   have the same eigenvalues and eigenvectors , except that k1 has been replaced by 0.
-Thus, the 2nd dominant eigenvalue of A is now the 1st dominant eigenvalue of A'

-"DFL" uses this method to calculate all the eigenvalues and eigenvectors , provided they are real and verify   | k1 | >  | k2 | > ....... >  | kn | > 0

  1°) "DFL" calculates k1 and U1  ( with flag F02 clear )  and the program stops ( line 87 )

     >>>>     k1 is in X-register
     >>>>     U1  is in registers R01 thru Rnn

  2°)  If you press R/S , the HP-41 sets flag F02 and computes  k1 once again & U1'
         the matrix A is replaced by A' in registers Rn+1 thru Rn2+n  and when the program stops again:

     >>>>     k2 is in X-register
     >>>>     U2  is in registers R01 thru Rnn

   ... etc ...
 

Data Registers:

                               R00 =  n
                               R01 = a11          Rn2-n+1 = a1n
                              ...............................................
   INPUTS
                               Rnn = an1             Rn2 = ann

   ---------

                      R00 =  n
                    R01 = u1       Rn+1 =  a'11          Rn2+1 = a'1n

                      .............              A will be replaced with                                           in X-register.   U ( u1 , ........ , un ) = a unit eigenvector
  OUTPUTS                          A' = A - k ( U.tU' ) / ( U.U' )  when F02 is set.

                    Rnn = un         R2n = a'n1           Rn2+n  = a'nn
 

Flag:   F02
Subroutine:  "M*"
 
 
      STACK        INPUTS      OUTPUTS
           X             /        k1 , k2 , ...

 
Example:   Find all the eigenvalues and the corresponding unit eigenvectors of the matrix

             1  2  4
     A =  4  3  5
             7  4  7

1°)    n = 3  therefore     3  STO 00

   Store the 9 numbers

   1  2  4              R01  R04  R07
   4  3  5     in      R02  R05  R08     respectively.
   7  4  7             R03  R06  R09

   XEQ "DFL"                  the HP-41 displays the successive k1-approximations ( with F02 clear )

- Eventually,  k1 = 12.90692995 is in X-register and U1 is in registers R01 R02 R03 :    U1 ( 0.348663346 ; 0.530674468 ; 0.772540278 )

2°) To obtain the second eigenvalue simply press  R/S

     the HP-41 sets flag F02 and displays the successive k1-approximations converging to 12.90692995 again
     then flag F02 is cleared and the HP-41 displays the successive k2-approximations and

          k2 = -2.092097594 in X-register and  U2 in  R01R02R03 :    U2 ( 0.800454175 ; -0.041651079 ; -0.597945065 )

3°) Once again R/S ( similar displays ) and finally,

          k3 =   0.185167649 in X-register and U3 in R01R02R03 :    U3 ( -0.094824730 ; 0.897989404 ; -0.429678136 )

Notes:

-"DFL" uses a random vector as first guess for the eigenvector, so you may get slightly different results than above.
-If  k1 and  k2 are real but  k3 ...  are complex,  the first 2 eigenvalues and eigenvectors will be correctly computed, but not the other ones.
 

Jacobi's Method
 

  "JCB" computes the eigenvalues & eigenvectors of a matrix A, provided all the eigenvalues are real ( and distinct if A is not symmetric ).

*** A non-symmetric real matrix may have complex eigenvalues, but all the eigenvalues of a real symmetric matrix A are necessarily real !

-In the Jacobi's method, the eigenvalues are the elements of a diagonal matrix P equal to the infinite product   ...Ok-1.Ok-1-1....O2-1.O1-1.A.O1.O2....Ok-1.Ok....
 where the Ok are rotation matrices. The eigenvectors are the columns of  O1.O2....Ok-1.Ok....

-"JCB" finds the greatest off-diagonal element ai j  ( lines 22 to 46 )
-Then,  O1  is determined so that it makes a pair of off-diagonal elements zero in O1-1.A.O1
-It happens if the rotation angle µ is choosen so that  Tan 2µ = 2.ai j / ( ai i - aj j )
-The process is repeated until the greatest off-diagonal element is smaller than E-8 in absolute value ( line 48 )

-The successive greatest | ai j |  ( i > j ) are displayed ( line 47 ) when the routine is running.
 

*** Though it's not really orthodox, we can try to apply the same program to non-symmetric matrices:  of course, it cannot work if some eigenvalues are complex.

 But if all the eigenvalues are real, the infinite product   T = [ ti j ] =  ...Ok-1.Ok-1-1....O2-1.O1-1.A.O1.O2....Ok-1.Ok....  is an upper triangular matrix,
 the diagonal elements of T ( i-e  ti i = ki ) are the n eigenvalues and the first column of  U = O1.O2....Ok-1.Ok....   is an eigenvector corresponding to k1

-If all the eigenvalues are distinct, we can create an upper triangular matrix W = [ wi j ] defined by

   wi j = 0    if  i > j
   wi i = 1
   wi j = Sumk = i+1, i+2 , .... , j      ti k wk j / ( tj j - ti i )      if  i < j

-Then, the n columns of  U.W  are the eigenvectors corresponding to the n eigenvalues  t11 = k1 , ............... , tnn = kn

*** When "JCB" stops:

   the n eigenvalues are stored in registers R01  thru Rnn
  the first eigenvector is stored in registers Rn+1 thru R2n
   ..................................................................................
  the nth eigenvector is stored in registers Rn2+1 thru Rn2+n  as shown below.
 

Data Registers:
                               R00 =  n
                               R01 = a11  .................        Rn2-n+1 = a1n
                              ..................................................................              ( the n2 elements of A in column order )
   INPUTS
                               Rnn = an1  .................        Rn2 = ann

   ---------
 

   DURING            R01 thru Rn2 = the "infinite" product   ...Ok-1.Ok-1-1....O2-1.O1-1.A.O1.O2....Ok-1.Ok....
       THE                Rn2+1 thru R2n2 =  O1.O2........Ok.....
 ITERATIONS      R2n2+1 thru R2n2+n: temp  ( for non-symmetric matrices only )

  ---------

                             R00 = n
                             R01 = k1       Rn+1 = V1(1)       R2n+1 = V1(2)  ....................    Rn2+1 = V1(n)
                             .............
  OUTPUTS
                             Rnn = kn        R2n =  Vn(1)         R3n =  Vn(2)    .....................   Rn2+n = Vn(n)
 

                        ( n eigenvalues  ;  1st eigenvector ;  2nd eigenvector ; .................. ;  nth eigenvector )
 

Flag:  F02            Set flag F02  for a symmetric matrix
                             Clear flag F02 for a non-symmetric matrix
Subroutines:  /

     ( SIZE 2n2+1 if A is symmetric / SIZE 2n2+n+1 if A is not symmetric )
 
 
      STACK        INPUTS      OUTPUTS
           X             /            k1

 
Example1:     Let's compute all the eigenvalues and the eigenvectors of the matrix

            1  2  4
    A =  2  7  3
            4  3  9

    n = 3  therefore   3  STO 00  and we store the 9 numbers

     1  2  4               R01  R04  R07
     2  7  3     into    R02  R05  R08     respectively.
     4  3  9               R03  R06  R09

   SF 02  ( A is symmetric )    XEQ "JCB"  yields

                 k1 =  12.81993499  in R01
                 k2 =  4.910741214  in R02
                 k3 = -0.730676199  in R03

and the 3 corresponding eigenvectors:

             V1 (  0.351369026 ; 0.521535689 ;  0.777521917 )   in   ( R04 ; R05 ; R06 )
             V2 ( -0.101146468 ; 0.846760701 ; -0.522269766 )  in   ( R07 ; R08 ; R09 )
             V3 ( -0.930757326 ; 0.104865823 ;  0.350276976 )   in   ( R10 ; R11 ; R12 )

Example2:

             1  2  4
     A =  4  3  5
             7  4  7

     n =  3   STO 00  and we store the 9 numbers

     1  2  4                R01  R04  R07
     4  3  5     into     R02  R05  R08     respectively.
     7  4  7                R03  R06  R09
 

    CF 02  ( A is not symmetric )

    XEQ "JCB"  >>>>    k1 = 12.90692994   and we have

         k1 =  12.90692994  in R01
         k2 =  0.185167649  in R02
         k3 = -2.092097593  in R03

and the 3 corresponding eigenvectors:

          V1 (  0.348663347 ; 0.530674467 ;  0.772540277 )   in   ( R04 ; R05 ; R06 )
          V2 ( -0.095420104 ; 0.903627529 ; -0.432375917 )  in   ( R07 ; R08 ; R09 )
          V3 ( -0.830063031 ; 0.043191757 ;  0.620063097 )   in   ( R10 ; R11 ; R12 )

-The eigenvectors are not unit vectors except the first one.
-Divide V2 & V3 by  || V2 || & || V3 || respectively if you want unit eigenvectors.

Notes:

-If A is non-symmetric, this program only works if all the eigenvalues are real and distinct.
-If all the eigenvalues are real but not distinct, they will be correctly computed but not the eigenvectors ( DATA ERROR line 219 )
 

Jacobi's Method for Complex Matrices
 

-"JCBZ" uses a variant of the Jacobi's algorithm to compute the eigenvalues & eigenvectors of a complex matrix:

*** The eigenvalues of A are the diagonal-elements of an upper triangular matrix T equal to the infinite product   ...Uk-1.Uk-1-1....U2-1.U1-1.A.U1.U2....Uk-1.Uk....
 where the Uk are unitary matrices. The eigenvectors are the columns of  U = U1.U2....Uk-1.Uk.... if A is Hermitian ( i-e if A equals its transconjugate )
 Actually, T is diagonal if A is Hermitian.

-"JCBZ" finds the greatest element ai j below the main diagonal ( lines 23 to 60 )
-Then,  U1  is determined so that it places a zero in position ( i , j )  in U1-1.A.U1

  U1  has the same elements as the Identity matrix, except that  ui i = uj j = x  and  ui j = y + i.z ,  uj i = -y + i.z

     with  y + i.z  = C.x  ,  x = ( 1 + | C |2 ) -1/2  and   C = 2.ai j / [ ai i - ai j + ( ( ai i - ai j )2 + 4 ai j aj i )1/2 ]      ( line 156 R-P avoids a test if the denominator = 0 )

-The process is repeated until the greatest sub-diagonal element is smaller than E-8 in magnitude ( line 62 )

-The successive greatest | ai j |  ( i > j ) are displayed ( line 61 ) when the routine is running.

*** If the matrix is non-Hermitian, and if all the eigenvalues are distinct,  we define an upper triangular matrix W = [ wi j ]  by

   wi j = 0    if  i > j
   wi i = 1
   wi j = Sumk = i+1, i+2 , .... , j      ti k wk j / ( tj j - ti i )      if  i < j

-Then, the n columns of  U.W  are the eigenvectors corresponding to the n eigenvalues  t11 = k1 , ............... , tnn = kn
 

Data Registers:

                               R00 =  n
                               R01,R02 = a11     .................        R2n2-2n+1,R2n2-2n+2 = a1n
                              .................................................................................................                         ( the n2 elements of A in column order )
   INPUTS
                               R2n-1,R2n = an1  .................        R2n2-1,R2n2 = ann

   ---------        odd registers = real part of the coefficients , even registers = imaginary part of the coefficients
 

   DURING            R01 thru R2n2 = the "infinite" product   ...Uk-1.Uk-1-1....U2-1.U1-1.A.U1.U2....Uk-1.Uk....
       THE                R2n2+1 thru R4n2 =  U1.U2........Uk.....
 ITERATIONS      R4n2+1 thru R4n2+2n: temp  ( for non-Hermitian matrices only )

  ---------

                             R00 = n
                             R01,R02 = k1          R2n+1,R2n+2 = V1(1)      ....................    R2n2+1,R2n2+2 = V1(n)
                             .......................            ........................                ....................     ..................................
  OUTPUTS
                             R2n-1,R2n = kn        R4n-1,R4n =    Vn(1)      .....................   R2n2+2n-1,R2n2+2n = Vn(n)
 

                            ( n eigenvalues     ;        1st eigenvector      ;  .................................. ;    nth eigenvector )
 

Flag:  F02          Set flag F02  for a Hermitian matrix
                           Clear flag F02 for a non-Hermitian matrix
Subroutines:  /
 

  ( SIZE 4n2+1 if A is Hermitian / SIZE 4n2+2n+1 if A is non-Hermitian )
 
 
      STACK        INPUTS      OUTPUTS
           Y             /             y1
           X             /            x1

   where  k1 = x1 + y1  = the first eigenvalue.

Example1:

             1      4 -7.i   3-4.i
  A =  4+7.i      6      1-5.i             A is equal to its transconjugate so A is Hermitian
          3+4.i   1+5.i      7

             1  0   4 -7   3 -4           R01 R02    R07 R08    R13 R14
  Store   4  7   6  0   1 -5   into   R03 R04    R09 R10    R15 R16    respectively     and  n = 3  STO 00
             3  4   1  5   7  0            R05 R06    R11 R12    R17 R18

  SF 02  ( the matrix is Hermitian )   XEQ "JCBZ"  yields

           k1 =  15.61385271 + 0.i  = ( X , Y ) = ( R01 , R02 )
           k2 =  5.230678474 + 0.i  = ( R03 , R04 )
           k3 = -6.844531162 + 0.i  = ( R05 , R06 )   and the 3 corresponding eigenvectors

    V1( 0.481585119 + 0.108201123 i , 0.393148652 + 0.527305194 i , -0.142958847 + 0.550739892 i )    =  ( R07 R08 , R09 R10 , R11 R12 )
    V2( -0.436763638 + 0.075843695 i , 0.210271354 - 0.463555612 i , -0.516799072 + 0.526598642 i )   =  ( R13 R14 , R15 R16 , R17 R18 )
    V3( -0.706095475 + 0.247553486 i , 0.326530385 + 0.449069516 i , 0.363126603 - 0. i )                      =  ( R19 R20 , R21 R22 , R23 R24 )

-Due to roundoff errors, registers R02 R04 R06 contain very small numbers instead of 0
  but actually, the eigenvalues of a Hermitian matrix are always real.
 

Example2:

           1+2.i   2+5.i   4+7.i
   A =  4+7.i   3+6.i   3+4.i        A is non-Hermitian
           3+4.i   1+7.i   2+4.i

              1  2   2  5   4  7            R01 R02    R07 R08    R13 R14
   Store   4  7   3  6   3  4   into   R03 R04    R09 R10    R15 R16    respectively     and  n = 3  STO 00
              3  4   1  7   2  4            R05 R06    R11 R12    R17 R18

  CF 02  ( the matrix is non-Hermitian )   XEQ "JCBZ"  produces

           k1 = 7.656606009 + 15.61073835 i  = ( X , Y ) = ( R01 , R02 )
           k2 = 1.661248138 - 1.507335315 i   = ( R03 , R04 )
           k3 = -3.317854151 - 2.103403073 i  = ( R05 , R06 )   and the 3 corresponding eigenvectors

    V1( 0.523491429 + 0.008555748 i , 0.650242685 - 0.046353000 i , 0.546288477 + 0.049882590 i )      =  ( R07 R08 , R09 R10 , R11 R12 )
    V2( -0.4777850326 - 0.196368365 i , 0.660547554 - 0.265888145 i , -0.356842635 + 0.318267519 i )  =  ( R13 R14 , R15 R16 , R17 R18 )
    V3( -0.567221101 - 0.574734664 i , 0.141603481 + 0.549379028 i , 0.480533312 - 0.090337847 i )     =  ( R19 R20 , R21 R22 , R23 R24 )
 

Characteristic Polynomial
 

-If there exist a number k and a non-zero vector V such that  A.V = k.V    k is called an eigenvalue and V an eigenvector of the matrix A.
-The characteristic polynomial   p(x) = xn + c1.xn-1 + c2.xn-2 + ......... + cn-2.x2 + cn-1.x + cn is defined by p(x) = det ( x.I - A )  where I is the n x n unit matrix.
  thus, its roots are actually the n eigenvalues of A.

-"PCAR" uses the following formulae: we define a sequence of matrices Mk by:

      M0 = I   and  Mk+1 = Mk.A - trace(Mk.A) I /k      ( the trace of a matrix is the sum of its diagonal elements )
    then we have:       ck  = - trace(Mk-1.A)/k

-The elements of A are to be stored into contiguous registers from R01 to Rn2 , n being stored in R00
-When the program stops,  R01 = c1 , R02 = c2 , ............... , Rnn = cn            ( c1 = -trace(A)  and cn = (-1)n det(A) )
 

Data Registers:

                               R00 =  n
                               R01 = a11    ......................      Rn2-n+1 = a1n
                              ...................................................................                          ( the n2 elements of A )
   INPUTS
                               Rnn = an1    ......................      Rn2 = ann

   ---------

                          R00 =  n
                          R01 = c1       Rn+1 =  a11          Rn2+1 = a1n
                          .............           ..........................................
  OUTPUTS
                          Rnn = cn         R2n = an1           Rn2+n  = ann
 

Flag:  F10
Subroutine:  "M*"

  ( SIZE 3n2+n+1 )
 
 
      STACK        INPUT      OUTPUT
           X             /         1.n+1

 
Example:        Find the characteristic polynomial of the matrix

             1  2  4
     A =  4  3  5
             7  4  7

-First we store these 9 elements in registers R01 thru R09

       R01  R04  R07         1  2  4
       R02  R05  R08   =    4  3  5    respectively  and  n = 3  STO 00
       R03  R06  R09         7  4  7

-Then  XEQ "PCAR"   >>>>  1.004 = control number of the characteristic polynomial.

   R01 = 1  ,  R02 = -11  ,  R03 = -25  ,  R04 = 5

       Thus p(x) = x3 - 11 x2 - 25 x + 5

Notes:

-A is saved in registers Rn+2 thru Rn2+n+1 and n is saved in R00
-The elements of A may be stored either in column order or in row order since A and tA ( the transpose of A ) have the same characteristic polynomial.

-Any root-finding program yields the 3 eigenvalues    k1 = 12.90692994  ,    k2 = -2.092097589  ,   k3 =  0.185167648
-Then, subtracting ( for instance ) k1 to the diagonal elements of A leads to a singular system with an infinite set of solutions ( the eigenvectors )
-One solution is (  0.451320606 ; 0.686921424 ; 1 )    which is an eigenvector corresponding to the first eigenvalue.

-This method could be used to obtain all the eigenvalues and eigenvectors of a n x n matrix A.
 

Minimal Polynomial
 

-The characteristic polynomial P verifies P(A) = 0 ( Cayley-Hamilton theorem )
-The minimal polynomial p is the unitary least degree polynomial that satisfies  p(A) = 0.
-It is a divisor of the characteristic polynomial.

-We start with a random vector V and we compute  AV , A2V , ........... , AnV
-Then, "PMIN" seeks the smallest k such that  V , AV , A2V , ........... , AkV  are linearly dependent.

-A relation  a0 V + a1 AV + ........... + ak-1 Ak-1V + AkV = 0  is found, whence p(x)
 

Data Registers:                R00 =  n                              ( Registers R00 thru Rn2 are to be initialized before executing "PMIN" )

                                           R01 = a11    ......................      Rn2-n+1 = a1n
                                           ...................................................................                          ( the n2 elements of A )

                                           Rnn = an1    ......................      Rn2 = ann                            Rn2+1 to R2n2+n: temp

          >>>>   When the program stops,  R01 thru Ree = the coefficients of p(x)

Flags:  /
Subroutines:  "M*"  "MCO"   "LS3"

   ( 164 bytes / SIZE 2n2+n+1 )
 
 
      STACK        INPUT      OUTPUT
           X             /         1.eee

 
Example:         Find the minimal polynomial of the matrix

           [ [  3  -1   1 ]
      A = [ -1   3   1 ]
             [  1   1   3 ] ]

-Store
                3   -1    1                     R01  R04  R07
               -1    3    1        into        R02  R05  R08        respectively     and      n = 3   STO 00
                1    1    3                      R03  R06  R09
 

   XEQ "PMIN"  >>>>   1.003     --- Execution time = 36s ---     and we get  R01 = 1   R02 = -5  R03 = 4      ( with some roundoff-errors )

-So,  p(x) = x2 - 5 x + 4       ( the characteristic polynomial is  P(x) = x3 - 9 x2 + 24 x - 16 = ( x - 4 ) p(x)  )
 

Notes:

-Since we use a pseudo-random vector V, there is a small risk of getting an eigenvector, which would lead to a wrong result.
-In practice however, that seems highly improbable.
-If the coefficients of A are integers, the coefficients of p are integers too.
 

Determinant of order 3
 

Data Registers:         R00 = unused                 ( Registers R01 thru R09 are to be initialized before executing "D3" )

                                  R01 = a11     R04 = a12     R07 = a13
                                  R02 = a21     R05 = a22     R08 = a23
                                  R03 = a31     R06 = a32     R09 = a33
Flags: /
Subroutines: /
 
 
      STACK        INPUT      OUTPUT
           X             /     Deteminant

 
Example:      Calculate

            |  4  9  2  |      Store these      R01  R04  R07
    D =  |  3  5  7  |       9 numbers      R02  R05  R08      respectively
            |  8  1  6  |           into            R03  R06  R09

 and   XEQ "D3"  >>>>   Det = 360
 

Determinant of order 4
 

-Determinants of order 4 are computed by a sum of products of determinants of order 2.
 

Data Registers:         R00 = determinant at the end                ( Registers R01 thru R16 are to be initialized before executing "D4" )

                                  R01 = a11     R05 = a12     R09 = a13     R13 = a14
                                  R02 = a21     R06 = a22     R10 = a23     R14 = a24
                                  R03 = a31     R07 = a32     R11 = a33     R15 = a34
                                  R04 = a41     R08 = a42     R12 = a43     R16 = a44
Flags: /
Subroutines: /
 
 
      STACK        INPUT      OUTPUT
           X             /     Deteminant

 
Example:    Calculate

           |   1    8   13  12  |        Store these         R01   R05   R09   R13
   D =  |  14  11   2    7   |       16 numbers         R02   R06   R10   R14      respectively
           |   4    5   16   9   |             into               R03   R07   R11   R15
           |  15  10   3    6   |                                  R04   R08   R12   R16
 

   XEQ "D4"  >>>>  Det =  0
 

Determinant of order 5
 

-Determinants of order 5 are computed by a sum of products of determinants of order 2 by determinants of order 3.
-Registers R00 to R25 are temporarily disturbed during the calculations, but their contents are finally restored.
 

Data Registers:         R00 = det A   at the end                        ( Registers R01 thru R25 are to be initialized before executing "D5" )

                                  R01 = a11     R06 = a12     R11 = a13     R16 = a14     R21 = a15
                                  R02 = a21     R07 = a22     R12 = a23     R17 = a24     R22 = a25
                                  R03 = a31     R08 = a32     R13 = a33     R18 = a34     R23 = a35
                                  R04 = a41     R09 = a42     R14 = a43     R19 = a44     R24 = a45
                                  R05 = a51     R10 = a52     R15 = a53     R20 = a54     R25 = a55
Flags: /
Subroutine:  "D0"
 
 
      STACK        INPUT      OUTPUT
           X             /     Deteminant

 
Example:    Calculate

           |   1    7   13   19   25  |                    R01   R06   R11   R16   R21
           |  14  20  21    2     8   |                    R02   R07   R12   R17   R22
   D =  |  22   3    9    15   16  |          =        R03   R08   R13   R18   R23       respectively
           |  10  11  17   23    4   |                    R04   R09   R14   R19   R24
           |  18  24   5     6    12  |                    R05   R10   R15   R20   R25
 

   XEQ "D5"  >>>>  Det =  -4680000
 

Determinant of order n
 

 "DET"  simply uses "LS3" to compute the determinant of a square matrix of order n.
 

Data Registers:         R00 = det A at the end                          ( Registers R01 thru Rn2 are to be initialized before executing "DET" )

                                  R01 = a11     Rn+1 = a12  .....................    Rn2-n+1 = a1n
                                  R02 = a21     Rn+2 = a22  .....................    Rn2-n+2 = a2n
                                    ........................................................................................

                                  Rnn = an1     R2n = an2    .....................     Rn2 = ann

Flags:    CF 00 = a pivot p is regarded as zero if  | p | < 10-7                 CF 01 = partial pivoting
               SF 00 = a pivot p is regarded as zero if    p = 0                       SF 01 = no pivoting

Subroutines: /
 
 
      STACK        INPUT      OUTPUT
           X             n     Deteminant

   where n is the order of the square matrix.

Example:    Calculate

            |  4  9  2  |                  R01  R04  R07
    D =  |  3  5  7  |      into      R02  R05  R08      respectively
            |  8  1  6  |                  R03  R06  R09

    3   XEQ "DET"  >>>>   Det = 360
 

Linear Systems
 

 "LS" allows you to solve linear systems, including overdetermined and underdetermined systems.
  You can also invert up to a 12x12 matrix.
  It is essentially equivalent to the RREF function of the HP-48 ( a "little" slower, I grant you that! ):
  Its objective is to reduce the matrix on the upper left to a diagonal form with only ones in the diagonal.
  The determinant of this matrix is also computed and stored in register R00.
  ( if there are more rows than columns, R00 is not always equal to the determinant of the upper left matrix because of rows exchanges )

  If flag F01 is clear, Gaussian elimination with partial pivoting is used.
  If flag F01 is set, the pivots are the successive elements of the diagonal:
  It's sometimes useful for matrices like Pascal's matrices of high order.
  They are extremely troublesome and many roundoff errors can occur.
  But if you set flag F01, all the coefficients will be computed with no roundoff error at all, because all the pivots will be ones!

  One advantage of this program is that you can choose the beginning data register - except R00 - ( this feature is used to solve non-linear systems too ):
  You store the first coefficient into Rbb , ... , up to the last one into Ree ( column by column )             ( bb > 00 )

  Then you key in  bbb.eeerr  ENTER^
                                 p          XQ "LS"    and the system will be solved.    ( bbb.eeerr ends up in L-register )

  where r is the number of rows of the matrix
  and  p  is a small non-negative number that you choose to determine tiny elements: if a pivot is smaller than p it will be considered to be zero.

  Don't interrupt "LS" because registers P and Q are used ( there is no risk of crash, but their contents could be altered )
 
 
       STACK         INPUTS         OUTPUTS
            Y         bbb.eeerr               1
            X               p       determinant
            L               /         bbb.eeerr

 
Example1:     Find the solution of the system:

    2x + 3y + 5z + 4t = 39
   -4x + 2y + z + 3t  = 15
    3x  -  y + 2z + 3t = 19
    5x + 7y - 3z + 2t = 18

 If you choose to store the first element into R01, you have to store these 20 numbers:

          2  3  5  4  39                  R01  R05  R09  R13  R17
        -4  2  1  3  15        in        R02  R06  R10  R14  R18      respectively
         3 -1  2  3  19                  R03  R07  R11  R15  R19
         5  7 -3  2  18                  R04  R08  R12  R16  R20

The first register is R01, the last register is R20 and there are 4 rows, therefore the control number of the matrix is  1.02004
If you choose p = 10-7  key in:

                CF 01
  1.02004 ENTER^
      E-7    XEQ "LS"   and  you obtain  840 in X-register and in R00:  it is the determinant of the 4x4 matrix on the left.

Registers R01 thru R16 now contains the unit matrix and registers R17 thru R20 contains the solution  x = 1 ; y = 2 ; z = 3 ; t = 4
Thus, the array has been changed into:

         1  0  0  0  1
         0  1  0 0  2              the solution is the last column
         0  0  1  0  3              of the new matrix.
         0  0  0  1  4

Example2:    Solve the system:

          5x + y + z  = 8
          4x - y + 2z = 13
           x + 2y - z  = -5
         7x - 4y + 5z = 31

-Suppose we need to store the first element into  R11 , then we store these 16 numbers

          5   1  1   8                       R11  R15  R19  R23
          4  -1  2  13         in          R12  R16  R20  R24          respectively
          1   2  -1  -5                     R13  R17  R21  R25
          7  -4  5  31                      R14  R18  R22  R26

Here, the control number of this array is   11.02604   and if we take p = 10-7 again,

              11.02604  ENTER^
                    E-7     XEQ "LS"   gives  0  in X-register and in R00 = the determinant of the 4x4 matrix.

and registers R11 thru R27 are now:

        1  0   0.333333333     2.333333333
        0  1  -0.666666667   -3.666666669
        0  0      5.10-10              -2.10-9
        0  0          0                    4.10-9

Thus, the system is equivalent to:

      x + z /3  =  7/3              and there is an infinite set of solutions:
      y - 2z /3 = -11/3           z may be chosen at random and x and y are determined by
                0  =  0                 x = -z /3 + 7/3
                0  =  0                 y = 2z /3 - 11/3
 

Note:    If the last number of the initial matrix is 32 instead of 31 the array is changed into:

       1  0   0.333333333   0        i-e         x + z /3 = 0
       0  1  -0.666666667   0                     y - 2z /3 = 0
       0  0        5.10-10        0                                 0 = 0
       0  0           0              1                                 0 = 1        and there is no solution at all !

Example3:   Invert Pascal's matrix:

      1  1  1   1   1
      1  2  3   4   5
      1  3  6  10 15
      1  4 10 20 35
      1  5 15 35 70

This is equivalent to solve 5 linear systems simultaneously. If you begin at register R01, you have to store the 50 numbers

        1  1  1   1  1    1  0  0  0  0                  R01  R06  R11  R16  R21  R26  R31  R36  R41  R46
        1  2  3   4  5    0  1  0  0  0                  R02  R07  R12  R17  R22  R27  R32  R37  R42  R47
        1  3  6  10 15  0  0  1  0  0         in      R03  R08  R13  R18  R23  R28  R33  R38  R43  R48        ( the right half is the unit matrix )
        1  4 10 20 35  0  0  0  1  0                  R04  R09  R14  R19  R24  R29  R34  R39  R44  R49
        1  5 15 35 70  0  0  0  0  1                  R05  R10  R15  R20  R25  R30  R35  R40  R45  R50

The control number of this rectangular array is 1.05005 and we can choose p = 0 because Pascal's matrix is non-singular:

    1.05005  ENTER^
           0      XEQ "LS"  yields  1 in X-register and in R00  ( the determinant of Pascal's matrices is always 1 )

 The array is now:

       1  0  0  0  0   5  -10  10  -5   1
       0  1  0  0  0 -10  30 -35  19 -4
       0  0  1  0  0  10 -35  46 -27  6           ( the unit matrix is now on the left and the inverse matrix is on the right,
       0  0  0  1  0  -5   19 -27  17 -4              in registers R21 thru R50 )
       0  0  0  0  1   1   -4    6   -4   1
 

-This program can be used to invert a n x n matrix, but it requires 2n2 registers and thus, it cannot invert a 13x13 matrix.
-Use "MINV" if you want to find the inverse a 17x17 matrix.
 

Linear Systems ( variant )
 

-In this variant, the matrix is the right-part of the array so that, when the program stops,
  the solution ( x1 , x2 , ............. , xn ) is in registers  R01  R02  ..............  Rnn
-Here, the attempt to diagonalization starts by the lower right corner of the matrix.
-Flags F00 & F01 play the same role as above.
 

   ( SIZE  n.m+1 )
 
 
       STACK         INPUTS         OUTPUTS
            T              /           n.nnn
            Z              /              /
            Y              n              /
            X              m          det A

     n  =  number of rows
     m = number of columns

 T-output is useful to retrieve n

Example:

          2x + 3y + 5z + 4t = 39                                       39 = 2x + 3y + 5z + 4t
         -4x + 2y + z + 3t  = 15           is re-written          15 = -4x + 2y + z + 3t
          3x  -  y + 2z + 3t = 19                                       19 = 3x  -  y + 2z + 3t
          5x + 7y - 3z + 2t = 18                                       18 = 5x + 7y - 3z + 2t

  and we store these 20 numbers:

         39    2  3  5  4                        R01  R05  R09  R13  R17
         15  -4  2  1  3         in            R02  R06  R10  R14  R18      respectively
         19   3 -1  2  3                        R03  R07  R11  R15  R19
         18   5  7 -3  2                        R04  R08  R12  R16  R20

-There are 4 rows and 5 columns,

      CF 00   CF 01
      4  ENTER^
      5  XEQ "LS3"  >>>>  Det A = 840 = R00

-Registers R05 thru R16 now contain the unit matrix and registers R01 thru R04 contain the solution  x = 1 , y = 2 , z = 3 , t = 4
-Thus, the array has been changed into:

           1  1  0  0  0
           2  0  1  0  0                the solution is the first column
           3  0  0  1  0                of the new matrix.
           4  0  0  0  1

-When the program stops, R00 = det A
-If you have to invert a matrix, the inverse will be in registers R01 thru Rn2 at the end
 

Underdetermined Linear Systems
 

-"LS3-" computes the Moore-Penrose solution given by the formula:   X = tA ( A.tA ) -1 B

   ( SIZE 2n.m-n+1 )
 
 
       STACK         INPUTS       OUTPUTS
            Y              n              /
            X              m         1.rrr,rr

    n = number of rows
   m = number of columns      ( n < m-1 )
   1.rrr,rr  is the control number of the solution ( in fact r = m-1 )
 

Example:    Let's find the Moore-Penrose solution of the system:

      1 = 2x + 3y + 7z + 4t                 1  2  3   7   4           R01  R04  R07  R10  R13
      4 = 3x + 2y - 5z + 8t      store    4  3  2  -5  8   into   R02  R05  R08  R11  R14   respectively
      7 = 4x + 5y + 6z + t                   7  4  5   6   1           R03  R06  R09  R12  R15

-There are 3 rows and 5 columns so:

      3  ENTER^
      5  XEQ "LS3-"  >>>>   1.00404

-The solution is in registers   R01 R02 R03 R04  namely:    x = 1.095088162 , y = 1.075776659 , z = -0.389378674 , t = -0.422963896

-When the program stops, R00 = det ( A.tA )   [ in this example, R00 = 128628 ]
-If R00 = 0 ( suggesting A.tA is singular), the results are probably meaningless.
 

Overdetermined Linear Systems
 

-An overdetermined system  A.X = B  has often no solution at all!
-But "LS3+" calculates the least-squares solution:
-It minimizes || A.X - B ||

Formula:   X = ( tA.A ) -1 ( tA.B )
 

  ( SIZE n.m+m2 -m+1 )
 
 
       STACK         INPUTS         OUTPUTS
            Y              n             det
            X              m            1.rrr

     n = number of rows
     m = number of columns      ( n >= m )
    1.rrr  is the control number of the solution ( in fact r = m-1 )

Example:      The system:

       8  =  5x + y + z              still has no solution.
      13 =  4x - y + 2z
      -5 =   x + 2y - z
      32 = 7x - 4y + 5z
    -20 =  2x + 5y - 9z

-We store these 20 coefficients into registers R01 thru R20   ( in column order )

          8    5   1   1                     R01  R06  R11  R16
         13   4  -1  2                     R02  R07  R12  R17
         -5   1   2  -1         =         R03  R08  R13  R18          respectively
         32   7  -4  5                     R04  R09  R14  R19
        -20   2   5  -9                    R05  R10  R15  R20

-There are 5 rows and 4 columns so

  5  ENTER^
  4  XEQ "LS3+"  >>>>   1.003

-The "least-squares solution" is in registers    R01 R02 R03  namely:   x = 2.071207430  ,  y = -3.201238398  ,  z = 0.904024763

-When the program stops, Y = R00 = det ( tA.A )   [ In this example  R00 = 55233 ]
-If det ( tA.A )  = 0  ( suggesting tA.A is singular ) , the results are probably meaningless.
 

4-Dimensional Magic Hypercubes of odd orders
 

-"MG4DHC" creates 4-D magic (hyper-) cubes of odd orders  n
-The magic constant =  n ( n4+1 )/2
 

Data Registers:       R00 = n        R01 thru R06:  temp

Flag:  F00   SF 00 to display the successive elements
                    CF 00 to store the elements into R01 to Rn^4  ( possible only if n = 1 or 3 )

Subroutines:  /
 
 
      STACK        INPUTS      OUTPUTS
           X             n             /

  Where n is an odd integer > 0

Example:     SF 00    3   XEQ "MG4DHC"   >>>>    N1=46   N2=59   N3=18   ...................   N79=64   N80=23   N81=36

The whole hypercube is:

   46   59   18        62   12   49        15   52   56
   17   48   58        51   61   11        55   14   54                   1st layer = cube n°1
   60   16   47        10   50   63        53   57   13

   08   39   76        42   79   02        73   05   45
   78   07   38        01   41   81        44   75   04                   2nd layer = cube n°2
   37   77   09        80   03   40        06   43   74

   69   25   29        19   32   72        35   66   22
   28   68   27        71   21   31        24   34   65                   3rd layer = cube n°3
   26   30   67        33   70   20        64   23   36

-Here the magic constant is 123 and the central element is 41.

-There are 8 great diagonals:

   36 + 41 + 46 = 123           67 + 41 + 15 = 123
   22 + 41 + 60 = 123           29 + 41 + 53 = 123
   35 + 41 + 47 = 123           26 + 41 + 56 = 123
   64 + 41 + 18 = 123           69 + 41 + 13 = 123
 

Magic Cubes
 

-"MGCB" creates normal ( Andrews ) magic cubes, i-e their elements are  1  2  3  ...............  n3
-An Andrews magic cube is a cubic array where the sums of the elements equal the magic constant n(n3+1)/2  in the 3 directions and the 4 space diagonals.
 

Data Registers:        R00 = n

                 >>>     R01 thru Rn3 = the elements of the cube  if flag F00 is clear

Flag:   F00

    SF 00 = The elements are successively displayed
    CF 00 = The elements are stored into R01  R02  .............

Subroutines: /
 
 
      STACK        INPUTS      OUTPUTS
           X          n # 2             /

 
Example:      With   n = 3

-SF 00   3   XEQ "MGCB"  the HP-41 displays   N1=1   N2=17   N3=24   .................   N25=4   N26=11   N27=27

-If  CF 00, we get in registers R01 thru R27

   R01   R04   R07       01   23   18
   R02   R05   R08  =   17   03   22
   R03   R06   R09       24   16   02

             R10   R13   R16       15   07   20
             R11   R14   R17  =   19   14   09
             R12   R15   R18       08   21   13

                                R19   R22   R25      26   12   04
                                R20   R23   R26  =  06   25   11
                                R21   R24   R27      10   05   27

-The magic constant = 42
 

Magic Squares
 

-"MGSQ"  creates normal magic squares of order n # 2

-The elements of a normal magic square are  1 , 2 , ............. , n2  and it has the following property:
 All its row sums, column sums and both diagonal sums are equal to the magic constant   M = n.(n2 + 1)/2
 

Data Registers:        R00 = n

                 >>>     R01 thru Rn2 = the elements of the square  if flag F00 is clear

Flags:   F00

    SF 00 = The elements are successively displayed
    CF 00 = The elements are stored into R01  R02  .............

Subroutines: /
 
 
      STACK        INPUTS      OUTPUTS
           X          n # 2             /

 
Example:     n = 4

-SF 00   4   XEQ "MGSQ"   the HP-41 displays successively   N1=16   N2=5  N3=9  N4=4  ....................  N13=13  N14=8  N15=12  N16=1

-The whole magic square is:

   16   02   03   13
   05   11   10   08
   09   07   06   12
   04   14   15   01

-Here, the magic constant = 34
 

Pandiagonal & Panmagic Cubes
 

-"PMGCB"  creates a perfect pandiagonal magic cube if n = 7 & panmagic cubes if n is a prime > 10  ( look at line 02 )

-The cube is perfect pandiagonal if all the orthogonal sections are panmagic squares.
-In a panmagic cube , all the orthogonal or diagonal sections - broken or unbroken - are panmagic squares.
 

Data Registers:        R00 = n

                 >>>     R01 thru Rn3 = the elements of the cube  if flag F00 is clear

Flag:   F00

    SF 00 = The elements are successively displayed
    CF 00 = The elements are stored into R01  R02  .............

Subroutines: /
 
 
      STACK        INPUTS      OUTPUTS
           X       n prime > 6             /

 
Example:     n = 7

  SF 00    7     XEQ "PMGCB"   >>>>   N1=1   N2=334   N3=268   N4=209   N5=192   N6=133   N7=67   ...............   N343=343

Note:

-Since n > 6 , only the first elements will be stored if  CF 00
 

Panmagic Squares
 

-"PMGSQ"  creates panmagic squares of order n  where n mod 4 # 2 & n mod 3 # 0

-A panmagic square ( also called pandiagonal square ) is a magic square
 where all the broken diagonal sums are also equal to the magic constant  M = n.(n2 + 1)/2
 

Data Registers:        R00 = n

                       >>>     R01 thru Rn2 = the elements of the square  if flag F00 is clear

Flag:   F00

    SF 00 = The elements are successively displayed
    CF 00 = The elements are stored into R01  R02  .............

Subroutines: /
 
 
      STACK        INPUTS      OUTPUTS
           X             n             /

   where n mod 4 # 2 & n mod 3 # 0

Example:    n = 5

-SF 00   5   XEQ "PMGSQ"   the HP-41 displays   N1=1   N2=7   N3=13   N4=19   N5=25  ...................  N25=12

-The panmagic square is

   01   14   22   10   18
   07   20   03   11   24
   13   21   09   17   05
   19   02   15   23   06
   25   08   16   04   12

-Magic Constant = 65

Note:

-There is no normal panmagic square of order 2 , 6 , 10 , ....
 

Bairstow's Method
 

-"BRST" factorizes the polynomial   p(x) = an.xn+an-1.xn-1+ ... + a1.x+a0  into quadratic factors and solves  p(x) = 0      ( n > 1 )
-If deg(p) is odd, we have   p(x) =    (an.x+b).(x2+u1.x+v1)...............(x2+um.x+vm)      where m = (n-1)/2
-If deg(p) is even                p(x) = (anx2+u1.x+v1)(x2+u2.x+v2)..........(x2+um.x+vm)      where m  =  n/2

-The coefficients u and v are found by the Newton method for solving 2 simultaneous equations.
-Then  p is divided by (x2+u.x+v)  and  u & v  are stored into Ree-1 & Ree respectively

-The process is repeated until all quadratic factors are found, and when the program stops:

  If deg(p) is odd          Rbb = an , Rbb+1 = b , Rbb+2 = u1 , Rbb+3 = v1 , ........ , Ree-1 = um , Ree = vm
  If deg(p) is even         Rbb = an , Rbb+1 = u1 , Rbb+2 = v1 , ............................ , Ree-1 = um , Ree = vm

-Then, press R/S to get the successive roots.
 

Data Registers:  R00 to R08: temp.   ( R06 = u ; R07 = v ; R08 = bbb.eee )

                           Rbb = an    Rbb+1 = an-1 , ....... ,   Ree = a0       ( These n+1 registers are to be initialized before executing "BRST")    bb > 08
Flag:  F00
Subroutines:  "DIV" & "P2"
 
 
 
      STACK         INPUT     OUTPUT0     OUTPUTS1  .........................    OUTPUTSk
           Y             /            /           y1  .........................           yk
           X         bbb.eee       bbb.eee           x1  .........................           xk

  where  bbb.eee  is the control number of the polynomial ( bbb > 008 )

-When the program stops,  Rbb thru Ree contain the factorization
-Then, the successive outputs are the different pairs of roots:     xk  ,  yk   if flag F00 is clear  or   xk+i.yk  ,   xk-i.yk  if F00 is set

 However , when deg(p) is odd , the first root is always real:  F00 is set but y = 0

-The last pair of roots is indicated by a BEEP, the others by a TONE 9.

Example:     Find all the roots of   2.x5+7.x4+20.x3+81.x2+190.x+150

-For instance,  2 STO 11   7 STO 12   20 STO 13   81 STO 14   190 STO 15   150 STO 16   ( control number = 11.016 ) and if we choose  u0 = v0 = 0

  FIX 8    11.016  XEQ "BRST"  the successive uk are displayed and sometime later, we hear a BEEP and X = 11.016

  R11 = 2   R12 = 3     //     R13 = -2   R14 = 10     //    R15 = 4   R16 = 5            whence   p(x) = (2x+3).(x2-2x+10).(x2+4x+5)    then:

       R/S  >>> ( TONE 9 )  -1.5    X<>Y   0   with  F00 set     the first root is   -1.5    ( deg(p) is odd , the 1st root is real )
       R/S  >>> ( TONE 9 )    1      X<>Y   3    with  F00 set     1+3.i  and  1-3.i
       R/S  >>>   ( BEEP )     -2      X<>Y   1    with  F00 set    -2+i  and  -2-i           the 5 roots are  -1.5 ; 1+3.i ; 1-3.i ; -2+i ; -2-i
 

Notes:

-"BRST" always choose  u0 = 2 and  v0 = PI
-The accuracy is controlled by the display setting: the loop stops when the rounded increment equals 0

-If you get "DATA ERROR" or if the process seems to diverge, stop the program, change registers R06 & R07 and press  XEQ 01

-If you want to see the roots again, press XEQ 10
-If the displayed u-values seem to oscillate between closed numbers without stopping the loop - it may happen with multiple roots -
  Stop the program, set F21 , press  R/S.  "BRST" will stop at the next VIEW 06.
  Clear X-register an R/S
 

Cohn's Irreducibility Criterion
 

-This test concerns a poynomial  p(x) = an xn + an-1 xn-1 + ... + a1 x + a0  in Z[X]   with  ak >= 0  for all k's  ( a0 # 0 )

    If there exists an integer m > max { ak }  such that  f(m) is a prime, then p(x) is irreducible in Z[X]
 

Data Registers:              R00 to R02:  are used by "PR?"            ( Registers Rbb thru Ree are to be initialized before executing "COHN" )
                                         R03 = bbb.eee    R04 = m

                                        Rbb = an    Rbb+1 = an-1  ........................    Ree = a0 # 0            bbb > 004
Flags: /
Subroutines:   "PL"    "PR?"
 
 
      STACK        INPUT      OUTPUT
           X        bbb.eee            1

   where   bbb.eee = control number of p(x)            X-output = 1    if Cohn's criterion is satisfied,                 bbb > 004

Example1:       p(x) = 2 x4 + 3 x2 + 5 x + 1

  With  [ 2  0  3  5  1 ]  =  [ R05  R06  R07  R08  R09 ]

  5.009  XEQ "COHN"  >>>>   1         p(x) is irreducible in Z[X]

  R04 = 6  and  p(6) = 2731 is a prime
 

Example2:       p(x) = 2 x6 +  3 x5 + 4 x3 + 5 x2 + 2 x + 1

  With  [ 2  3  0  4  5  2  1 ]  =  [ R05  R06  R07  R08  R09  R10  R11 ]

  5.011  XEQ "COHN"  >>>>   1         p(x) is irreducible in Z[X]

  R04 = 8  and  p(8) = 624977 is a prime
 

Note:

-If no prime is found, the routine continues until f(m) exceeds 1010 and the HP-41 displays "OUT OF RANGE" and in this case, we cannot conclude.
 

Polynomial Euclidean Division
 

-If   a(x) = an.xn+an-1.xn-1+ ... + a1.x+a0  and   b(x) = bm.xm+bm-1.xm-1+ ... + b1.x+b0  are 2 polynomials,
  there are only 2 polynomials q(x) and r(x)  such that   a = b.q + r  with  deg(r) < deg(b)

-"DIV" calculates q and r.
 

Data Registers:  R00 to R03: temp                     ( Registers Rbb thru Ree & Rbb' thru Ree' are to be initialized before executing "DIV" )

                            Rbb = an     Rbb+1 = an-1  , ....... ,    Ree = a0        bb > 03
                            Rbb' = bm    Rbb'+1 = bm-1 , ....... ,   Ree' = b0       bb' > 03
Flags: /
Subroutines: /
 
 
 
      STACK        INPUTS      OUTPUTS
           Y       (bbb.eee)a       (bbb.eee)r
           X       (bbb.eee)b       (bbb.eee)q

where         (bbb.eee)a = the control number of the dividend                        (bbb.eee)r  = the control number of the remainder        ( all  bbb > 003 )
                  (bbb.eee)b = the control number of the divisor                           (bbb.eee)q = the control number of the quotient

Example:     a(x) =  2.x5+5.x4-21.x3+23.x2+3.x+5      b(x) = 2.x2-3.x+1     Find   q(x) and r(x)

For instance:        2 STO 04   5 STO 05   -21 STO 06   23 STO 07   3 STO 08   5 STO 09     ( control number = 4.009 )
             and        2  STO 11   -3 STO 12   1 STO 13                                                                ( control number = 11.013 )

     4.009  ENTER^
   11.013  XEQ "DIV"  >>>>     4.007  X<>Y  8.009        ( q and r have replaced a  ;  b is unchanged )

 R04 = 1   R05 = 4   R06 = -5   R07 = 2   whence  q(x) = x3+4.x2-5.x+2
                R08 = 14   R09 = 3                   whence  r(x) = 14.x + 3

Notes:

-When b(x) is constant, (bbb.eee)r = 1+eee.eee
-Roundoff errors may produce small leading coefficients instead of zero in the remainder.

       "DIV" doesn't work if  deg(a) < deg(b)    but in this case,  q = 0 and r = a
 
 

Factorization over Z[X]
 

-"FCTR" uses "BRST" to factorize in R[X] a polynomial p(x)  = an.xn+an-1.xn-1+ ... + a1.x+a0  whose coefficients are integers
-The products of their constant coefficients are tested to see if we get integers after rounding the results.

-Then, p(x) is divided by the possible factors to check if the remainder really equals zero.

-So, the display setting plays a central role.

-To simplify the program, we assume that the leading coefficient  an = +/-1  though it may sometimes work in other cases too.
 

Data Registers:       R00 to R12+4n: temp                                           ( Registers R01 thru Rnn+1 are to be initialized before executing "FCTR" )

                                 R04 = counter   R05 = eee+1   R06 = address of the first free register
                                 R07 = control number of the list of the irreducible factors in R[X]
                                 R08 = control number of a potential divisor
                                 R09 = control number of f(x)
                                 R10 = counter = 2m-1 - 1 , 2m-1 - 2 ,  .......  , 0

                   R01 = an    R02 = an-1  ........................    Rnn+1 = a0          and several other registers to store the products of polynomials.

Flags:  F00-F01-F26
Subroutines:   "BRST"   "P2"  "¨PPRD"    "MCO"
 
 
      STACK        INPUT      OUTPUTk
           Y         1.n+1            /
           X            m        bbb.eeek

  where   1.n+1 = control number of p(x)                                      bbb.eeek = control number of the successive factors.
    and     m is used to round numbers in FIX m                            with  n = deg p

Example:          p(x) = x8 - 128 x6 - x5 + 5120 x4 + 64 x3 - 65538 x2 - 512 x + 262144
 

  With  [ R01  R02  R03  R04  R05  R06  R07  R08  R09 ]  =  [ 1  0  -128  -1  5120  64  -65538  -512  262144 ]

 >>>  FIX 8  if you want to execute "BRST" in this format, and if you want to execute the rest of "FCTR" in FIX 4

  1.009  ENTER^
     4      XEQ "FCTR"  >>>>   ( TONE 9 )    60.064

  [ R60  R61  R62  R63  R64 ]  =  [ 1  0  -64  -2  512 ]

               R/S    >>>>    ( BEEP )   44.048     ( in a few seconds )

  [ R44  R45  R46  R47  R48 ]  =  [ 1  0  -64   1   512 ]

-And    p(x) = ( x4 - 64 x2 - 2 x + 512 )( x4 - 64 x2 + x + 512 )

Notes:

-With   1.009   ENTER^   5    XEQ "FCTR"   NO factorization would have been found !
-Before rounding the coefficients, the product was  [ 1  0.000000999  -63.99999993  -2.000006500  511.9999998 ]
  and  -2.000006500  would have been rounded to  -2.00001  which is not an integer!

-If deg f is large and/or if there are many factors in R[X], we'll have to execute"FCTR" with m = 1
 and even then, we could miss some factors because of roundoff errors...
-So the result is not completely guaranteed. Use "COHN" and/or "IRR?" to check that the factors are really irreducible...

-If the process seems to diverge while "BRST" is executing, stop the subroutine, store other guesses in R06 & R07 and press GTO 01  R/S
-But do not press  XEQ 01 :  that would clear the return stack !

-If the leading coefficient  c = an # +/-1 , "FCTR" can fail.
-Apply the program to h(x) = cn-1 p(x/c)
-After factorizing h(x), the factorization of p(x) will be easy to deduce.
 

A Better Irreducibility Criterion
 

-Unlike Cohn's criterion, the following one also works if the polynomial has negative coefficients.
-Let  p(x) = an xn + an-1 xn-1 + ... + a1 x + a0  in Z[X]  deg p > 1    ( a0 # 0 )

      m1 = max { | ai/an | ,  i = 0 , 1 , .... , n-1 }
      m2 = max { | ai/an |1/(n-i) ,  i = 0 , 1 , .... , n-1 }

      H = min { 1+m11/(r+1) , 2 m2 }   where  r is such that  an-1 = an-2 = ..... = an-r = 0 ,     ( r = 0 if an-1 # 0 )

    If there exists an integer m > H  such that  p(m) is a prime or +/-1, then p(x) is irreducible in Z[X]
 

Data Registers:              R00 to R02:  temp                               ( Registers Rbb thru Ree are to be initialized before executing "IRR?" )
                                         R03 = bbb.eee    R04 = m

                                       Rbb = an    Rbb+1 = an-1  ........................    Ree = a0 # 0                bbb > 004
Flags: /
Subroutines:   "PL"   "PR?"
 
 
      STACK        INPUT      OUTPUT
           X        bbb.eee   1 or 0 or E10

  where   bbb.eee = control number of p(x)                                   bbb > 004

     if  X-output = 1   p is irreducible
     if  X-output = 0   p is reducible: we've found an integer root
     if  X-output = E10 the HP-41 displays "OUT OF RANGE" and we cannot conclude

Example:       f(x) =  2 x7 -  3 x4 + 4 x3 - 9 x2 + 6 x - 1

  With  [ 2  0  0  -3  4  -9  6  -1 ]  =  [ R05  R06  R07  R08  R09  R10  R11  R12 ]

   5.012  XEQ "IRR?"   >>>>   1    thus,  f is irreducible

   R04 = m = 6  and  f(6) = 556559  is a prime
 

Notes:

-In some cases - when all the coefficients are non-negative and an is small - Cohn's criterion may be better ( for instance f(x) = x2 + 1 )
 

Quadratic Equations
 

-"P2" solves the equation:   a.x2+b.x+c = 0    with a # 0

Data Registers:  /
Flags:  F00  ( indicates complex roots )
Subroutines:  /
 
 
 
      STACK        INPUTS      OUTPUTS
           Z          a # 0           c/a
           Y             b             y
           X             c             x

-If  CF 00  the 2 solutions are  x , y
-If  SF 00   ------------------  x+i.y , x-i.y

Examples:     Solve  2.x2 + 3.x - 4 = 0   and   2.x2 + 3.x + 4 = 0

             2  ENTER^
             3  ENTER^
           -4  XEQ "P2"  >>>>   -2.350781059
                    X<>Y                   0.850781060      Flag  F00 is clear, therefore the 2 solutions are  -2.350781059  and  0.850781060

            2  ENTER^
            3  ENTER^
            4     R/S         >>>>   -0.75
                 X<>Y                    1.198957881      Flag  F00 is set, therefore the 2 solutions are  -0.75 + 1.198957881.i  and   -0.75 - 1.198957881.i
 

Cubic Equations
 

-"P3" finds the 3 roots of   a.x3+b.x2+c.x+d  by the Cardan's ( or Tartaglia's ) formulae:  ( with a # 0 )
-First, the term in x2 is removed by a change of argument, leading to  x3+p.x+q = 0
-Then, x = u+v  with u.v = -p/3 leads to a quadratic equation in u3
 

Data Registers:  /
Flags:  F01             ( indicates complex roots )
Subroutine:  none  if d # 0 , "P2" if d = 0
 
 
 
      STACK        INPUTS      OUTPUTS
           T          a # 0             /
           Z             b             z
           Y             c             y
           X             d             x

-If  CF 01  the 3 solutions are  x ; y ; z
-If  SF 01  the 3 solutions are  x ; y+i.z ; y-i.z

Example:      Solve  2.x3-5.x2-7.x+1 = 0  and   2.x3-5.x2+7.x+1 = 0

        2   ENTER^
       -5  ENTER^
       -7  ENTER^
        1   XEQ "P3"  >>>>   3.467727065   RDN   0.131206041   RDN   -1.098933107   which are the 3 real solutions since flag F01 is clear.

        2  ENTER^
       -5  ENTER^
        7   ENTER^
        1     R/S         >>>>   -0.130131632   RDN   1.315065816   RDN   1.453569820

-Flag F01 is set , therefore the 3 solutions are  -0.130131633  ;  1.315065816 - 1.453569820.i  ;  1.315065816 + 1.453569820.i
 

Quartic Equations
 

-"P4" solves the equation   x4+a.x3+b.x2+c.x+d = 0   ( if the leading coefficient is not 1 , divide all the equation by this coefficient ).
-First, the term in x3 is removed by a change of argument, leading to  x4+p.x2+q.x+r = 0  (E')
-Then, the resolvant  z3+p.z2/2+(p2-4r).z/16-q2/64 = 0 is solved by "P3" and if we call  z1 , z2 , z3  the 3 roots of this equation, the zeros of (E') are

          x = z11/2 sign(-q) +/- ( z21/2+z31/2 )    ;    x = -(z11/2) sign(-q) +/- ( z21/2-z31/2 )

Data Registers:  R00 thru R04: temp
Flags:   F01 F02
Subroutine:    "P3"
 
 
 
      STACK        INPUTS      OUTPUTS
           T             a             t
           Z             b             z
           Y             c             y
           X             d             x

-If  CF 01  &  CF 02  the 4 solutions are    x ; y ; z ; t
-If  SF 01  &  CF 02   ------------------    x+i.y ; x-i.y ; z ; t
-If  SF 01  &  SF 02   ------------------    x+i.y ; x-i.y ; z+i.t ; z-i.t

Example1:    Solve  x4-2.x3-35.x2+36.x+180 = 0

       -2   ENTER^
      -35  ENTER^
       36   ENTER^
      180  XEQ "P4"  >>>>   6  RDN   3   RDN   -2   RDN   -5    Flags F01 & F02 are clear , the 4 solutions are   6 ; 3 ; -2 ; -5

Example2:    Solve  x4-5.x3+11.x2-189.x+522 = 0

        -5   ENTER^
        11   ENTER^
     -189   ENTER^
       522      R/S     >>>>   -2  RDN    5   RDN   3   RDN   6  Flag F01 is set & F02 is clear , the 4 solutions are  -2+5.i ; -2-5.i ; 3 ; 6

Example3:    Solve  x4-8.x3+26.x2-168.x+1305 = 0

       -8   ENTER^
       26   ENTER^
     -168  ENTER^
     1305    R/S      >>>>   -2  RDN    5   RDN   6   RDN   3    Flags F01 & F02 are set , the 4 solutions are  -2+5.i ; -2-5.i ; 6+3.i ; 6-3.i
 

Evaluation of a Polynomial P(x)
 

-"PL"  calculates  p(x) = an.xn+an-1.xn-1+ ... + a1.x+a0  for any given x-value.
 

Data Registers:        Rbb = an    Rbb+1 = an-1 , ....... ,   Ree = a0         ( These n+1 registers are to be initialized before executing "PL" )
Flags: /
Subroutines:  /
 
 
 
      STACK        INPUTS      OUTPUTS
           Y        bbb.eee             /
           X             x           p(x)
           L             /             x

 
Example:     p(x) = 2.x3+3.x2-4.x+7     Find  p(5)

-If we store these 4 coefficients into R01 to R04  ( 2  STO 01  3  STO 02  -4  STO 03  7  STO 04 )

   1.004  ENTER^
       5     XEQ "PL"  >>>>  p(5) = 312
 

Product of 2 Polynomials
 

-Let  a(x) = an.xn+an-1.xn-1+ ... + a1.x+a0  and   b(x) = bm.xm+bm-1.xm-1+ ... + b1.x+b0   2 polynomials,

 "PPRD" computes and stores the coefficients of  p(x) = a(x).b(x) into contiguous registers. ( you have to specify the first of these registers )
 

Data Registers:  R00 to R03: temp                             ( Registers Rbb thru Ree & Rbb' thru Ree' are to be initialized before executing "PPRD" )

                            Rbb = an     Rbb+1 = an-1  , ....... ,    Ree = a0         bb > 03
                            Rbb' = bm    Rbb'+1 = bm-1 , ....... ,   Ree' = b0        bb' > 03
Flags: /
Subroutines: /
 
 
 
      STACK        INPUTS      OUTPUTS
           Z       (bbb.eee)a             /
           Y       (bbb.eee)b             /
           X            bbbp       (bbb.eee)p

  where       (bbb.eee)a = the control number of  a(x)
    and         (bbb.eee)b = the control number of  b(x)                           (bbb.eee)p = the control number of the product        ( all  bbb > 003 )

Example:     a(x) =  2.x5+5.x4-21.x3+23.x2+3.x+5      b(x) = 2.x2-3.x+1         Find  p(x) = a(x).b(x)

For instance:        2 STO 04   5 STO 05   -21 STO 06   23 STO 07   3 STO 08   5 STO 09     ( control number = 4.009 )
                           2  STO 11   -3 STO 12   1 STO 13                                                                ( control number = 11.013 )

  and if we want to have p(x) in registers R21  R22  ... etc ...

     4.009  ENTER^
   11.013  ENTER^
       21     XEQ "PPRD"  yields  21.028

-We have   R21 = 4   R22 = 4   R23 = -55   R24 = 114   R25 = -84   R26 = 24   R27 = -12   R28 = 5

     whence  p(x) = 4.x7+4.x6-55.x5+114.x4-84.x3+24.x2-12.x+5
 

Note:    Neither (bbb.eee)a  &   (bbb.eee)p   nor   (bbb.eee)b  &   (bbb.eee)p   can overlap.
 

Primality Test ( Non-negative integers )
 

-"PR?" takes a non-negative integer n in X-register,
  and returns 1 if n is a prime or if n = 1   and 0 otherwise
-The smallest divisor d is in L-register.
 

Data Registers:   R00 = 2   R01 = 4  R02 = 6
Flags: /
Subroutines: /
 
 
      STACK        INPUTS      OUTPUTS
           Y             /             n
           X             n        1  or  0
           L             /            d

 
Examples:

 n =     999983       XEQ "PR?"  >>>>   1      so  999983 is a prime
 n =  9999865009        R/S         >>>>   0      so  n  is not a prime   and the smallest non-trivial divisor = L = 10007
 
 

References:

[1]    F. Scheid -  "Numerical Analysis" - Mc Graw Hill   ISBN 07-055197-9
[2]    J-M  Ferrard - "Mastering your HP-48" - D3I Diffusion   ISBN 2-908791-04-8
[3]    R. Thangadurai - Mathematics Newsletter - Vol 17 #2 - "Irreducibility of Polynomials whose Coefficients are Integers"
[4]    Marian Trenkler - "An Algorithm for making Magic Cubes"
[5]    Marian Trenkler - "An Algorithm for Magic Tesseracts"
[6]    W. S. Andrews - "Magic Squares and Cubes"
[7]    Harm Derksen, Christian Eggermont, Arno van den Essen - "Multimagic Squares"
[8]    W. H. Benson & O. Jacoby - "Magic Cubes, New Recreations" - Dover - ISBN 0-486-24140-8
[9]    http://www.multimagie.com
[10]  Pierre Tougne - "Pour la Science" - pages 121 to 127 - Issue#46 - Aout 1981  ( in French )