MatrixOp

# 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 aij 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  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

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+n2)

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+n2)

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

-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+ = Id3

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  ck that verify  tck ck ~ 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 || = ( SUMi,j a2i,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 = [ 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.

 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  bbb2 <  bbb1

9°) Power of a Square Matrix

a) Program#1

-The program hereafter 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

 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 3n2+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  A7
2  1  8          R13  R16  R19

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

7851276   8652584   31076204            R20  R23  R26
A7 =   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 3n2+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  A7
2  1  8          R13  R16  R19

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

7851276   8652584   31076204          R20  R23  R26
A7 =   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 + A2/2! + A3/3! + ..... + Ak/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 4n2+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 = [ aij ] is a diagonal matrix too where the diagonal elements are exp(aii)

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/2 + ( A - I )3/3 -  ( A - I )4/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 4n2+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:

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*"  "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 4n2+1 )

 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"   the HP-41 displays the successive || Mk - 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 4n2+1 )

 STACK INPUTS OUTPUTS X / 1.eee,nn

( eee = n2 )

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

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

( 410 bytes / SIZE 5n2+2 )

 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 "SQRTM"   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            --- 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  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           --- Execution time = 29mn22s ---
1  4  2  9                           -0.224295471     0.139205813   0.171840655    1.642919190

14°) Matrix Polynomial

-Given a polynomial  p(x) = a0 xm + a1 xm-1 + ............... + am     and a  nxn matrix A,
"PMAT" calculates  p(A) = a0 Am + a1 Am-1 + ............... + am 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 = a0  .......   •  Ree = am
•  Rbb'   thru   •  Ree' =  the coefficients of the matrix A, in column order.

the 2n2  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 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             --- 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