Matrix Operations for the HP-41

Overview
 

 1°)  Addition & Subtraction
 2°)  Multiplication
   a)  Program#1
   b)  Program#2
 3°)  Matrix Inversion
   a)  Gauss-Jordan Elimination ( 2 programs )
   b)  Pseudo-Inverse: Greville's Method
 4°)  Lie Product
 5°)  Norm of a Matrix
 6°)  Trace of a Square Matrix
 7°)  Transpose
   a)  General case
   b)  Square Matrix
 8°)  Copying a Matrix
 9°)  Power of a Square Matrix
   a)  Program#1
   b)  Program#2 ( faster for large exponents )
10°) Exponential of a Square Matrix
11°) Logarithm of a Square Matrix
12°) Matrix Square-Root
13°) Matrix Pth Root
14°) Matrix Polynomial
15°) Storing a Matrix
 


>>> Latest Update: 2nd program in paragraph 3°a)

-All these programs deal with real matrices only.
-Some of them use synthetic registers.
-Don't interrupt a routine that uses registers P and/or Q  ( no risk of "crash" but their contents could be altered )

-The elements a_ij 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

-A routine is listed §15 to help you to store matrices in this way.

Notes:

-These routines do not check the control numbers that you enter in the stack.
-If you have an Advantage Module, most of the following programs will be probably unuseful for you...
 

1°) Addition & Subtraction
 

 CF 02  for an addition
 SF 02  for a subtraction
 
 

01  LBL "M+-"   02  STO M 03  STO N 04  CLX 05   E3 06  * 07  INT 08   E3

09  ST/ Y 10  LBL 01 11  CLX 12  RCL IND Z 13  RCL IND Y 14  FS? 02 15  CHS 16  +

17  STO IND M 18  CLX 19  SIGN 20  ST+ M 21  ST+ Z 22  ISG Y 23  GTO 01 24  RCL M

25  X<>Y 26  - 27  R^ 28   E3 29  * 30  FRC            31  + 32   E3

33  / 34  RCL N          35  + 36  CLA 37  END

 
   ( 62 bytes )
 
 

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

 CF 02   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

-Likewise for A-B after setting flag F02

-Unlike most of the following routines, the destination block of registers may be the same as the one of matrix A or B.
 

2°) Multiplication
 

     a) Program#1 (  product A.B )
 
 
 

01  LBL "M*"     02  STO O 03  STO P 04  RDN 05  STO N 06  X<>Y 07  STO M 08  FRC 09  ISG X 10  INT 11  STO Q

12  LBL 01 13  CLX 14  STO IND P 15  RCL M 16  RCL N 17  LBL 02 18  RCL IND Y 19  RCL IND Y 20  * 21  ST+ IND P 22  CLX

23  SIGN 24  + 25  ISG Y 26  GTO 02 27  LASTX        28  ST+ M 29  ST+ P 30  DSE Q 31  GTO 01 32  RCL M 33  FRC

34  ISG X 35  INT 36  STO Q 37  ST- M 38  ISG N 39  GTO 01        40  ENTER^ 41  SIGN 42  % 43  RCL P 44  LASTX

45  - 46  + 47   E3 48  / 49  RCL O        50  + 51  CLA 52  END

 
  ( 86 bytes )
 
 

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.
 

     b) Program#2  ( product ^tA.B )
 

-It is sometimes useful to multiply the transpose of a matrix A by another matrix B directly.
-This routine is used to solve underdetermined and overdetermined systems ( cf "Linear and non-Linear Systems for the HP-41" )
  and to perform greville's algorithm in §3-b).
 
 

01  LBL "TM*"   02  STO O 03  STO P 04  RDN 05  STO N 06  X<>Y 07  STO M 08  INT 09  LASTX 10  FRC 11  ENTER^ 12  ISG X 13  INT 14  X<>Y 15  PI 16  INT 17  10^X

18  * 19  INT 20  ISG X           21  CLX 22  R^ 23  - 24  RCL Y 25  ST- Y 26  / 27  PI 28  INT 29  10^X 30  / 31  STO Q 32  LBL 01 33  CLX 34  RCL M

35  RCL N 36  LBL 02 37  ENTER^ 38  CLX 39  STO IND P 40  RDN 41  LBL 03 42  RCL IND Y 43  SIGN 44  CLX 45  RCL IND Y 46  ST* L 47  X<> L 48  ST+ IND P 49  CLX 50  SIGN 51  ST+ Z

52  + 53  DSE Z 54  GTO 03 55  LASTX 56  ST+ P 57  RDN 58  ENTER^ 59  FRC 60  ISG X 61  INT 62  STO T          63  - 64  ISG Q 65  GTO 02 66  X<> Q 67  FRC 68  STO Q

69  RDN 70  ISG N 71  GTO 01 72  LASTX 73  INT 74  ENTER^ 75  SIGN 76  % 77  RCL P          78  DSE X 79  + 80   E3 81  / 82  RCL O 83  + 84  CLA 85  END

 
  ( 126 bytes )
 
 

STACK	INPUTS	OUTPUTS
Z	bbb.eeerr1	rr3
Y	bbb.eeerr2	rr3
X	bbb3	bbb.eeerr3

-We must have  rr₁  =  rr₂  , 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
 

3°) Matrix Inversion
 

     a) Gauss-Jordan Elimination
 

 "INV" can invert up to a 16x16 matrix, or even a 17x17 matrix if it's executed directly from extended memory.
  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 "INV"
 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.

 The XTOA  AROT  ATOX and REGSWAP functions of the X-Functions module are used.
 If you don't have an X-Functions module, you can delete lines 148 to 177 and lines 25 to 69 but your matrix will have to be "completely regular"
 I mean that every diagonal pivot must be non-zero. However, if the first element is 0 , you can swap row n°1 and row n°2  for instance ,
 and after the whole calculation , you will have to swap column n°1 and column n°2

 In fact, the HP-41 has to remember the different exchanges of rows in order to do the same exchanges of columns in the reverse order
 after the calculation is achieved, and  XTOA  is used for that purpose.

-Line 147 is a three-byte  GTO 14
 
 

