The Museum of HP Calculators

# Matrix Operations for the HP-41

This program is Copyright © 2006 by Jean-Marc Baillard and is used here by permission.

This program is supplied without representation or warranty of any kind. Jean-Marc Baillard and The Museum of HP Calculators therefore assume no responsibility and shall have no liability, consequential or otherwise, of any kind arising from the use of this program material or any part thereof.

Overview

2°)  Multiplication
a)  Program#1
b)  Program#2
3°)  Matrix Inversion
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°) Storing a Matrix

-All these programs deal with real matrices only.
-Some of them use synthetic registers.
-Don't interrupt a routine which 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 into    R12   R15   R18   R21   respectively   and its control number =  11.02203
5   3   5   8                                    R13   R16   R19   R22

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

Notes:

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

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    is in registers  R02  R04  R06   and  B =   8  2  8   is in registers  R08  R10  R12  respectively

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

CF 02   1.00602  ENTER^
7.01202  ENTER^                                                         5  8  5                          R15  R17  R19
15       XEQ "M+-"  >>>>  15.02002  and   A+B =   9  7 17   is in registers   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      ( synthetic )
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

3  1
2  7  1  3                         4  2
Example:  Calculate  C = A.B where  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 -------------------- R15 thru R22         as the first register of C

R01  R04  R07  R10                2  7  1  3                                   R15  R19         3  1
In other words:      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" )

01  LBL "TM*"
02  STO O
03  STO P          ( synthetic )
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.

2  1  4                           3  1
7  9  6                           4  2
Example:  Calculate  C = tA.B where  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 -------------------- R15 thru R22         as the first register of C

R01  R05  R09                2  1  4                                   R15  R19         3  1
In other words:      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

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

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    ( a synthetic three-byte GTO )
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       must be stored into       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

Then you 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          will be in            R08  R13  R18  R23  R28         respectively     ( execution time = 1mn 36s if CF01
-5   19 -27  17 -4                                   R09  R14  R19  R24  R29                                  -------------- =  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 )

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

1  2  4                         R03  R06  R09                                                              1  4  1                        R12  R15  R18
Example:    A =   3  5  7   is in registers   R04  R07  R10   ( control number  3.01103 )   and  B = 5  9  2  is in registers  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^                                                                           15  11  -23            R21  R24  R27
21        XEQ "LIE"  >>>>  21.02903   and  [A,B] = AB - BA  =  24  19  -65    is in   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 ||

2  7  1  3             R01  R04  R07  R10
Example:     A =  1  9  4  2    is in   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)

1  2  4                         R03  R06  R09
Example:    A =   3  5  7   is in registers   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

2  7  1  3             R01  R04  R07  R10
Example:     A =  1  9  4  2    is in   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^                                                     2  1  4                         R21  R25  R29
21      XEQ "TRN"  >>>>   21.03204  and  tA =  7  9  6   is in registers   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

3  1  4                                     R01  R04  R07
Example:      A =   1  5  9    is stored in registers   R02  R05  R08
2  6  5                                     R03  R06  R09

R01  R04  R07                          3  1  2
1.00903  XEQ "TRN2"  >>>>   1.00903   and     R02  R05  R08   now contain     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

2  7  1  3             R01  R04  R07  R10
Example:     A =  1  9  4  2    is in   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:    A =  3  5  7    is in   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     is in     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.

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  /                                   If you don't have an HP-41CX, replace line 20 by
20  CLRGX                       ENTER^  SIGN  CLX  LBL 04  STO IND L  ISG L  GTO 04  X<>Y
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 )

1  4  9             R11  R14  R17
Example:    A =  3  5  7    is in   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   is in   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 )

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                                  If you have an X-Functions Module:
39  INT                                        -replace lines 38 to 40 by  RCL 06   REGMOVE   RCL 03
40  XEQ "MCO"                          -add  E3  /  RCL 04  +  RCL 06  E6  /  +  STO 06  after line 21
41   E3                                         -add  STO 06  after line 09
42  *
43  INT                                        The example below is then found in  8mn  instead of  9mn17s  ( choose bbb > 006 )
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?                                 The iterations stop when 2 consecutive partial sums are equal
75  GTO 01
76  RCL 02
77  END

( 123 bytes / SIZE 4n2+bbb )

 STACK INPUTS OUTPUTS X bbb.eeerr1 bbb.eeerr2

with  bbb > 005

1   2   3                          R11   R14   R17
Example:    If we store  A =   0   1   2    into registers   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
and  Exp(A) =  8.534640269    32.26024414   33.27906416    is in registers   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 )

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                                  If you have an X-Functions Module:
52  INT                                       -replace lines 51 to 53 by  RCL 06   REGMOVE   RCL 03
53  XEQ "MCO"                          -add  E3  /  RCL 04  +  RCL 06  E6  /  +  STO 06  after line 43
54  INT                                        -add  STO 06  after line 19
55  STO 05
56  RCL 02                                  The example below is then solved in 6mn03s  instead of  7mn06s  ( choose bbb > 006 )
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?                                     The iterations stop when 2 consecutive partial sums are equal
79  GTO 12
80  RCL 02
81  END

( 126 bytes / SIZE 4n2+bbb )

 STACK INPUTS OUTPUTS X bbb.eeerr1 bbb.eeerr2

with  bbb > 005

1.2   0.1   0.3                          R11   R14   R17
Example:    If we store  A =   0.1   0.8   0.1    into registers   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
and  Ln(A) =  0.097783005   -0.240971674   0.103424021    is in registers   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°) 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.

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              or      RCLFLAG    if you have an X-Functions Module
12  FIX 0
13  CF 28
14  CF 29
15  " A"
16  ARCL L
17  "~,"                  I mean   append  ,
18  ARCL Y
19  "~=?"               append  =  ?
20  STO d             or      STOFLAG    if you have an X-Functions Module
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

3   1   4   1
Example:    You want to store   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