The Museum of HP Calculators
The Museum of HP Calculators
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
1°) Addition & Subtraction
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...
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 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
Go back to the HP-41 software library
Go back to the general software library
Go
back to the main exhibit hall