01  LBL "INV"   02  STO 01   03  .1   04  %   05  ST+ 01   06  ST0 02   07  X<>Y   08  ST* Y   09   E-5   10  STO T   11  *   12  X<>Y   13  6.005   14  ST+ 02   15  +   16  STO 05   17  +   18  STO 03         19  +   20  STO 04   21  SIGN   22  STO 00   23  CLA   24  LBL 14   25  FS? 01   26  GTO 01   27  CLST   28  RCL 04   29  INT   30  RCL 02   31  FRC   32  +   33  RDN   34  +   35  LBL 00   36  CLX

37  RCL IND Z   38  ABS   39  X<=Y?   40  GTO 00   41  X<>Y   42  ENTER^   43  +   44  LBL 00   45  ISG Z   46  GTO 00   47  R^   48  INT   49  RCL 04   50  INT   51  -   52  RCL 03   53  STO Z   54  +   55  XTOA   56  X=Y?   57  GTO 01   58  LBL 09   59  RCL IND X   60  X<> IND Z   61  STO IND Y   62  ISG Y   63  RDN   64  ISG Y   65  GTO 09   66  RCL 00   67  CHS   68  STO 00   69  LBL 01   70  RCL 02   71  INT   72  RCL 05

73  INT   74  X=Y?   75  GTO 03   76  RCL 03   77  INT   78  X=Y?   79  GTO 02   80  RCL IND X   81  RCL IND T   82  *   83  RCL IND 04   84  /   85  ST- IND Z   86  LBL 02   87  ISG 02   88  GTO 04   89  RCL 01   90  INT   91  ST- 02   92  ST+ 03   93  GTO 04   94  LBL 03   95  RCL 01   96  INT   97  ST+ 03   98  ST+ 05   99  DSE 05 100  LBL 04 101  ISG 05 102  GTO 01 103  LBL 05 104  RCL 02 105  INT 106  RCL 04 107  INT 108  X=Y?

109  GTO 06 110  RCL IND X 111  ST/ IND Z 112  LBL 06 113  ISG 02 114  GTO 05 115  RCL 01 116  INT 117  X^2 118  ST- 03 119  LBL 07 120  RCL 03 121  INT 122  RCL 04 123  INT 124  X=Y? 125  GTO 08 126  RCL IND X 127  CHS 128  ST/ IND Z 129  LBL 08 130  ISG 03 131  GTO 07 132  RCL IND 04 133  ST* 00 134  1/X 135  STO IND 04 136  RCL 01 137  FRC 138  ST+ 02 139  LASTX 140  INT 141  X^2 142  ST- 03 143  ST- 05 144  SIGN

145  ST+ 03 146  ISG 04 147  GTO 14 148  FS? 01 149  GTO 11 150  RCL 01 151  INT 152  STO 04 153  STO 05 154  LBL 10 155  RCL 05 156  1 157  - 158  AROT 159  ATOX 160  6 161  - 162  RCL 04 163  ST* Z 164  * 165  6 166  ST+ Z 167  + 168  RCL 01 169  FRC 170  + 171   E3 172  / 173  + 174  REGSWAP 175  DSE 05 176  GTO 10 177  LBL 11 178  RCL 00 179  END

 
( 260 bytes / SIZE 6+n²)
 

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

-Put the order of the matrix into X  and execute "INV"

          5  XEQ "INV"     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     ( execution time = 1mn 36s if CF01
        -5   19 -27  17 -4                     R09  R14  R19  R24  R29                                  execution time =  1mn22s if SF01 )
         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 !
 ( Execution time = 42mn ) 

-In the following proogram, there is no pivoting.
-So, the elements are stored in registers R01 .... Rn^2

-Line 100 is a three-byte  GTO 01
 
 

01 LBL "INV"  02 STO M  03 .1  04 %  05 ST+ M  06 STO N  07 X<>Y  08 ST* Y  09  E-5  10 STO T  11 *  12 X<>Y  13 1  14 STO 00            15 ST+ N  16 +  17 STO Q  18 +  19 STO O  20 +  21 STO P

22 LBL 01  23 RCL N  24 INT  25 RCL Q  26 INT  27 X=Y?  28 GTO 03  29 RCL O  30 INT  31 X=Y?  32 GTO 02  33 RCL IND X  34 RCL IND T  35 *  36 RCL IND P  37 /  38 ST- IND Z  39 LBL 02  40 ISG N  41 GTO 04  42 RCL M

43 INT  44 ST- N  45 ST+ O  46 GTO 04  47 LBL 03  48 RCL M  49 INT  50 ST+ O  51 ST+ Q  52 DSE Q  53 LBL 04  54 ISG Q  55 GTO 01  56 LBL 05  57 RCL N            58 INT  59 RCL P  60 INT  61 X=Y?  62 GTO 06  63 RCL IND X

64 ST/ IND Z  65 LBL 06  66 ISG N  67 GTO 05  68 RCL M  69 INT  70 X^2  71 ST- O  72 LBL 07  73 RCL O  74 INT  75 RCL P            76 INT  77 X=Y?  78 GTO 08  79 RCL IND X  80 CHS  81 ST/ IND Z  82 LBL 08  83 ISG O

84 GTO 07  85 RCL IND P  86 ST* 00  87 1/X  88 STO IND P  89 RCL M            90 FRC  91 ST+ N  92 LASTX  93 INT  94 X^2  95 ST- O  96 ST- Q  97 SIGN  98 ST+ O  99 ISG P 100 GTO 01 101 RCL 00 102 CLA 103 END

 
    ( 171 bytes / SIZE 1+n²)
 

Example:   the same 5x5 Pascal's matrix M

          1  1  1   1  1                 R01  R06  R11  R16  R21
          1  2  3   4   5                R02  R07  R12  R17  R22
          1  3  6  10 15      =       R03  R08  R13  R18  R23        respectively
          1  4 10 20 35               R04  R09  R14  R19  R24
          1  5 15 35 70               R05  R10  R15  R20  R25

