hp41programs

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

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

     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<Y?
223  GTO 07
224  3
225  RCL 00
226  X^2
227  ST* Y
228   E6
229  /
230  +
231  1.001
232  +
233  REGMOVE
234  CLD
235  END

 
 ( 410 bytes / SIZE 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