-Put the order of the matrix into X  and execute "INV"

          5  XEQ "INV"   >>>>   Det M = 1 =  R00                                         ---Execution time = 75s---

 and the inverse matrix:

          5  -10  10  -5   1                    R01  R06  R11  R16  R21
       -10  30 -35  19 -4                     R02  R07  R12  R17  R22
        10 -35  46 -27  6         =          R03  R08  R13  R18  R23         respectively
        -5   19 -27  17 -4                     R04  R09  R14  R19  R24                          
         1   -4    6   -4   1                     R05  R10  R15  R20  R25

Note:

-Synthetic registers P & Q are used.
-So do not stop this program ( no risk of crash but the result would be wrong )
 

     b) Pseudo-Inverse: Greville'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  a_k = the k-th column of A
        A_k = the matrix built with the k first columns of A
       A⁺_k  = the pseudoinverse at step k

-We start with  A⁺₁ = ^ta₁ / ( ^ta₁ a₁ )   if   a₁ # 0
                  or  A⁺₁ = ^ta₁                   if   a₁ = 0     and then,  for k = 2 , 3 , ..... , m

     Let    d_k = A⁺_k-1  a_k
              c_k = a_k - A_k-1 d_k
              b_k = ^tc_k / ( ^tc_k c_k )   if   c_k # 0      or      b_k = ^td_k A⁺_k-1 / ( 1 +  ^td_k d_k )   if  c_k = 0

   >>>   A⁺_k  = [ A⁺_k-1 - d_k b_k ]           In other words,  A⁺_k  is obtained by placing  the row   b_k  under the matrix  A⁺_k-1 - d_k b_k
                          [      b_k      ]

-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*"

-Line 208 is a three-byte  GTO 01
 
 

01  LBL "GRVL   02  STO 00   03  X<>Y   04  STO 01   05  INT   06  LASTX   07  FRC   08   E3   09  STO 08   10  *   11  INT   12  X<>Y   13  DSE X   14  -   15  RCL 01   16  FRC   17  ISG X   18  INT   19  STO 03   20  /   21  STO 02   22  RCL 00   23  +   24  DSE X   25  STO 05   26  RCL 01   27  RCL 02   28  RCL 03   29  ST+ T   30  ST* Y   31  -   32  RCL 08   33  /   34  -   35  STO 04

36  X<>Y   37  XEQ "TRN"   38  STO 06   39  RCL 04   40  RCL 05   41  XEQ "M*"   42  RCL IND X   43  SQRT   44   E-7    45  X>Y?   46  GTO 01   47  RCL 06   48  LASTX   49  LBL 00   50  ST/ IND Y   51  ISG Y   52  GTO 00   53  LBL 01   54  RCL 03   55  .1   56  %   57  +   58  RCL 04   59  STO 07   60  +   61  STO 04   62  RCL 06   63  DSE 02   64  X=0?   65  RTN    66  X<>Y   67  RCL 00   68  INT   69  XEQ "M*"   70  STO 00

71  RCL 07   72  FRC   73  RCL 01   74  INT   75  +   76  X<>Y   77  RCL 05   78  INT   79  XEQ "M*"   80  STO 05   81  RCL 08   82  *   83  INT   84  RCL 04   85  RCL 08   86  ST/ Z   87  *   88  INT   89  RCL 08   90  /   91  X<>Y   92  LBL 02   93  RCL IND Y   94  RCL IND Y   95  -   96  STO IND Y   97  ISG Y   98  RDN   99  ISG Y 100  GTO 02 101  RCL 05 102  RCL 05 103  7 104  XEQ "TM*" 105  GTO 04

106  LBL 03 107  RCL 00 108  RCL 06 109  RCL 05 110  INT 111  XEQ "TM*" 112  STO 05 113  RCL 00 114  RCL 00 115  7 116  XEQ "TM*" 117  1 118  ST+ 07 119  LBL 04 120  RCL 05 121  RCL 08 122  * 123  INT 124  RCL 08 125  / 126  STO 05 127  RCL 07 128  SQRT 129   E-7    130  X>Y? 131  GTO 03 132  RCL 05 133  LASTX 134  LBL 05 135  ST/ IND Y 136  ISG Y 137  GTO 05 138  RCL 00 139  RCL 05 140  RCL 06

141  FRC 142  RCL 08 143  * 144  INT 145  1 146  + 147  XEQ "M*" 148  RCL 08 149  * 150  INT 151  RCL 06 152  RCL 08 153  ST/ Z 154  * 155  INT 156  RCL 08 157  / 158  LBL 06 159  RCL IND Y 160  ST- IND Y 161  ISG Y 162  RDN 163  ISG Y 164  GTO 06 165  RCL 06 166  FRC 167  ISG X 168  INT 169  STO 07 170  RCL 05 171  FRC 172  RCL 08 173  * 174  INT 175  RCL 06

176  FRC 177  RCL 08 178  * 179  INT 180  RCL X 181  RCL 03 182  + 183  LBL 07 184  RCL IND Z 185  STO IND Y 186  CLX 187  RCL 07 188  STO 09 189  RDN 190  DSE Z 191  DSE X 192  LBL 08 193  RCL IND Y 194  STO IND Y 195  RDN 196  DSE Y 197  DSE X 198  DSE 09 199  GTO 08 200  X#Y? 201  GTO 07 202  RCL 03 203  RCL 08 204  / 205   E-5 206  + 207  ST+ 06 208  GTO 01 209  END

 
  ( 300 bytes )
 
 

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       --- execution time = 88s ---

                   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⁺ = Id₃

Notes:

-If A is a matrix with 4 rows and 6 columns, execution time ~  3mn40s  ( SIZE at least 083 )
-If A is a matrix with 6 rows and 4 columns, execution time ~  2mn10s  ( SIZE at least 079 )
-The coefficients of A are unchanged.

-If you want that the Moore-Penrose inverse is stored from register RBB ( R22 in our example ), replace line 65 by GTO 09
  and add   LBL 09  RCL 00  INT  XEQ "MCO"   after line 208
-Lines 44 and 129 contain small numbers to identify the columns  c_k that verify  ^tc_k c_k ~ 0
-In some cases,  E-7  may be too small - or too large - for instance if the program produces huge coefficients for A⁺ when the coefficients of A are of order of 1
-Change these lines by another "small" number or replace them by RCL 10 and store your own small-number into R10 before executing "GRVL"
-Of course, in this case, choose bbb > 010
 

4°) Lie Product
 

-This program calculates  [A,B] = AB - BA  where A & B are 2 square-matrices.
-It uses CX functions.

-"M*" is called as a subroutine.
 
 

01  LBL "LIE" 02  RDN 03  STO 02 04  X<>Y 05  STO 01 06  R^ 07  XEQ "M*" 08  STO 00

09  FRC 10  ISG X 11  INT 12  X^2 13  "T" 14  CRFLD 15  RCL 00 16  SAVERX

17  RCL 01 18  RCL 02 19  RCL 00 20  INT 21  XEQ "M*" 22   E3 23  * 24  INT

25   E3 26  / 27  SIGN 28  CLX 29  SEEKPT 30  LBL 01 31  GETX 32  ST- IND L

33  ISG L 34  GTO 01 35  "T" 36  PURFL 37  RCL 00      38  CLA 39  END

 
  ( 66 bytes )
 
 

STACK	INPUTS	OUTPUTS
Z	bbb.eeerr1	rr3
Y	bbb.eeerr2	rr3
X	bbb3	bbb.eeerr3

-Every  bbb  > 002

Example:

            1  2  4         R03  R06  R09                                                               1  4  1          R12  R15  R18
   A =   3  5  7    =   R04  R07  R10   ( control number  3.01103 )   and  B = 5  9  2    =    R13  R16  R19      ( c.n. = 12.02003 )
            7  9  8         R05  R08  R11                                                               6  5  3          R14  R17  R20

-If you want to store [A,B] in registers R21 ...etc...

  3.01103   ENTER^
 12.02003  ENTER^
      21        XEQ "LIE"  >>>>  21.02903

                                         15  11  -23         R21  R24  R27
and  [A,B] = AB - BA  =  24  19  -65    =   R22  R25  R28      ( execution time = 28 seconds )
                                         58  85  -34         R23  R26  R29
 

5°) Norm of a Matrix
 

-This short routine computes  || A || = ( SUM_i,j a²_i,j )^1/2
 
 

01  LBL "MABS"  02   E3  03  *  04  INT  05   E3  06  /  07  0  08  LBL 01  09  RCL IND Y  10  X^2  11  +  12  ISG Y  13  GTO 01  14  SQRT  15  END

 
  ( 29 bytes )
 
 

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 "MABS"  >>>>  || A || = 14.8996643
 

6°) Trace of a Square Matrix
 

-The trace of a square matrix equals the sum of its diagonal elements.
 
 

01  LBL "TRACE"  02   E-5  03  +  04  0  05  LBL 01  06  RCL IND Y  07  +  08  ISG Y  09  GTO 01  10  END

 
  ( 25 bytes )
 
 

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
 

7°) Transpose
 

     a) General Case
 

-The transpose of a nxm matrix A = [ a_{i j} ] is the mxn matrix ^tA = [ b_{i j} ] defined by  b_{i j} = a_{j i}
-"TRN" stores the transpose of a matrix  in a different block of registers.
-The 2 blocks cannot overlap.
 
 

01  LBL "TRN"   02  STO O 03  RCL Y 04  STO T 05  FRC 06  ISG X 07  INT 08  STO M

09  STO N 10  RDN 11  LBL 01 12  RCL IND Y 13  STO IND Y 14  CLX 15  SIGN 16  +

17  ISG Y 18  GTO 01       19  RDN 20  SIGN 21  + 22  ENTER^ 23  R^ 24  DSE N

25  GTO 01 26  STO Y 27  RCL O 28  - 29  RCL M        30  / 31  DSE Y 32   E2

33  / 34  + 35   E3 36  / 37  RCL O        38  + 39  CLA 40  END

 
   ( 67 bytes )
 
 

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
 

     b) Square Matrix
 

-Unlike "TRN" , the following program replaces a square matrix by its transpose.
 
 

01  LBL "TRN2"  02  STO L  03  LBL 01  04  ENTER^  05  ENTER^  06  LBL 02  07  RCL IND X  08  X<> IND Z  09  STO IND Y  10  CLX  11  1  12  ST+ Y  13  RDN  14  ISG Y  15  GTO 02  16  R^  17  ST+ T  18  R^  19  ISG X  20  GTO 01  21  LASTX  22  END

 
   ( 41 bytes )
 
 

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
 

8°) Copying a Matrix
 
 
 

01  LBL "MCO" 02  RCL Y 03   E3 04  * 05  INT 06   E3 07  /

08  SIGN 09  RDN 10  ENTER^ 11  ENTER^ 12  LBL 01 13  CLX 14  RCL IND L

15  STO IND Y 16  ISG Y 17  CLX 18  ISG L 19  GTO 01 20  CLX 21  SIGN

22  - 23   E3 24  / 25  + 26  X<>Y          27  FRC 28  ISG X

29  INT           30   E5 31  / 32  + 33  END

 
( 52 bytes )
 
 

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

-The 2 blocks of registers cannot overlap unless  bbb₂ <  bbb₁
 

9°) Power of a Square Matrix
 

     a) Program#1
 

-The program hereafter calculates A^p  where A is a square matrix and p is a positive integer   ( A⁰ = 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
 
 

01  LBL "MPOW" 02  STO 00 03  X<>Y 04  STO 01 05  ENTER^ 06  ENTER^ 07  FRC

08  ISG X 09  INT 10  X^2 11  STO 03 12  + 13  INT 14  XEQ "MCO"

15  STO 02 16  INT 17  ST+ 03 18  LBL 01 19  RCL 01           20  RCL 02 21  DSE 00

22  X=0? 23  RTN 24  RCL 03 25  XEQ "M*"       26  DSE 00 27  X=0? 28  RTN

29  RCL 01 30  RCL 02 31  INT 32  XEQ "M*"       33  GTO 01 34  END

 
  ( 58 bytes / SIZE 3n²+bbb )
 
 

STACK	INPUTS	OUTPUTS
Y	bbb.eeerr1	/
X	p > 0	bbb.eeerr2

    ( bbb > 003 )

                            1  4  9            R11  R14  R17
Example:

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

  11.01903  ENTER^
       7          XEQ "MPOW"  >>>>  20.02803 = the control number of A⁷    ( in 1mn20s )

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

-This program requires a "reasonable" time if p is small
-But if A is a 3x3 matrix and p = 100 , execution time = 21mn29s
-The following routine is faster in this case.
 

     b) Program#2
 

-"MPOW2" minimizes the number of multiplications.
-It also uses "M*" and "MCO" as subroutines.

-If you don't have an HP-41CX, replace line 20 by    ENTER^  SIGN  CLX  LBL 04  STO IND L  ISG L  GTO 04  X<>Y
 
 

01  LBL "MPOW2" 02  STO 00 03  X<>Y 04  STO 01 05  ENTER^ 06  FRC 07  ISG X 08  INT 09  X^2 10  STO 03 11  1.001 12  * 13  +

14  STO 02 15   E3 16  * 17  INT 18   E3 19  /     20  CLRGX            21  INT 22  ST+ 03 23  RCL 02 24   E-5 25  + 26  1

27  LBL 00 28  STO IND Y      29  ISG Y 30  GTO 00 31  LBL 01 32  RCL 00 33  2 34  / 35  STO 00 36  FRC 37  ST- 00 38  X=0? 39  GTO 02

40  RCL 01 41  RCL 02 42  RCL 03 43  XEQ "M*" 44  RCL 02 45  INT 46  XEQ "MCO"    47  LBL 02 48  RCL 00 49  X=0? 50  GTO 03 51  RCL 01 52  RCL 01

53  RCL 03 54  XEQ "M*" 55  RCL 01 56  INT 57  XEQ "MCO"    58  GTO 01 59  LBL 03 60  RCL 02 61  END

 
   ( 103 bytes / SIZE 3n²+bbb )
 
 

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⁷
            2  1  8          R13  R16  R19

  11.01903  ENTER^
       7          XEQ "MPOW2"  >>>>  20.02803 = the control number of A⁷    ( in 1mn38s )

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

-This program also works if p = 0.
-"MPOW2" is actually slower than "MPOW" for n = 3 ( the order of the matrix ) and p = 7
-But if n = 3 and p = 100 , execution time is only 2mn31s ( instead of 21mn29s with "MPOW"! )
 

10°) Exponential of a Square Matrix
 

-The exponential of a nxn matrix A is defined as

   exp(A) = I + A + A²/2! + A³/3! + ..... + A^k/k! + ....                ( I = the Unit matrix = the Identity matrix:  1 in the diagonal, 0 elsewhere )

- If you have an X-Functions Module:

   replace lines 38 to 40 by  RCL 06   REGMOVE   RCL 03
   add  E3  /  RCL 04  +  RCL 06  E6  /  +  STO 06  after line 21
   add  STO 06  after line 09

-The example below is then solved in  8mn  instead of  9mn17s  ( choose bbb > 006 )
 
 

01  LBL "EXPM" 02  STO 01 03  ENTER^ 04  ENTER^ 05  FRC 06  ISG X 07  INT 08  X^2 09  STO 04 10  + 11  INT 12  XEQ "MCO" 13  ENTER^ 14  STO 02 15  RCL 04 16  +

17  INT 18  XEQ "MCO" 19  STO 03 20  INT 21  ST+ 04 22  RCL 02 23   E-5 24  + 25  1 26  STO 00 27  LBL 00 28  ST+ IND Y 29  ISG Y 30  GTO 00 31  LBL 01 32  ISG 00

33  CLX 34  RCL 01 35  RCL 03 36  RCL 04 37  XEQ "M*" 38  RCL 03   39  INT    40  XEQ "MCO" 41   E3   42  * 43  INT   44   E3 45  / 46  RCL 00 47  LBL 02 48  ST/ IND Y

49  ISG Y 50  GTO 02 51  RCL 02 52   E3 53  * 54  INT 55   E3 56  / 57  RCL 03 58  INT 59  STO 05 60  CLX 61  LBL 03 62  RCL IND 05 63  RCL IND Z 64  +

65  STO IND Z  66  LASTX 67  - 68  ABS 69  + 70  ISG 05 71  CLX 72  ISG Y 73  GTO 03 74  X#0?    75  GTO 01 76  RCL 02 77  END

 
  ( 123 bytes / SIZE 4n²+bbb )
 
 

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 = [ a_ij ] is a diagonal matrix too where the diagonal elements are exp(a_ii)
 

11°) Logarithm of a Square 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 + ( A - I )³/3 -  ( A - I )⁴/4 + ......                   ( I = the Unit matrix = the Identity matrix:  1 in the diagonal, 0 elsewhere )

-If you have an X-Functions Module:

  replace lines 51 to 53 by  RCL 06   REGMOVE   RCL 03
  add  E3  /  RCL 04  +  RCL 06  E6  /  +  STO 06  after line 43
  add  STO 06  after line 19

-The example below is then solved in 6mn03s  instead of  7mn06s  ( choose bbb > 006 )
 
 

01  LBL "LNM" 02  STO 01 03   E-5 04  + 05  1 06  CHS 07  STO 00 08  LBL 00 09  ST+ IND Y   10  ISG Y 11  GTO 00 12  RCL 01 13  ENTER^ 14  ENTER^ 15  FRC 16  ISG X 17  INT

18  X^2 19  STO 04 20  + 21  INT 22  XEQ "MCO" 23  STO 02 24  RCL 01 25   E3 26  * 27  INT 28   E3 29  / 30  RCL 00 31  LBL 01 32  ST* IND Y 33  ISG Y 34  GTO 01

35  RCL 01 36  RCL 02 37  RCL 04 38  + 39  INT 40  XEQ "MCO" 41  STO 03 42  INT 43  ST+ 04 44  LBL 12 45  DSE 00 46  CLX 47  RCL 01 48  RCL 03 49  RCL 04 50  XEQ "M*" 51  RCL 03

52  INT    53  XEQ "MCO"  54  INT    55  STO 05 56  RCL 02   57   E3 58  * 59  INT 60   E3 61  / 62  0 63  LBL 02 64  RCL IND 05 65  RCL 00 66  / 67  RCL IND Z 68  +

69  STO IND Z   70  LASTX 71  - 72  ABS 73  + 74  ISG 05 75  CLX 76  ISG Y 77  GTO 02 78  X#0?   79  GTO 12 80  RCL 02 81  END

 
  ( 126 bytes / SIZE 4n²+bbb )
 
 

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)     ( in 7mn06s )

                   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
 

12°) Matrix Square-Root
 

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

    M₀ = A      M_k+1 = (1/2) [ I + ( M_k + M_k^-1 ) / 2 ]      where  I  is the Identity matrix
    X₀ = A      X_k+1 = X_k ( I + M_k^-1 ) / 2

-If A has no negative real eigenvalue,  X_k quadratically tends to  A^1/2  and  M_k 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, A^1/2  has replaced A in registers R01 to Rn²
 

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

                                  •  R01   thru   •  Rn² =  the coefficients of the matrix A, in column order.                Rn²+1 to R4n²  temp

Flag:  F07
Subroutines:  "MCO"  "M*"  "LS2"  ( cf "Linear and non-linear systems for the HP-41" )

-Lines 20-21-30-51-52-76 are not essential: they only save time and they could be deleted
-If you don't have an HP-41CX, replace line 39 by  ENTER^   ENTER^   CLX   LBL 09   STO IND Z   ISG Z   GTO 09   CLX
-Line 75 is a three-byte GTO but it doesn't really matter
-Line 177 is a three-byte  GTO 01
-Lines 186 to 190 can be deleted: the output will be  1.eee01  instead of  1.eee,nn
 
 

01  LBL "SQRTM"   02  RCL 00   03  X^2   04   E3   05  /   06  RCL 00   07  X^2   08  ST+ X   09  1   10  ST+ Z   11  +   12  XEQ "MCO"   13  RCL 00   14  X^2   15  3   16  *   17  1   18  +   19  XEQ "MCO"   20  SF 07     21  GTO 00   22  LBL 01   23  3.004   24  RCL 00   25  X^2   26  *   27  1   28  ST+ Y   29  XEQ "MCO"   30  LBL 00   31  RCL 00   32  X^2

33  ENTER^   34  ST+ X   35   E3   36  /   37  +   38  ISG X   39  CLRGX    40  RCL 00   41  1   42  +   43   E5   44  /   45  +   46  SIGN   47  LBL 02   48  STO IND L       49  ISG L   50  GTO 02   51  FS?C 07   52  GTO 00   53  RCL 00   54  X^2   55  ST+ X   56  X<> 00   57  X^2   58  0   59  LBL 03   60  RCL IND Y   61  RCL IND 00   62  -   63  ABS   64  +

65  DSE 00   66  DSE Y   67  GTO 03   68  RCL 00   69  SQRT   70  STO 00   71  CLX   72   E-9   73  VIEW Y   74  X>Y?   75  GTO 08     76  LBL 00   77  RCL 00   78  ENTER^   79  ST+ X   80  XEQ "LS2"        81  4   82  R^   83  SQRT   84  INT   85  STO 00   86  X^2   87  ST* Y   88  ENTER^   89  ST+ Y   90   E3   91  /   92  +   93  LBL 04   94  RCL IND Y   95  RCL IND Y   96  +

97  4   98  /   99  STO IND Z 100  DSE Z 101  RDN 102  DSE X 103  GTO 04 104  CLX 105   E3 106  / 107  RCL 00 108  X^2 109  ST+ X 110  + 111  1 112  ST+ Y 113  XEQ "MCO"   114  RCL 00 115  DSE X 116   E5 117  / 118  + 119  RCL 00 120  X^2 121  ENTER^ 122  ST+ Y 123   E3 124  / 125  + 126  RCL 00 127  1 128  +

129   E5 130  / 131  + 132  SIGN 133  LBL 05 134  ST+ IND L     135  DSE L 136  GTO 05 137  CLX 138  RCL 00 139  X^2 140  .1 141  % 142  + 143  + 144  RCL 00 145  X^2 146  ST+ X 147  1 148  + 149  XEQ "M*" 150   E3 151  * 152  INT 153   E3 154  / 155  2 156  LBL 06 157  ST/ IND Y 158  ISG Y 159  GTO 06 160  3.004

161  RCL 00 162  X^2 163  * 164  RCL 00 165  1 166  ST+ Z 167  + 168   E5 169  / 170  + 171  X<>Y 172  1/X 173  LBL 07 174  ST+ IND Y     175  ISG Y 176  GTO 07 177  GTO 01 178  LBL 08 179  2.003 180  RCL 00 181  X^2 182  * 183  1 184  ST+ Y 185  XEQ "MCO" 186  RCL 00 187  DSE X 188   E5 189  / 190  + 191  CLD 192  END

 
   ( 307 bytes / SIZE 4n²+1 )
 
 

STACK	INPUTS	OUTPUTS
X	/	1.eee,nn

   with  eee = n²
 
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"   the HP-41 displays the successive || M_k - I ||          ( here, the norm is simply the sum of the magnitude of the elements )

                             5.0000
                             0.7716
                             0.0380
                             0.0001
                             6.9269 E-10

-Finally, the control number  1.00903  is returned in X-register.      --- Execution time = 4mn38s ---

                        R01  R04  R07              1.794981016   0.656367530   0.540425260
   sqty(A)   =    R02  R05  R08      =      0.772143191    1.061374174   1.582989341
                        R03  R06  R09              0.501888909    0.231634627   1.833600669
 

Notes:

-This routine doesn't work for n = 1. Simply replace lines 115 & 187 by  1  -   ( but SQRT is much faster... )
-A few bytes may be saved if you prefer "LS3" to "LS2"
-Use the following listing:

-If you don't have an HP-41CX, replace line 27 by  ENTER^   ENTER^   CLX   LBL 09   STO IND Z   ISG Z   GTO 09   CLX
-Lines 65 & 168 are three-byte GTOs
-Lines 177 to 181 can be deleted: the output will be  1.eee01  instead of  1.eee,nn
 
 

01  LBL "SQRTM"   02  RCL 00   03  X^2   04   E3   05  /   06  RCL 00   07  X^2   08  ST+ X   09  1   10  ST+ Z   11  +   12  XEQ "MCO"     13  RCL 00   14  X^2   15  3   16  *   17  1   18  +   19  XEQ "MCO"   20  SF 07   21  LBL 01   22  RCL 00   23  X^2   24   E3   25  /   26  ISG X   27  CLRGX      28  RCL 00   29  1   30  +   31   E5   32  /   33  +   34  SIGN   35  LBL 02   36  STO IND L   37  ISG L

38  GTO 02   39  FS?C 07   40  GTO 00   41  RCL 00   42  ST+ X   43  X^2   44  X<> 00   45  X^2   46  0   47  LBL 03   48  RCL IND Y   49  RCL IND 00     50  -   51  ABS   52  +   53  DSE 00   54  DSE Y   55  GTO 03   56  RCL 00   57  3   58  /   59  SQRT   60  STO 00   61  CLX   62   E-9   63  VIEW Y   64  X>Y?   65  GTO 08   66  LBL 00   67  3.004   69  RCL 00   69  X^2   70  ST*Y   71  1   72  ST+ Z   73  +   74  XEQ "MCO"

75  RCL 00   76  ENTER^   77  ST+ X   78  XEQ "LS3"   79  R^   80  INT   81  STO 00   82  ST+ X   83  X^2   84  RCL 00   85  X^2   86  LBL 04   87  RCL IND Y      88  RCL IND Y   89  +   90  4   91  /   92  STO IND Z   93  DSE Z   94  RDN   95  DSE X   96  GTO 04   97  CLX   98   E3   99  / 100  RCL 00 101  X^2 102  ST+ X 103  + 104  RCL 00 105  X^2 106  1 107  ST+ Z 108  + 109  XEQ "MCO" 110  RCL 00 111  DSE X

112   E5 113  / 114  + 115  RCL 00 116  X^2 117  LASTX 118  1 119  + 120   E5 121  / 122  + 123  SIGN 124  LBL 05 125  ST+ IND L     126  DSE L 127  GTO 05 128  CLX 129  RCL 00 130  X^2 131  .1 132  % 133  + 134  - 135  RCL 00 136  X^2 137  ST+ X 138  1 139  + 140  XEQ "M*" 141   E3 142  * 143  INT 144   E3 145  / 146  2 147  LBL 06 148  ST/ IND Y

149  ISG Y 150  GTO 06 151  3.004 152  RCL 00 153  X^2 154  * 155  RCL 00 156  1 157  ST+ Z 158  + 159   E5 160  / 161  + 162  X<>Y 163  1/X 164  LBL 07 165  ST+ IND Y     166  ISG Y 167  GTO 07 168  GTO 01 169  LBL 08 170  2.003 171  RCL 00 172  X^2 173  * 174  1 175  ST+ Y 176  XEQ "MCO" 177  RCL 00  178  DSE X 179   E5 180  / 181  + 182  CLD 183  END

 
  ( 293 bytes / SIZE 4n²+1 )
 
 

STACK	INPUTS	OUTPUTS
X	/	1.eee,nn

     ( eee = n² )

-With the same example, it yields:

                      1.794981016   0.656367529   0.540425261
    sqrt(A)  =   0.772143191   1.061374174   1.582989342                --- Execution time = 4mn33s ---
                      0.501888909   0.231634627   1.833600670

-There are minor differences in the last places for some coefficients.
 

13°) Matrix Pth Root  ( X-Functions Module required )
 

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

    M₀ = A      M_k+1 = M_k [ ( 2.I + ( p - 2 ) M_k ) ( I + ( p - 1) M_k ) ^-1 ] ^p      where  I  is the Identity matrix
    X₀ =  I       X_k+1  = X_k  ( 2.I + ( p - 2 ) M_k ) ^-1 ( I + ( p - 1) M_k )

    M_k  tends to  I        as k tends to infinity
    X_k   tends to A ^1/p  as k tends to infinity

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

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

                                  •  R01   thru   •  Rn² =  the coefficients of the matrix A, in column order.                Rn²+1 to R5n²+1 temp  ( R5n²+1 = p )

          >>>>   When the program stops,  A^1/p  is in registers  R01 thru Rn²

Flag:  F10
Subroutines:    "M*"   ( cf §2a )  &   "LS3"  ( cf "Linear and non-linear systems for the HP-41" )

-Line 223 is a three-byte  GTO 07
 
 

01  LBL "MROOT"   02  RCL 00   03  X^2   04  5   05  *   06  1   07  ST- Z   08  +   09  X<>Y   10  STO IND Y       11  .004001   12  RCL 00   13  X^2   14  *   15  1.001   16  +   17  REGMOVE   18  3.004   19  RCL 00   20  X^2   21  *   22  ISG X   23  XEQ 05   24  GTO 07   25  LBL 01   26   E3   27  /   28  +   29  RCL IND Y   30  FC? 10   31  DSE X   32  LBL 02   33  ST* IND Y   34  DSE Y   35  GTO 02   36  CLX   37  RCL 00   38  1   39  +   40   E5   41  /   42  +   43  RCL 00   44  X^2   45  +   46  SIGN   47  FC? 10

48  ST+ X   49  LBL 03   50  ST+ IND L   51  DSE L   52  GTO 03   53  FC?C 10   54  RTN   55  LBL 04   56  RCL 00   57  X^2   58   E3   59  /   60  ISG X   61  XEQ 05   62  RCL 00   63  ENTER^   64  ST+ X   65  XEQ "LS3"        66  R^   67  INT   68  STO 00   69  X^2   70  RTN   71  LBL 05   72  CLRGX   73  RCL 00   74  1   75  +   76   E5   77  /   78  +   79  SIGN   80  LBL 06   81  STO IND L   82  ISG L   83  GTO 06   84  RTN   85  LBL 07   86  4.001001   87  RCL 00   88  X^2   89  *   90  1.001   91  +   92  REGMOVE   93  RCL 00   94  X^2

95  ENTER^   96  ST+ Y   97  ST+ Z   98  SF 10   99  XEQ 01 100  4.002001 101  RCL 00 102  X^2 103  * 104  1.001 105  + 106  REGMOVE     107  3 108  RCL 00 109  X^2 110  ST+ Z 111  ST* Y 112  ST+ X 113  XEQ 01 114  LASTX 115  RCL 00 116  X^2 117   E3 118  / 119  + 120  ISG X 121  LASTX 122  ISG X 123  RCL 00 124   E-5 125  ST- T 126  * 127  + 128  RCL 00 129  X^2 130  1 131  + 132  XEQ "M*" 133  STO 04 134  SIGN 135  RCL 00 136  X^2 137  5 138  * 139  STO 03 140  + 141  RCL IND X

142  STO 02 143  2.004001 144  RCL 00 145  X^2 146  * 147  1.001 148  + 149  X<> 03 150  .1 151  % 152  + 153  RCL 00 154  X^2 155  - 156  ISG X 157  RCL 00            158   E5 159  / 160  + 161  STO 01 162  ISG 02 163  LBL 08 164  RCL 01 165  RCL 04 166  RCL 00 167  X^2 168  ST+ X 169  1 170  + 171  XEQ "M*" 172  RCL 03 173  REGMOVE 174  DSE 02 175  GTO 08 176  XEQ 04 177  4.004 178  X<>Y 179  ST* Y 180  ST- Y 181   E3 182  / 183  RCL 00 184   E5 185  / 186  1 187  + 188  ST+ Z

189  + 190  RCL 00 191  X^2 192  ST+ X 193  1 194  + 195  XEQ "M*" 196  2.003001 197  RCL 00 198  X^2 199  STO 01 200  * 201  1.001 202  + 203  REGMOVE     204  5 205  ST* 01 206  2.001 207  RCL 00 208  X^2 209  * 210  0 211  LBL 09 212  RCL IND Y 213  RCL IND 01 214  - 215  ABS 216  + 217  DSE 01 218  DSE Y 219  GTO 09 220   E-8 221  VIEW Y 222  X<Y? 223  GTO 07 224  3 225  RCL 00 226  X^2 227  ST* Y 228   E6 229  / 230  + 231  1.001 232  + 233  REGMOVE 234  CLD 235  END

 
 ( 410 bytes / SIZE 5n²+2 )
 
 

STACK	INPUTS	OUTPUTS
X	p	/

 
Example:      Calculate  A^1/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 "SQRTM"   the HP-41 displays the successive || M_k - I ||  and eventually,
 

                           R01  R04  R07              1.437771414   0.396760708   0.231139046
    A^(1/3)   =     R02  R05  R08     =       0.421053139   0.977706016   0.944423239            --- Execution time = 8mn04s ---
                           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  X_k+1 = (1/p) [ ( p - 1 ) X_k + A/X_k^p-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           --- Execution time = 29mn22s ---
                    1  4  2  9                           -0.224295471     0.139205813   0.171840655    1.642919190 
 

14°) Matrix Polynomial
 

-Given a polynomial  p(x) = a₀ x^m + a₁ x^m-1 + ............... + a_m     and a  nxn matrix A,
 "PMAT" calculates  p(A) = a₀ A^m + a₁ A^m-1 + ............... + a_m 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 = a₀  .......   •  Ree = a_m
                              •  Rbb'   thru   •  Ree' =  the coefficients of the matrix A, in column order.

                                  the 2n²  registers from  RBB: temp

Flag:  F07
Subroutines:  "MCO"  "M*"

-If you don't have an HP-41CX, replace line 19 by  ENTER^   ENTER^   CLX   LBL 00   STO IND Z   ISG Z   GTO 00   CLX
 
 

01  LBL "PMAT" 02  STO T 03  RDN 04  STO 01 05  FRC 06  ISG X 07  INT 08  STO 02 09  X^2

10  X<>Y 11  STO 03         12  RDN 13  + 14  STO 00 15  DSE X 16   E3 17  / 18  +

19  CLRGX  20  RCL 02         21   E5 22  / 23  + 24  STO 02 25  LBL 01 26   E-5 27  +

28  RCL IND 03 29  LBL 02 30  ST+ IND Y 31  ISG Y 32  GTO 02 33  RCL 01 34  RCL 02 35  ISG 03 36  X=0?

37  RTN 38  RCL 00 39  XEQ "M*" 40  RCL 02 41  INT 42  XEQ "MCO" 43  GTO 01 44  END

 
  ( 74 bytes )
 
 

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 x⁴ - x³ + 3 x² - 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             --- execution time = 64 seconds ---

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

15°) Storing a Matrix
 

-You can use for instance the "IO" routine listed in "Linear and non-Linear Systems for the HP-41" ( paragraph 3 )
  but the following program is more specific to matrices.

-Synthetic RCL d  & STO d  may be replaced by  RCLFLAG  &  STO FLAG   if you have an X-Functions Module
-The append character is denoted  ~
 
 

01  LBL "MSTO" 02  SF 10 03  1.001 04  ST* Y 05  LBL 01 06  ENTER^ 07  FRC 08   E3 09  ST* Y 10  CLX 11  RCL d  12  FIX 0 13  CF 28 14  CF 29 15  " A" 16  ARCL L

17  "~,"    18  ARCL Y 19  "~=?"   20  STO d  21  R^ 22  R^ 23  PROMPT 24  FC?C 22 25  GTO 02 26  STO IND Z  27  CLX 28   E-5 29  FC? 10 30  CLX 31  ST+ Z 32  SIGN

33  ST+ Z 34  + 35  FS? 10 36  GTO 01 37  RCL Y 38  FRC 39  ISG X 40  INT 41  RCL Y 42  INT 43  - 44  X<0? 45  GTO 03        46  RDN 47  GTO 01 48  LBL 02

49  FC?C 10 50  GTO 04 51  ENTER^ 52  LBL 03 53  CLX 54   E-3 55  + 56  FRC 57  ISG X 58  GTO 01        59  LBL 04 60  SIGN 61  - 62  INT 63  LASTX 64  FRC

65   E3 66  * 67  FRC 68  ST+ Y 69  X<> L 70  INT 71  X<>Y          72   E3 73  / 74  + 75  CLA 76  END

 
( 132 bytes )
 
 

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