Overview
-The elements aij of a matrix 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
-"MSTO" may help you to store matrices in this way.
-Likewise, the control number of a polynomial is bbb.eee
-If you have used "MSTO" to store its coefficients, subtract
rr E-5.
| XROM | Function | Desciption |
| 08,00
08,01 08,02 08,03 08,04 08,05 08,06 08,07 08,08 08,09 08,10 08,11 08,12 08,13 08,14 08,15 08,16 08,17 08,18 08,19 08,20 08,21 08,22 08,23 08,24 08,25 08,26 08,27 08,28 08,29 08,30 08,31 08,32 08,33 08,34 08,35 08,36 08,37 08,38 08,39 08,40 08,41 08,42 08,43 08,44 08,45 08,46 08,47 08;48 08,49 08,50 08,51 08,52 08,53 08,54 08,55 08,56 |
-JMB MATRX
"CRMAT" "MCO" "MSTO" "RANM" "RANSYM" "TRN" "TRN2" -MTR FNS "EXPM" "GRVL" "LNM" "M-" "M+" "M*" "MINV" "MNORM" "MPOW" "MROOT" "PMAT" "SQRTM" "TM*" "TRACE" -EIGENVAL "DFL" "JCB" "JCBZ" "PCAR" "PMIN" -LNR SYS "D0" "D3" "D4" "D5" "DET" "LS" "LS3" "LS3-" "LS3+" -MAGIC MTR "MG4DHC" "MGCB" "MGSQ" "PMGCB" "PMGSQ" -POLYNOM "BRST" "COHN" "DIV" "FCTR" "IRR?" "P2" "P3" "P4" "PL" "PPRD" "PR?" |
Section Header
Creating a Matrix Copying a Matrix Storing the Elements of a Matrix & Control Number Random Matrices Random Symmetric Matrices Transpose of a Matrix Transpose of a Symmetric Matrix Section Header Exponential of a Matrix Pseudo-Inverse of a Matrix ( Gréville's Method ) Logarithm of a Matrix Subtraction of 2 Matrices Addition of 2 Matrices Multiplication of 2 Matrices Inverse of a Matrix Euclidean Norm of a Matrix Power of a Matrix P-th Root of a Matrix Matrix Polynomial Square Root of a Matrix Product of a Matrix Transpose by another Matrix Trace of a Square Matrix Section Header Deflation Method Jacobi's Method Jacobi's Method for Complex Matrices Characteristic Polynomial Minimal Polynomial Section Header A Subroutine used by "D5" Determinant of order 3 Determinant of order 4 Determinant of order 5 Determinant of order n Linear Systems Linear Systems ( variant ) Underdetermined Linear Systems Overdetermined Linear Systems Section Header 4-Dimensional Magic Hypercubes of odd orders Magic Cubes Magic Squares Pandiagonal & Panmagic Cubes Panmagic Squares Section Header Bairstow's Method Cohn's Irreducibility Criterion Polynomial Euclidean Division Factorization over Z[X] A Better Irreducibility Criterion Quadratic Equations Cubic Equations Quartic Equations Evaluation of a Polynomial P(x) Product of 2 Polynomials Primality Test ( Non-negative integers ) |
Data Registers: • R00 = function name ( Register R00 is to be initialized before executing "CRMAT" )
and the registers containing the matrix elements ( no other register is
used )
Flags: /
Subroutine: The
function f ( assuming i is in X-register and j is
in Y-register upon entry )
If needed, you have also
1 in L-register,
n = the number of rows in
N-register and
k = the pointer value of
the current element in M-register ( from 1 to n.m , counted in column order
)
| STACK | INPUTS | OUTPUTS |
| X | bbb.eeerr | / |
where bbb.eeerr is the control number of the array. rr = number of rows.
Example: Store the elements of the 3x4 matrix [ai,j] defined by ai,j = i2 + j3 into registers R01 thru R12.
01 LBL "T" defines the function
f(i,j) = i2 + j3
02 X^2
03 X<>Y
04 3
05 Y^X
06 +
07 RTN
"T" ASTO 00
1.01203 XEQ "CRMAT" and it yields
2 9 28
65
R01 R04 R07 R10
5 12 31 68
= R02 R05 R08
R11 respectively
10 17 36 73
R03 R06 R09 R12
Note:
-If your matrix is actually a vector ( only one row ) you can
key in bbb.eee instead of bbb.eee01
-You want to copy the elements of a matrix ( control number = bbb.eeerr1
) into a block of registers that starts with Rbb2
| 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
Warning:
-The 2 blocks of registers cannot overlap unless bbb2
< bbb1
Storing the Elements of a Matrix & Control Number
| 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
-"RANM" stores pseudo-random integers ( between 1 and N ) into registers
Rbb , ........ , Ree
-No other register is used.
| STACK | INPUTS | OUTPUTS |
| Y | bbb.eeerr | / |
| X | N | / |
where bbb.eeerr is the control number of the array
Example: Store random integers between 1 and 12 into registers R01 , R02 , ........ , R07
1.007 ENTER^
12 XEQ "RANM" gave
R01 = 10 R02 = 5 R03 = 8 R04 = 6
R05 = 6 R06 = 8 R07 = 10 ... when I
did it !
-Since the date & time are used to initialize the random number
generator, you will probably get different values...
-"RANSYM" directly stores pseudo-random integers ( between 1 and N )
in registers R01 thru Rn2
-No other register is used.
| STACK | INPUTS | OUTPUTS |
| Y | n | / |
| X | N | / |
Example: Store a random symmetric matrix of order 5 - elements between 1 and 12 - in registers R01 thru R25
5 ENTER^
12 XEQ "RANSYM" gave ( when I executed "RANSYM"
)
R01 R06 R11 R16
R21
8 2 12 9
11
R02 R07 R12 R17
R22
2 12 11 8
11
R03 R08 R13 R18
R23 =
12 11 12 10 6
R04 R09 R14 R19
R24
9 8 10 12
6
R05 R10 R15 R20
R25
11 11 6 6
4
-Like "RANM", "RANSYM" uses the date & time to initialize the random
number generator, so you will probably get different values...
-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.
| 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
Transpose of a Symmetric Matrix
-Unlike "TRN" , "TRN2" replaces a square matrix by its
transpose.
| 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
-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 )
| 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)
Pseudo-Inverse of a Matrix ( Gréville'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*"
| 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
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
-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 )
| 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)
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
"M-" calculates the difference between 2 matrices A & B and
returns the control number of A-B
| 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 .....
1.00602 ENTER^
7.01202 ENTER^
15
XEQ "M-" >>>> 15.02002 and
1 -6 3 R15
R17 R19
A-B = -7 3 1
= R16 R18 R20 respectively
-The destination block of registers may be the same as the one of matrix
A or B.
"M+" calculates the sum of 2 matrices A & B and returns the
control number of A+B
| 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 .....
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
-The destination block of registers may be the same as the one of matrix
A or B.
"M*" calculates the product of 2 matrices A & B and returns
the control number of A*B
| 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.
"MINV" can invert up to a 17x17 matrix.
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 "MINV"
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.
| STACK | INPUTS | OUTPUTS |
| X | n | det A |
where n is the order of the matrix A
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
5 XEQ "MINV" 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
-5 19 -27
17 -4
R09 R14 R19 R24 R29
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 !
- "MNORM" computes || A || = ( SUMi,j a2i,j
)1/2
| 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 "MORM" >>>> || A ||
= 14.8996643
-"MPOW" 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
| STACK | INPUTS | OUTPUTS |
| Y | bbb.eeerr1 | / |
| X | p > 0 | bbb.eeerr2 |
( bbb > 003 )
Example:
1 4
9 R11 R14
R17
A = 3 5 7 =
R12 R15 R18 and you want to compute
A^7
2 1
8 R13 R16
R19
11.01903 ENTER^
7
XEQ "MPOW" >>>> 20.02803 = the control number of A7
7851276 8652584 31076204
R20 R23 R26
A7 = 8911228 9823060
35267932 = R21 R24
R27
5829472 6422156 23076808
R22 R25 R28
-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 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*"
& "LS3"
| 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 "MROOT" 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
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
1 4 2 9
-0.224295471 0.139205813
0.171840655 1.642919190
-The turbo mode is recommended...
-Given a polynomial p(x) = am xm + am-1
xm-1 + ............... + a0
and a nxn matrix A,
"PMAT" calculates p(A) = am Am +
am-1 Am-1 + ............... + a0 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 = am ....... • Ree = a0
• Rbb' thru • Ree' = the coefficients
of the matrix A, in column order.
the 2n2 registers from RBB: temp
Flag: F07
Subroutines: "MCO" "M*"
| 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
3548 1887 4705
R24 R27 R30
p(A) = 3727
1987 4962 = R25
R28 R31
2539 1351 3385
R26 R29 R32
-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*" "LS3"
| 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" >>>> X = 1.00903 = control number and we get:
1.794981016 0.656367529 0.540425261
sqrt(A) = 0.772143191
1.061374174 1.582989342
0.501888909 0.231634627 1.833600670
Product of a Matrix Transpose by another Matrix
-It is sometimes useful to multiply the transpose of a matrix A by another
matrix B directly. ( product
tA.B )
| 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
-The trace of a square matrix equals the sum of its diagonal elements.
| 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
-If k1 is an eigenvalue , U1 the corresponding
unit eigenvector of A and U1' the corresponding unit eigenvector
of tA
then A and the new matrix A' = A - k1
( U1 tU1' ) / ( U1.U1'
) have the same eigenvalues and eigenvectors , except that
k1 has been replaced by 0.
-Thus, the 2nd dominant eigenvalue of A is now the 1st
dominant eigenvalue of A'
-"DFL" uses this method to calculate all the eigenvalues and eigenvectors , provided they are real and verify | k1 | > | k2 | > ....... > | kn | > 0
1°) "DFL" calculates k1 and U1 ( with flag F02 clear ) and the program stops ( line 87 )
>>>> k1
is in X-register
>>>> U1
is in registers R01 thru Rnn
2°) If you press R/S , the HP-41 sets flag F02 and
computes k1 once again & U1'
the matrix A is replaced
by A' in registers Rn+1 thru Rn2+n and when the program
stops again:
>>>> k2
is in X-register
>>>> U2
is in registers R01 thru Rnn
... etc ...
Data Registers:
• R00 = n
• R01 = a11
• Rn2-n+1 = a1n
...............................................
INPUTS
• Rnn = an1
• Rn2 = ann
---------
R00 = n
R01 = u1 Rn+1 =
a'11 Rn2+1
= a'1n
.............
A will be replaced with
k in X-register. U ( u1 , ........
, un ) = a unit eigenvector
OUTPUTS
A' = A - k ( U.tU' ) / ( U.U' ) when F02 is set.
Rnn = un
R2n = a'n1
Rn2+n = a'nn
Flag: F02
Subroutine: "M*"
| STACK | INPUTS | OUTPUTS |
| X | / | k1 , k2 , ... |
Example: Find all the eigenvalues and the
corresponding unit eigenvectors of the matrix
1 2 4
A = 4 3 5
7 4 7
1°) n = 3 therefore 3 STO 00
Store the 9 numbers
1 2 4
R01 R04 R07
4 3 5 in
R02 R05 R08 respectively.
7 4 7
R03 R06 R09
XEQ "DFL" the HP-41 displays the successive k1-approximations ( with F02 clear )
- Eventually, k1 = 12.90692995 is in X-register and U1 is in registers R01 R02 R03 : U1 ( 0.348663346 ; 0.530674468 ; 0.772540278 )
2°) To obtain the second eigenvalue simply press R/S
the HP-41 sets flag F02 and displays the successive
k1-approximations converging to 12.90692995 again
then flag F02 is cleared and the HP-41 displays
the successive k2-approximations and
k2 = -2.092097594 in X-register and U2 in R01R02R03 : U2 ( 0.800454175 ; -0.041651079 ; -0.597945065 )
3°) Once again R/S ( similar displays ) and finally,
k3 = 0.185167649 in X-register and U3 in R01R02R03 : U3 ( -0.094824730 ; 0.897989404 ; -0.429678136 )
Notes:
-"DFL" uses a random vector as first guess for the eigenvector, so you
may get slightly different results than above.
-If k1 and k2 are real but
k3 ... are complex, the first 2 eigenvalues and
eigenvectors will be correctly computed, but not the other ones.
"JCB" computes the eigenvalues & eigenvectors of a matrix A, provided all the eigenvalues are real ( and distinct if A is not symmetric ).
*** A non-symmetric real matrix may have complex eigenvalues, but all the eigenvalues of a real symmetric matrix A are necessarily real !
-In the Jacobi's method, the eigenvalues are the elements of a diagonal
matrix P equal to the infinite product ...Ok-1.Ok-1-1....O2-1.O1-1.A.O1.O2....Ok-1.Ok....
where the Ok are rotation matrices. The eigenvectors
are the columns of O1.O2....Ok-1.Ok....
-"JCB" finds the greatest off-diagonal element ai j
( lines 22 to 46 )
-Then, O1 is determined so that it makes a pair
of off-diagonal elements zero in O1-1.A.O1
-It happens if the rotation angle µ is choosen so that
Tan 2µ = 2.ai j / ( ai i - aj j
)
-The process is repeated until the greatest off-diagonal element is
smaller than E-8 in absolute value ( line 48 )
-The successive greatest | ai j | ( i > j ) are displayed
( line 47 ) when the routine is running.
*** Though it's not really orthodox, we can try to apply the same program to non-symmetric matrices: of course, it cannot work if some eigenvalues are complex.
But if all the eigenvalues are real, the infinite product
T = [ ti j ] = ...Ok-1.Ok-1-1....O2-1.O1-1.A.O1.O2....Ok-1.Ok....
is an upper triangular matrix,
the diagonal elements of T ( i-e ti i = ki
) are the n eigenvalues and the first column of U = O1.O2....Ok-1.Ok....
is an eigenvector corresponding to k1
-If all the eigenvalues are distinct, we can create an upper triangular matrix W = [ wi j ] defined by
wi j = 0 if i > j
wi i = 1
wi j = Sumk = i+1, i+2 , .... , j
ti k wk j / ( tj j - ti i )
if i < j
-Then, the n columns of U.W are the eigenvectors corresponding to the n eigenvalues t11 = k1 , ............... , tnn = kn
*** When "JCB" stops:
the n eigenvalues are stored in registers R01 thru
Rnn
the first eigenvector is stored in registers Rn+1 thru R2n
..................................................................................
the nth eigenvector is stored in registers Rn2+1
thru Rn2+n as shown below.
Data Registers:
• R00 = n
• R01 = a11 .................
• Rn2-n+1 = a1n
..................................................................
( the n2 elements of A in column order )
INPUTS
• Rnn = an1 .................
• Rn2 = ann
---------
DURING
R01 thru Rn2 = the "infinite" product ...Ok-1.Ok-1-1....O2-1.O1-1.A.O1.O2....Ok-1.Ok....
THE
Rn2+1 thru R2n2 = O1.O2........Ok.....
ITERATIONS R2n2+1 thru
R2n2+n: temp ( for non-symmetric matrices only )
---------
R00 = n
R01 = k1 Rn+1 = V1(1)
R2n+1 = V1(2) ....................
Rn2+1 = V1(n)
.............
OUTPUTS
Rnn = kn R2n =
Vn(1)
R3n = Vn(2) .....................
Rn2+n = Vn(n)
( n eigenvalues ; 1st eigenvector ; 2nd eigenvector ;
.................. ; nth eigenvector )
Flag: F02
Set flag F02 for a symmetric matrix
Clear flag F02 for a non-symmetric matrix
Subroutines: /
( SIZE 2n2+1 if A is
symmetric / SIZE 2n2+n+1 if A is not symmetric )
| STACK | INPUTS | OUTPUTS |
| X | / | k1 |
Example1: Let's compute all the
eigenvalues and the eigenvectors of the matrix
1
2 4
A = 2 7 3
4 3 9
n = 3 therefore 3 STO 00 and we store the 9 numbers
1 2 4
R01 R04 R07
2 7 3
into R02 R05 R08
respectively.
4 3 9
R03 R06 R09
SF 02 ( A is symmetric ) XEQ "JCB" yields
k1 = 12.81993499 in R01
k2 = 4.910741214 in R02
k3 = -0.730676199 in R03
and the 3 corresponding eigenvectors:
V1 ( 0.351369026 ; 0.521535689 ; 0.777521917
) in ( R04 ; R05 ; R06 )
V2 ( -0.101146468 ; 0.846760701 ; -0.522269766
) in ( R07 ; R08 ; R09 )
V3 ( -0.930757326 ; 0.104865823 ; 0.350276976
) in ( R10 ; R11 ; R12 )
Example2:
1 2 4
A = 4 3 5
7 4 7
n = 3 STO 00 and we store the 9 numbers
1 2 4
R01 R04 R07
4 3 5
into R02 R05 R08
respectively.
7 4 7
R03 R06 R09
CF 02 ( A is not symmetric )
XEQ "JCB" >>>> k1 = 12.90692994 and we have
k1 =
12.90692994 in R01
k2 =
0.185167649 in R02
k3 = -2.092097593
in R03
and the 3 corresponding eigenvectors:
V1
( 0.348663347 ; 0.530674467 ; 0.772540277
) in ( R04 ; R05 ; R06 )
V2
( -0.095420104 ; 0.903627529 ; -0.432375917 )
in ( R07 ; R08 ; R09 )
V3
( -0.830063031 ; 0.043191757 ; 0.620063097 )
in ( R10 ; R11 ; R12 )
-The eigenvectors are not unit vectors except the first
one.
-Divide V2 & V3 by || V2
|| & || V3 || respectively if you want unit eigenvectors.
Notes:
-If A is non-symmetric, this program only works if all the eigenvalues
are real and distinct.
-If all the eigenvalues are real but not distinct, they will be correctly
computed but not the eigenvectors ( DATA ERROR line 219 )
Jacobi's Method for Complex Matrices
-"JCBZ" uses a variant of the Jacobi's algorithm to compute the eigenvalues & eigenvectors of a complex matrix:
*** The eigenvalues of A are the diagonal-elements of an upper triangular
matrix T equal to the infinite product ...Uk-1.Uk-1-1....U2-1.U1-1.A.U1.U2....Uk-1.Uk....
where the Uk are unitary matrices. The eigenvectors
are the columns of U = U1.U2....Uk-1.Uk....
if A is Hermitian ( i-e if A equals its transconjugate )
Actually, T is diagonal if A is Hermitian.
-"JCBZ" finds the greatest element ai j below the main diagonal
( lines 23 to 60 )
-Then, U1 is determined so that it places a
zero in position ( i , j ) in U1-1.A.U1
U1 has the same elements as the Identity matrix, except that ui i = uj j = x and ui j = y + i.z , uj i = -y + i.z
with y + i.z = C.x , x = ( 1 + | C |2 ) -1/2 and C = 2.ai j / [ ai i - ai j + ( ( ai i - ai j )2 + 4 ai j aj i )1/2 ] ( line 156 R-P avoids a test if the denominator = 0 )
-The process is repeated until the greatest sub-diagonal element is smaller than E-8 in magnitude ( line 62 )
-The successive greatest | ai j | ( i > j ) are displayed ( line 61 ) when the routine is running.
*** If the matrix is non-Hermitian, and if all the eigenvalues are distinct, we define an upper triangular matrix W = [ wi j ] by
wi j = 0 if i > j
wi i = 1
wi j = Sumk = i+1, i+2 , .... , j
ti k wk j / ( tj j - ti i )
if i < j
-Then, the n columns of U.W are the eigenvectors corresponding
to the n eigenvalues t11 = k1 , ...............
, tnn = kn
Data Registers:
• R00 = n
• R01,R02 = a11 .................
• R2n2-2n+1,R2n2-2n+2 = a1n
.................................................................................................
( the n2 elements of A in column order )
INPUTS
• R2n-1,R2n = an1 .................
• R2n2-1,R2n2 = ann
--------- odd
registers = real part of the coefficients , even registers = imaginary
part of the coefficients
DURING
R01 thru R2n2 = the "infinite" product ...Uk-1.Uk-1-1....U2-1.U1-1.A.U1.U2....Uk-1.Uk....
THE
R2n2+1 thru R4n2 = U1.U2........Uk.....
ITERATIONS R4n2+1 thru
R4n2+2n: temp ( for non-Hermitian matrices only )
---------
R00 = n
R01,R02 = k1
R2n+1,R2n+2 = V1(1)
.................... R2n2+1,R2n2+2
= V1(n)
.......................
........................
.................... ..................................
OUTPUTS
R2n-1,R2n = kn R4n-1,R4n
= Vn(1)
..................... R2n2+2n-1,R2n2+2n
= Vn(n)
( n eigenvalues ;
1st eigenvector ; ..................................
; nth eigenvector )
Flag: F02
Set flag F02 for a Hermitian matrix
Clear flag F02 for a non-Hermitian matrix
Subroutines: /
( SIZE 4n2+1 if A is Hermitian / SIZE
4n2+2n+1 if A is non-Hermitian )
| STACK | INPUTS | OUTPUTS |
| Y | / | y1 |
| X | / | x1 |
where k1 = x1 + y1 = the first eigenvalue.
Example1:
1 4 -7.i 3-4.i
A = 4+7.i 6
1-5.i
A is equal to its transconjugate so A is Hermitian
3+4.i
1+5.i 7
1 0 4 -7 3 -4
R01 R02 R07 R08 R13 R14
Store 4 7 6 0
1 -5 into R03 R04 R09 R10
R15 R16 respectively and
n = 3 STO 00
3 4 1 5 7 0
R05 R06 R11 R12 R17 R18
SF 02 ( the matrix is Hermitian ) XEQ "JCBZ" yields
k1
= 15.61385271 + 0.i = ( X , Y ) = ( R01 , R02 )
k2
= 5.230678474 + 0.i = ( R03 , R04 )
k3
= -6.844531162 + 0.i = ( R05 , R06 ) and the 3
corresponding eigenvectors
V1( 0.481585119 + 0.108201123
i , 0.393148652 + 0.527305194 i , -0.142958847 + 0.550739892
i ) = ( R07 R08 , R09 R10 , R11 R12 )
V2( -0.436763638 + 0.075843695
i , 0.210271354 - 0.463555612 i , -0.516799072 + 0.526598642
i ) = ( R13 R14 , R15 R16 , R17 R18 )
V3( -0.706095475 + 0.247553486
i , 0.326530385 + 0.449069516 i , 0.363126603 - 0.
i )
= ( R19 R20 , R21 R22 , R23 R24 )
-Due to roundoff errors, registers R02 R04 R06 contain very small numbers
instead of 0
but actually, the eigenvalues of a Hermitian matrix are always
real.
Example2:
1+2.i
2+5.i 4+7.i
A = 4+7.i 3+6.i 3+4.i
A is non-Hermitian
3+4.i
1+7.i 2+4.i
1 2 2 5 4 7
R01 R02 R07 R08 R13 R14
Store 4 7 3 6
3 4 into R03 R04 R09 R10
R15 R16 respectively and
n = 3 STO 00
3 4 1 7 2 4
R05 R06 R11 R12 R17 R18
CF 02 ( the matrix is non-Hermitian ) XEQ "JCBZ" produces
k1
= 7.656606009 + 15.61073835 i = ( X , Y ) = ( R01 ,
R02 )
k2
= 1.661248138 - 1.507335315 i = ( R03 , R04 )
k3
= -3.317854151 - 2.103403073 i = ( R05 , R06 )
and the 3 corresponding eigenvectors
V1( 0.523491429 + 0.008555748
i , 0.650242685 - 0.046353000 i , 0.546288477 + 0.049882590
i ) = ( R07 R08 , R09 R10 , R11 R12
)
V2( -0.4777850326 - 0.196368365
i , 0.660547554 - 0.265888145 i , -0.356842635 + 0.318267519
i ) = ( R13 R14 , R15 R16 , R17 R18 )
V3( -0.567221101 - 0.574734664
i , 0.141603481 + 0.549379028 i , 0.480533312 - 0.090337847
i ) = ( R19 R20 , R21 R22 , R23 R24 )
-If there exist a number k and a non-zero vector V such that A.V
= k.V k is called an eigenvalue and V an eigenvector
of the matrix A.
-The characteristic polynomial p(x) = xn + c1.xn-1
+ c2.xn-2 + ......... + cn-2.x2
+ cn-1.x + cn is defined by p(x) = det ( x.I - A
) where I is the n x n unit matrix.
thus, its roots are actually the n eigenvalues of A.
-"PCAR" uses the following formulae: we define a sequence of matrices Mk by:
M0 = I and
Mk+1 = Mk.A - trace(Mk.A) I /k
( the trace of a matrix is the sum of its diagonal elements )
then we have:
ck = - trace(Mk-1.A)/k
-The elements of A are to be stored into contiguous registers from R01
to Rn2 , n being stored in R00
-When the program stops, R01 = c1 , R02 = c2
, ............... , Rnn = cn
( c1 = -trace(A) and cn = (-1)n
det(A) )
Data Registers:
• R00 = n
• R01 = a11 ......................
• Rn2-n+1 = a1n
...................................................................
( the n2 elements of A )
INPUTS
• Rnn = an1 ......................
• Rn2 = ann
---------
R00 = n
R01 = c1 Rn+1 =
a11 Rn2+1
= a1n
.............
..........................................
OUTPUTS
Rnn = cn
R2n = an1
Rn2+n = ann
Flag: F10
Subroutine: "M*"
( SIZE 3n2+n+1 )
| STACK | INPUT | OUTPUT |
| X | / | 1.n+1 |
Example: Find
the characteristic polynomial of the matrix
1 2 4
A = 4 3 5
7 4 7
-First we store these 9 elements in registers R01 thru R09
R01 R04 R07
1 2 4
R02 R05 R08
= 4 3 5 respectively
and n = 3 STO 00
R03 R06 R09
7 4 7
-Then XEQ "PCAR" >>>> 1.004 = control number of the characteristic polynomial.
R01 = 1 , R02 = -11 , R03 = -25 , R04 = 5
Thus p(x) = x3 - 11 x2 - 25 x + 5
Notes:
-A is saved in registers Rn+2 thru Rn2+n+1 and n is saved
in R00
-The elements of A may be stored either in column order or in row order
since A and tA ( the transpose of A ) have the same characteristic
polynomial.
-Any root-finding program yields the 3 eigenvalues
k1 = 12.90692994 , k2
= -2.092097589 , k3 = 0.185167648
-Then, subtracting ( for instance ) k1 to the diagonal elements
of A leads to a singular system with an infinite set of solutions ( the
eigenvectors )
-One solution is ( 0.451320606 ; 0.686921424 ;
1 ) which is an eigenvector corresponding to the first
eigenvalue.
-This method could be used to obtain all the eigenvalues and eigenvectors
of a n x n matrix A.
-The characteristic polynomial P verifies P(A) = 0 ( Cayley-Hamilton
theorem )
-The minimal polynomial p is the unitary least degree
polynomial that satisfies p(A) = 0.
-It is a divisor of the characteristic polynomial.
-We start with a random vector V and we compute AV , A2V
, ........... , AnV
-Then, "PMIN" seeks the smallest k such that V , AV , A2V
, ........... , AkV are linearly dependent.
-A relation a0 V + a1 AV + ........... +
ak-1 Ak-1V + AkV = 0 is found, whence
p(x)
Data Registers: • R00 = n ( Registers R00 thru Rn2 are to be initialized before executing "PMIN" )
• R01 = a11 ......................
• Rn2-n+1 = a1n
...................................................................
( the n2 elements of A )
• Rnn = an1 ...................... • Rn2 = ann Rn2+1 to R2n2+n: temp
>>>> When the program stops, R01 thru Ree = the coefficients of p(x)
Flags: /
Subroutines: "M*" "MCO"
"LS3"
( 164 bytes / SIZE 2n2+n+1 )
| STACK | INPUT | OUTPUT |
| X | / | 1.eee |
Example:
Find the minimal polynomial of the matrix
[ [
3 -1 1 ]
A = [ -1 3 1
]
[ 1 1 3 ] ]
-Store
3 -1 1
R01 R04 R07
-1 3 1
into R02 R05 R08
respectively and
n = 3 STO 00
1 1 3
R03 R06 R09
XEQ "PMIN" >>>> 1.003 --- Execution time = 36s --- and we get R01 = 1 R02 = -5 R03 = 4 ( with some roundoff-errors )
-So, p(x) = x2 - 5 x + 4
( the characteristic polynomial is P(x) = x3 - 9 x2
+ 24 x - 16 = ( x - 4 ) p(x) )
Notes:
-Since we use a pseudo-random vector V, there is a small risk of getting
an eigenvector, which would lead to a wrong result.
-In practice however, that seems highly improbable.
-If the coefficients of A are integers, the coefficients of p are integers
too.
Data Registers: R00 = unused ( Registers R01 thru R09 are to be initialized before executing "D3" )
• R01 = a11 • R04 = a12
• R07 = a13
• R02 = a21 • R05 = a22
• R08 = a23
• R03 = a31 • R06 = a32
• R09 = a33
Flags: /
Subroutines: /
| STACK | INPUT | OUTPUT |
| X | / | Deteminant |
Example: Calculate
|
4 9 2 | Store these
R01 R04 R07
D = | 3 5 7 |
9 numbers R02 R05 R08
respectively
| 8 1 6 |
into
R03 R06 R09
and XEQ "D3" >>>> Det = 360
-Determinants of order 4 are computed by a sum of products of determinants
of order 2.
Data Registers: R00 = determinant at the end ( Registers R01 thru R16 are to be initialized before executing "D4" )
• R01 = a11 • R05 = a12
• R09 = a13 • R13 = a14
• R02 = a21 • R06 = a22
• R10 = a23 • R14 = a24
• R03 = a31 • R07 = a32
• R11 = a33 • R15 = a34
• R04 = a41 • R08 = a42
• R12 = a43 • R16 = a44
Flags: /
Subroutines: /
| STACK | INPUT | OUTPUT |
| X | / | Deteminant |
Example: Calculate
|
1 8 13 12 |
Store these R01
R05 R09 R13
D = | 14 11 2
7 | 16 numbers
R02 R06 R10 R14
respectively
|
4 5 16 9 |
into
R03 R07 R11 R15
|
15 10 3 6 |
R04 R08 R12 R16
XEQ "D4" >>>> Det = 0
-Determinants of order 5 are computed by a sum of products of determinants
of order 2 by determinants of order 3.
-Registers R00 to R25 are temporarily disturbed during the calculations,
but their contents are finally restored.
Data Registers: R00 = det A at the end ( Registers R01 thru R25 are to be initialized before executing "D5" )
• R01 = a11 • R06 = a12
• R11 = a13 • R16 = a14
• R21 = a15
• R02 = a21 • R07 = a22
• R12 = a23 • R17 = a24
• R22 = a25
• R03 = a31 • R08 = a32
• R13 = a33 • R18 = a34
• R23 = a35
• R04 = a41 • R09 = a42
• R14 = a43 • R19 = a44
• R24 = a45
• R05 = a51 • R10 = a52
• R15 = a53 • R20 = a54
• R25 = a55
Flags: /
Subroutine: "D0"
| STACK | INPUT | OUTPUT |
| X | / | Deteminant |
Example: Calculate
|
1 7 13 19 25
|
R01 R06 R11 R16 R21
|
14 20 21 2 8
|
R02 R07 R12 R17 R22
D = | 22 3 9
15 16 |
= R03 R08
R13 R18 R23
respectively
|
10 11 17 23 4 |
R04 R09 R14 R19 R24
|
18 24 5 6 12
|
R05 R10 R15 R20 R25
XEQ "D5" >>>> Det = -4680000
"DET" simply uses "LS3" to compute the determinant of a
square matrix of order n.
Data Registers: R00 = det A at the end ( Registers R01 thru Rn2 are to be initialized before executing "DET" )
• R01 = a11 • Rn+1 = a12
..................... • Rn2-n+1 = a1n
• R02 = a21 • Rn+2 = a22
..................... • Rn2-n+2 = a2n
........................................................................................
• Rnn = an1 • R2n = an2 ..................... • Rn2 = ann
Flags: CF 00 = a pivot p is regarded
as zero if | p | < 10-7
CF 01 = partial pivoting
SF 00 = a pivot p is regarded as zero if p = 0
SF 01 = no pivoting
Subroutines: /
| STACK | INPUT | OUTPUT |
| X | n | Deteminant |
where n is the order of the square matrix.
Example: Calculate
|
4 9 2 |
R01 R04 R07
D = | 3 5 7 |
into R02 R05 R08
respectively
| 8 1 6 |
R03 R06 R09
3 XEQ "DET" >>>> Det
= 360
"LS" allows you to solve linear systems, including overdetermined
and underdetermined systems.
You can also invert up to a 12x12 matrix.
It is essentially equivalent to the RREF function of the HP-48
( a "little" slower, I grant you that! ):
Its objective is to reduce the matrix on the upper left to a
diagonal form with only ones in the diagonal.
The determinant of this matrix is also computed and stored in
register R00.
( if there are more rows than columns, R00 is not always equal
to the determinant of the upper left matrix because of rows exchanges )
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:
It's sometimes useful for matrices like Pascal's matrices of
high order.
They are extremely troublesome and many roundoff errors can
occur.
But if you set flag F01, all the coefficients will be computed
with no roundoff error at all, because all the pivots will be ones!
One advantage of this program is that you can choose the beginning
data register - except R00 - ( this feature is used to solve non-linear
systems too ):
You store the first coefficient into Rbb , ... , up to the last
one into Ree ( column by column )
( bb > 00 )
Then you key in bbb.eeerr ENTER^
p XQ "LS"
and the system will be solved. ( bbb.eeerr ends
up in L-register )
where r is the number of rows of the matrix
and p is a small non-negative number that you choose
to determine tiny elements: if a pivot is smaller than p it will be considered
to be zero.
Don't interrupt "LS" because registers P and Q are used ( there
is no risk of crash, but their contents could be altered )
| STACK | INPUTS | OUTPUTS |
| Y | bbb.eeerr | 1 |
| X | p | determinant |
| L | / | bbb.eeerr |
Example1: Find the solution of
the system:
2x + 3y + 5z + 4t = 39
-4x + 2y + z + 3t = 15
3x - y + 2z + 3t = 19
5x + 7y - 3z + 2t = 18
If you choose to store the first element into R01, you have to store these 20 numbers:
2 3
5 4 39
R01 R05 R09 R13 R17
-4 2 1
3 15 in
R02 R06 R10 R14 R18
respectively
3 -1 2
3 19
R03 R07 R11 R15 R19
5 7 -3
2 18
R04 R08 R12 R16 R20
The first register is R01, the last register is R20 and there are 4
rows, therefore the control number of the matrix is 1.02004
If you choose p = 10-7 key in:
CF 01
1.02004 ENTER^
E-7 XEQ "LS"
and you obtain 840 in X-register and in R00: it is the
determinant of the 4x4 matrix on the left.
Registers R01 thru R16 now contains the unit matrix and registers R17
thru R20 contains the solution x = 1 ; y = 2 ; z = 3 ; t = 4
Thus, the array has been changed into:
1 0 0
0 1
0 1 0
0 2
the solution is the last column
0 0 1
0 3
of the new matrix.
0 0 0
1 4
Example2: Solve the system:
5x + y + z
= 8
4x - y + 2z
= 13
x + 2y
- z = -5
7x - 4y + 5z = 31
-Suppose we need to store the first element into R11 , then we store these 16 numbers
5
1 1 8
R11 R15 R19 R23
4 -1
2 13 in
R12 R16 R20 R24
respectively
1
2 -1 -5
R13 R17 R21 R25
7 -4
5 31
R14 R18 R22 R26
Here, the control number of this array is 11.02604 and if we take p = 10-7 again,
11.02604 ENTER^
E-7 XEQ "LS" gives 0 in
X-register and in R00 = the determinant of the 4x4 matrix.
and registers R11 thru R27 are now:
1 0 0.333333333
2.333333333
0 1 -0.666666667
-3.666666669
0 0
5.10-10
-2.10-9
0 0
0
4.10-9
Thus, the system is equivalent to:
x + z /3 = 7/3
and there is an infinite set of solutions:
y - 2z /3 = -11/3
z may be chosen at random and x and y are determined by
0 = 0
x = -z /3 + 7/3
0 = 0
y = 2z /3 - 11/3
Note: If the last number of the initial matrix is 32 instead of 31 the array is changed into:
1 0 0.333333333
0 i-e
x + z /3 = 0
0 1 -0.666666667
0
y - 2z /3 = 0
0 0
5.10-10 0
0 = 0
0 0
0
1
0 = 1 and there is no solution
at all !
Example3: Invert Pascal's matrix:
1 1 1 1
1
1 2 3 4
5
1 3 6 10 15
1 4 10 20 35
1 5 15 35 70
This is equivalent to solve 5 linear systems simultaneously. If you begin at register R01, you have to store the 50 numbers
1 1 1
1 1 1 0 0 0 0
R01 R06 R11 R16 R21 R26 R31 R36
R41 R46
1 2 3
4 5 0 1 0 0 0
R02 R07 R12 R17 R22 R27 R32 R37
R42 R47
1 3 6
10 15 0 0 1 0 0
in R03 R08 R13 R18
R23 R28 R33 R38 R43 R48
( the right half is the unit matrix )
1 4 10 20 35
0 0 0 1 0
R04 R09 R14 R19 R24 R29 R34 R39
R44 R49
1 5 15 35 70
0 0 0 0 1
R05 R10 R15 R20 R25 R30 R35 R40
R45 R50
The control number of this rectangular array is 1.05005 and we can choose p = 0 because Pascal's matrix is non-singular:
1.05005 ENTER^
0
XEQ "LS" yields 1 in X-register and in R00 ( the determinant
of Pascal's matrices is always 1 )
The array is now:
1 0 0 0
0 5 -10 10 -5 1
0 1 0 0
0 -10 30 -35 19 -4
0 0 1 0
0 10 -35 46 -27 6
( the unit matrix is now on the left and the inverse matrix is on the right,
0 0 0 1
0 -5 19 -27 17 -4
in registers R21 thru R50 )
0 0 0 0
1 1 -4 6 -4
1
-This program can be used to invert a n x n matrix, but it requires
2n2 registers and thus, it cannot invert a 13x13 matrix.
-Use "MINV" if you want to find the inverse a 17x17 matrix.
-In this variant, the matrix is the right-part of the array so that,
when the program stops,
the solution ( x1 , x2 , .............
, xn ) is in registers R01 R02 ..............
Rnn
-Here, the attempt to diagonalization starts by the lower right corner
of the matrix.
-Flags F00 & F01 play the same role as above.
( SIZE n.m+1 )
| STACK | INPUTS | OUTPUTS |
| T | / | n.nnn |
| Z | / | / |
| Y | n | / |
| X | m | det A |
n = number of rows
m = number of columns
T-output is useful to retrieve n
Example:
2x + 3y + 5z
+ 4t = 39
39 = 2x + 3y + 5z + 4t
-4x + 2y + z + 3t
= 15 is re-written
15 = -4x + 2y + z + 3t
3x -
y + 2z + 3t = 19
19 = 3x - y + 2z + 3t
5x + 7y - 3z
+ 2t = 18
18 = 5x + 7y - 3z + 2t
and we store these 20 numbers:
39
2 3 5 4
R01 R05 R09 R13 R17
15 -4
2 1 3 in
R02 R06 R10 R14 R18
respectively
19 3 -1
2 3
R03 R07 R11 R15 R19
18 5
7 -3 2
R04 R08 R12 R16 R20
-There are 4 rows and 5 columns,
CF 00 CF 01
4 ENTER^
5 XEQ "LS3" >>>> Det
A = 840 = R00
-Registers R05 thru R16 now contain the unit matrix and registers R01
thru R04 contain the solution x = 1 , y = 2 , z = 3 , t = 4
-Thus, the array has been changed into:
1
1 0 0 0
2
0 1 0 0
the solution is the first column
3
0 0 1 0
of the new matrix.
4
0 0 0 1
-When the program stops, R00 = det A
-If you have to invert a matrix, the inverse will be in registers R01
thru Rn2 at the end
Underdetermined Linear Systems
-"LS3-" computes the Moore-Penrose solution given by the formula: X = tA ( A.tA ) -1 B
( SIZE 2n.m-n+1 )
| STACK | INPUTS | OUTPUTS |
| Y | n | / |
| X | m | 1.rrr,rr |
n = number of rows
m = number of columns (
n < m-1 )
1.rrr,rr is the control number of the solution
( in fact r = m-1 )
Example: Let's find the Moore-Penrose solution of the system:
1 = 2x + 3y + 7z + 4t
1 2 3 7 4
R01 R04 R07 R10 R13
4 = 3x + 2y - 5z + 8t
store 4 3 2 -5 8
into R02 R05 R08 R11 R14
respectively
7 = 4x + 5y + 6z + t
7 4 5 6 1
R03 R06 R09 R12 R15
-There are 3 rows and 5 columns so:
3 ENTER^
5 XEQ "LS3-" >>>>
1.00404
-The solution is in registers R01 R02 R03 R04 namely: x = 1.095088162 , y = 1.075776659 , z = -0.389378674 , t = -0.422963896
-When the program stops, R00 = det ( A.tA ) [
in this example, R00 = 128628 ]
-If R00 = 0 ( suggesting A.tA is singular), the results
are probably meaningless.
-An overdetermined system A.X = B has often no solution
at all!
-But "LS3+" calculates the least-squares solution:
-It minimizes || A.X - B ||
Formula: X = ( tA.A ) -1
( tA.B )
( SIZE n.m+m2 -m+1 )
| STACK | INPUTS | OUTPUTS |
| Y | n | det |
| X | m | 1.rrr |
n = number of rows
m = number of columns
( n >= m )
1.rrr is the control number of the
solution ( in fact r = m-1 )
Example: The system:
8 = 5x + y + z
still has no solution.
13 = 4x - y + 2z
-5 = x + 2y - z
32 = 7x - 4y + 5z
-20 = 2x + 5y - 9z
-We store these 20 coefficients into registers R01 thru R20 ( in column order )
8
5 1 1
R01 R06 R11 R16
13 4
-1 2
R02 R07 R12 R17
-5 1
2 -1 =
R03 R08 R13 R18
respectively
32 7
-4 5
R04 R09 R14 R19
-20 2
5 -9
R05 R10 R15 R20
-There are 5 rows and 4 columns so
5 ENTER^
4 XEQ "LS3+" >>>> 1.003
-The "least-squares solution" is in registers R01 R02 R03 namely: x = 2.071207430 , y = -3.201238398 , z = 0.904024763
-When the program stops, Y = R00 = det ( tA.A )
[ In this example R00 = 55233 ]
-If det ( tA.A ) = 0 ( suggesting
tA.A
is singular ) , the results are probably meaningless.
4-Dimensional Magic Hypercubes of odd orders
-"MG4DHC" creates 4-D magic (hyper-) cubes of odd orders
n
-The magic constant = n ( n4+1 )/2
Data Registers: R00 = n R01 thru R06: temp
Flag: F00 SF 00 to display the
successive elements
CF 00 to store the elements into R01 to Rn^4 ( possible only if n
= 1 or 3 )
Subroutines: /
| STACK | INPUTS | OUTPUTS |
| X | n | / |
Where n is an odd integer > 0
Example: SF 00 3 XEQ "MG4DHC" >>>> N1=46 N2=59 N3=18 ................... N79=64 N80=23 N81=36
The whole hypercube is:
46 59 18
62 12 49
15 52 56
17 48 58
51 61 11
55 14 54
1st layer = cube n°1
60 16 47
10 50 63
53 57 13
08 39 76
42 79 02
73 05 45
78 07 38
01 41 81
44 75 04
2nd layer = cube n°2
37 77 09
80 03 40
06 43 74
69 25 29
19 32 72
35 66 22
28 68 27
71 21 31
24 34 65
3rd layer = cube n°3
26 30 67
33 70 20
64 23 36
-Here the magic constant is 123 and the central element is 41.
-There are 8 great diagonals:
36 + 41 + 46 = 123
67 + 41 + 15 = 123
22 + 41 + 60 = 123
29 + 41 + 53 = 123
35 + 41 + 47 = 123
26 + 41 + 56 = 123
64 + 41 + 18 = 123
69 + 41 + 13 = 123
-"MGCB" creates normal ( Andrews ) magic cubes, i-e their elements
are 1 2 3 ............... n3
-An Andrews magic cube is a cubic array where the sums of the elements
equal the magic constant n(n3+1)/2 in the 3 directions
and the 4 space diagonals.
Data Registers: R00 = n
>>> R01 thru Rn3 = the elements of the cube if flag F00 is clear
Flag: F00
SF 00 = The elements are successively displayed
CF 00 = The elements are stored into R01 R02
.............
Subroutines: /
| STACK | INPUTS | OUTPUTS |
| X | n # 2 | / |
Example: With
n = 3
-SF 00 3 XEQ "MGCB" the HP-41 displays N1=1 N2=17 N3=24 ................. N25=4 N26=11 N27=27
-If CF 00, we get in registers R01 thru R27
R01 R04 R07
01 23 18
R02 R05 R08 =
17 03 22
R03 R06 R09
24 16 02
R10 R13 R16
15 07 20
R11 R14 R17 = 19
14 09
R12 R15 R18
08 21 13
R19 R22 R25 26
12 04
R20 R23 R26 = 06 25
11
R21 R24 R27 10
05 27
-The magic constant = 42
-"MGSQ" creates normal magic squares of order n # 2
-The elements of a normal magic square are 1 , 2 , .............
, n2 and it has the following property:
All its row sums, column sums and both diagonal sums are equal
to the magic constant M = n.(n2 + 1)/2
Data Registers: R00 = n
>>> R01 thru Rn2 = the elements of the square if flag F00 is clear
Flags: F00
SF 00 = The elements are successively displayed
CF 00 = The elements are stored into R01 R02
.............
Subroutines: /
| STACK | INPUTS | OUTPUTS |
| X | n # 2 | / |
Example: n = 4
-SF 00 4 XEQ "MGSQ" the HP-41 displays successively N1=16 N2=5 N3=9 N4=4 .................... N13=13 N14=8 N15=12 N16=1
-The whole magic square is:
16 02 03 13
05 11 10 08
09 07 06 12
04 14 15 01
-Here, the magic constant = 34
-"PMGCB" creates a perfect pandiagonal magic cube if n = 7 & panmagic cubes if n is a prime > 10 ( look at line 02 )
-The cube is perfect pandiagonal if all the orthogonal sections are
panmagic squares.
-In a panmagic cube , all the orthogonal or diagonal sections - broken
or unbroken - are panmagic squares.
Data Registers: R00 = n
>>> R01 thru Rn3 = the elements of the cube if flag F00 is clear
Flag: F00
SF 00 = The elements are successively displayed
CF 00 = The elements are stored into R01 R02
.............
Subroutines: /
| STACK | INPUTS | OUTPUTS |
| X | n prime > 6 | / |
Example: n = 7
SF 00 7 XEQ "PMGCB" >>>> N1=1 N2=334 N3=268 N4=209 N5=192 N6=133 N7=67 ............... N343=343
Note:
-Since n > 6 , only the first elements will be stored if CF 00
-"PMGSQ" creates panmagic squares of order n where n mod 4 # 2 & n mod 3 # 0
-A panmagic square ( also called pandiagonal square ) is a magic square
where all the broken diagonal sums are also equal to the
magic constant M = n.(n2 + 1)/2
Data Registers: R00 = n
>>> R01 thru Rn2 = the elements of the square if flag F00 is clear
Flag: F00
SF 00 = The elements are successively displayed
CF 00 = The elements are stored into R01 R02
.............
Subroutines: /
| STACK | INPUTS | OUTPUTS |
| X | n | / |
where n mod 4 # 2 & n mod 3 # 0
Example: n = 5
-SF 00 5 XEQ "PMGSQ" the HP-41 displays N1=1 N2=7 N3=13 N4=19 N5=25 ................... N25=12
-The panmagic square is
01 14 22 10
18
07 20 03 11
24
13 21 09 17
05
19 02 15 23
06
25 08 16 04
12
-Magic Constant = 65
Note:
-There is no normal panmagic square of order 2 , 6 , 10 , ....
-"BRST" factorizes the polynomial p(x) = an.xn+an-1.xn-1+
... + a1.x+a0 into quadratic factors and solves
p(x) = 0 ( n > 1 )
-If deg(p) is odd, we have p(x) = (an.x+b).(x2+u1.x+v1)...............(x2+um.x+vm)
where m = (n-1)/2
-If deg(p) is even
p(x) = (anx2+u1.x+v1)(x2+u2.x+v2)..........(x2+um.x+vm)
where m = n/2
-The coefficients u and v are found by the Newton method for solving
2 simultaneous equations.
-Then p is divided by (x2+u.x+v) and u
& v are stored into Ree-1 & Ree respectively
-The process is repeated until all quadratic factors are found, and when the program stops:
If deg(p) is odd
Rbb = an , Rbb+1 = b , Rbb+2 = u1 , Rbb+3 = v1
, ........ , Ree-1 = um , Ree = vm
If deg(p) is even
Rbb = an , Rbb+1 = u1 , Rbb+2 = v1 , ............................
, Ree-1 = um , Ree = vm
-Then, press R/S to get the successive roots.
Data Registers: R00 to R08: temp. ( R06 = u ; R07 = v ; R08 = bbb.eee )
• Rbb = an • Rbb+1 = an-1 , .......
, • Ree = a0 (
These n+1 registers are to be initialized before executing "BRST")
bb > 08
Flag: F00
Subroutines: "DIV" & "P2"
| STACK | INPUT | OUTPUT0 | OUTPUTS1 | ......................... | OUTPUTSk |
| Y | / | / | y1 | ......................... | yk |
| X | bbb.eee | bbb.eee | x1 | ......................... | xk |
where bbb.eee is the control number of the polynomial ( bbb > 008 )
-When the program stops, Rbb thru Ree contain the factorization
-Then, the successive outputs are the different pairs of roots:
xk , yk if flag F00 is clear
or xk+i.yk , xk-i.yk
if F00 is set
However , when deg(p) is odd , the first root is always real: F00 is set but y = 0
-The last pair of roots is indicated by a BEEP, the others by a TONE 9.
Example: Find all the roots of 2.x5+7.x4+20.x3+81.x2+190.x+150
-For instance, 2 STO 11 7 STO 12 20 STO 13 81 STO 14 190 STO 15 150 STO 16 ( control number = 11.016 ) and if we choose u0 = v0 = 0
FIX 8 11.016 XEQ "BRST" the successive uk are displayed and sometime later, we hear a BEEP and X = 11.016
R11 = 2 R12 = 3 // R13 = -2 R14 = 10 // R15 = 4 R16 = 5 whence p(x) = (2x+3).(x2-2x+10).(x2+4x+5) then:
R/S >>> ( TONE 9 )
-1.5 X<>Y 0 with
F00 set the first root is -1.5
( deg(p) is odd , the 1st root is real )
R/S >>> ( TONE 9 )
1 X<>Y 3
with F00 set 1+3.i and 1-3.i
R/S >>> ( BEEP
) -2 X<>Y
1 with F00 set -2+i and
-2-i the
5 roots are -1.5 ; 1+3.i ; 1-3.i ; -2+i ; -2-i
Notes:
-"BRST" always choose u0 = 2 and v0
= PI
-The accuracy is controlled by the display setting: the loop
stops when the rounded increment equals 0
-If you get "DATA ERROR" or if the process seems to diverge, stop the program, change registers R06 & R07 and press XEQ 01
-If you want to see the roots again, press XEQ 10
-If the displayed u-values seem to oscillate between closed numbers
without stopping the loop - it may happen with multiple roots -
Stop the program, set F21 , press R/S. "BRST" will
stop at the next VIEW 06.
Clear X-register an R/S
Cohn's Irreducibility Criterion
-This test concerns a poynomial p(x) = an xn + an-1 xn-1 + ... + a1 x + a0 in Z[X] with ak >= 0 for all k's ( a0 # 0 )
• If there exists an integer m > max { ak }
such that f(m) is a prime, then p(x) is irreducible in Z[X]
Data Registers:
R00 to R02: are used by "PR?"
( Registers Rbb thru Ree are to be initialized before executing
"COHN" )
R03 = bbb.eee R04 = m
• Rbb = an • Rbb+1 = an-1
........................ • Ree = a0 # 0
bbb > 004
Flags: /
Subroutines: "PL"
"PR?"
| STACK | INPUT | OUTPUT |
| X | bbb.eee | 1 |
where bbb.eee = control number of p(x) X-output = 1 if Cohn's criterion is satisfied, bbb > 004
Example1: p(x) = 2 x4 + 3 x2 + 5 x + 1
With [ 2 0 3 5 1 ] = [ R05 R06 R07 R08 R09 ]
5.009 XEQ "COHN" >>>> 1 p(x) is irreducible in Z[X]
R04 = 6 and p(6) = 2731 is a prime
Example2: p(x) = 2 x6 + 3 x5 + 4 x3 + 5 x2 + 2 x + 1
With [ 2 3 0 4 5 2 1 ] = [ R05 R06 R07 R08 R09 R10 R11 ]
5.011 XEQ "COHN" >>>> 1 p(x) is irreducible in Z[X]
R04 = 8 and p(8) = 624977 is a prime
Note:
-If no prime is found, the routine continues until f(m) exceeds 1010
and the HP-41 displays "OUT OF RANGE" and in this case, we cannot conclude.
-If a(x) = an.xn+an-1.xn-1+
... + a1.x+a0 and b(x) = bm.xm+bm-1.xm-1+
... + b1.x+b0 are 2 polynomials,
there are only 2 polynomials q(x) and r(x) such that
a = b.q + r with deg(r) < deg(b)
-"DIV" calculates q and r.
Data Registers: R00 to R03: temp ( Registers Rbb thru Ree & Rbb' thru Ree' are to be initialized before executing "DIV" )
• Rbb = an • Rbb+1 = an-1
, ....... , • Ree = a0
bb > 03
• Rbb' = bm • Rbb'+1 = bm-1 , .......
, • Ree' = b0 bb'
> 03
Flags: /
Subroutines: /
| STACK | INPUTS | OUTPUTS |
| Y | (bbb.eee)a | (bbb.eee)r |
| X | (bbb.eee)b | (bbb.eee)q |
where (bbb.eee)a
= the control number of the dividend
(bbb.eee)r = the control number of the remainder
( all bbb > 003 )
(bbb.eee)b = the control number of the divisor
(bbb.eee)q = the control number of the quotient
Example: a(x) = 2.x5+5.x4-21.x3+23.x2+3.x+5 b(x) = 2.x2-3.x+1 Find q(x) and r(x)
For instance: 2 STO 04
5 STO 05 -21 STO 06 23 STO 07 3 STO
08 5 STO 09 ( control number = 4.009
)
and 2 STO 11
-3 STO 12 1 STO 13
( control number = 11.013 )
4.009 ENTER^
11.013 XEQ "DIV" >>>>
4.007 X<>Y 8.009
( q and r have replaced a ; b is unchanged )
R04 = 1 R05 = 4 R06 = -5 R07
= 2 whence q(x) = x3+4.x2-5.x+2
R08 = 14 R09 = 3
whence r(x) = 14.x + 3
Notes:
-When b(x) is constant, (bbb.eee)r = 1+eee.eee
-Roundoff errors may produce small leading coefficients instead of
zero in the remainder.
"DIV" doesn't work if
deg(a) < deg(b) but in this case, q
= 0 and r = a
-"FCTR" uses "BRST" to factorize in R[X] a polynomial p(x) = an.xn+an-1.xn-1+
... + a1.x+a0 whose coefficients are integers
-The products of their constant coefficients are tested to see if we
get integers after rounding the results.
-Then, p(x) is divided by the possible factors to check if the remainder really equals zero.
-So, the display setting plays a central role.
-To simplify the program, we assume that the leading coefficient
an = +/-1 though it may sometimes work in other cases
too.
Data Registers: R00 to R12+4n: temp ( Registers R01 thru Rnn+1 are to be initialized before executing "FCTR" )
R04 = counter R05 = eee+1 R06 = address of the
first free register
R07 = control number of the list of the irreducible factors in R[X]
R08 = control number of a potential divisor
R09 = control number of f(x)
R10 = counter = 2m-1 - 1 , 2m-1 - 2 , .......
, 0
• R01 = an • R02 = an-1 ........................ • Rnn+1 = a0 and several other registers to store the products of polynomials.
Flags: F00-F01-F26
Subroutines: "BRST" "P2"
"¨PPRD" "MCO"
| STACK | INPUT | OUTPUTk |
| Y | 1.n+1 | / |
| X | m | bbb.eeek |
where 1.n+1 = control number of p(x)
bbb.eeek = control number of the successive factors.
and m is used to round numbers
in FIX m
with n = deg p
Example:
p(x) = x8 - 128 x6 - x5 + 5120 x4
+ 64 x3 - 65538 x2 - 512 x + 262144
With [ R01 R02 R03 R04 R05 R06 R07 R08 R09 ] = [ 1 0 -128 -1 5120 64 -65538 -512 262144 ]
>>> FIX 8 if you want to execute "BRST" in this format, and if you want to execute the rest of "FCTR" in FIX 4
1.009 ENTER^
4 XEQ
"FCTR" >>>> ( TONE 9 ) 60.064
[ R60 R61 R62 R63 R64 ] = [ 1 0 -64 -2 512 ]
R/S >>>> ( BEEP ) 44.048 ( in a few seconds )
[ R44 R45 R46 R47 R48 ] = [ 1 0 -64 1 512 ]
-And p(x) = ( x4 - 64 x2 - 2 x + 512 )( x4 - 64 x2 + x + 512 )
Notes:
-With 1.009 ENTER^ 5
XEQ "FCTR" NO factorization would have been found !
-Before rounding the coefficients, the product was [ 1
0.000000999 -63.99999993 -2.000006500
511.9999998 ]
and -2.000006500 would have been rounded
to -2.00001 which is not an integer!
-If deg f is large and/or if there are many factors in R[X], we'll have
to execute"FCTR" with m = 1
and even then, we could miss some factors because of roundoff
errors...
-So the result is not completely guaranteed. Use "COHN" and/or "IRR?"
to check that the factors are really irreducible...
-If the process seems to diverge while "BRST" is executing, stop the
subroutine, store other guesses in R06 & R07 and press GTO 01
R/S
-But do not press XEQ 01 : that would clear the return
stack !
-If the leading coefficient c = an # +/-1 , "FCTR"
can fail.
-Apply the program to h(x) = cn-1 p(x/c)
-After factorizing h(x), the factorization of p(x) will be easy to
deduce.
A Better Irreducibility Criterion
-Unlike Cohn's criterion, the following one also works if the polynomial
has negative coefficients.
-Let p(x) = an xn + an-1 xn-1
+ ... + a1 x + a0 in Z[X] deg p > 1
( a0 # 0 )
m1 = max { | ai/an
| , i = 0 , 1 , .... , n-1 }
m2 = max { | ai/an
|1/(n-i) , i = 0 , 1 , .... , n-1 }
H = min { 1+m11/(r+1) , 2 m2 } where r is such that an-1 = an-2 = ..... = an-r = 0 , ( r = 0 if an-1 # 0 )
• If there exists an integer m > H such that
p(m) is a prime or +/-1, then p(x) is irreducible in Z[X]
Data Registers:
R00 to R02: temp
( Registers Rbb thru Ree are to be initialized before executing
"IRR?" )
R03 = bbb.eee R04 = m
• Rbb = an • Rbb+1 = an-1
........................ • Ree = a0 # 0
bbb > 004
Flags: /
Subroutines: "PL" "PR?"
| STACK | INPUT | OUTPUT |
| X | bbb.eee | 1 or 0 or E10 |
where bbb.eee = control number of p(x) bbb > 004
if X-output = 1 p is irreducible
if X-output = 0 p is reducible:
we've found an integer root
if X-output = E10 the HP-41 displays
"OUT OF RANGE" and we cannot conclude
Example: f(x) = 2 x7 - 3 x4 + 4 x3 - 9 x2 + 6 x - 1
With [ 2 0 0 -3 4 -9 6 -1 ] = [ R05 R06 R07 R08 R09 R10 R11 R12 ]
5.012 XEQ "IRR?" >>>> 1 thus, f is irreducible
R04 = m = 6 and f(6) = 556559 is a prime
Notes:
-In some cases - when all the coefficients are non-negative and an
is small - Cohn's criterion may be better ( for instance f(x) = x2
+ 1 )
-"P2" solves the equation: a.x2+b.x+c = 0 with a # 0
Data Registers: /
Flags: F00 ( indicates complex roots
)
Subroutines: /
| STACK | INPUTS | OUTPUTS |
| Z | a # 0 | c/a |
| Y | b | y |
| X | c | x |
-If CF 00 the 2 solutions are x , y
-If SF 00 ------------------ x+i.y , x-i.y
Examples: Solve 2.x2 + 3.x - 4 = 0 and 2.x2 + 3.x + 4 = 0
2 ENTER^
3 ENTER^
-4
XEQ "P2" >>>> -2.350781059
X<>Y
0.850781060 Flag F00 is clear,
therefore the 2 solutions are -2.350781059 and
0.850781060
2
ENTER^
3 ENTER^
4 R/S
>>>> -0.75
X<>Y
1.198957881 Flag F00 is set,
therefore the 2 solutions are -0.75 + 1.198957881.i
and -0.75 - 1.198957881.i
-"P3" finds the 3 roots of a.x3+b.x2+c.x+d
by the Cardan's ( or Tartaglia's ) formulae: ( with a # 0 )
-First, the term in x2 is removed by a change of argument,
leading to x3+p.x+q = 0
-Then, x = u+v with u.v = -p/3 leads to a quadratic equation
in u3
Data Registers: /
Flags: F01
( indicates complex roots )
Subroutine: none if d # 0 , "P2" if
d = 0
| STACK | INPUTS | OUTPUTS |
| T | a # 0 | / |
| Z | b | z |
| Y | c | y |
| X | d | x |
-If CF 01 the 3 solutions are x ; y ; z
-If SF 01 the 3 solutions are x ; y+i.z ; y-i.z
Example: Solve 2.x3-5.x2-7.x+1 = 0 and 2.x3-5.x2+7.x+1 = 0
2 ENTER^
-5 ENTER^
-7 ENTER^
1 XEQ "P3"
>>>> 3.467727065 RDN 0.131206041
RDN -1.098933107 which are the 3 real solutions
since flag F01 is clear.
2 ENTER^
-5 ENTER^
7 ENTER^
1
R/S >>>> -0.130131632
RDN 1.315065816 RDN 1.453569820
-Flag F01 is set , therefore the 3 solutions are -0.130131633
; 1.315065816 - 1.453569820.i ; 1.315065816
+ 1.453569820.i
-"P4" solves the equation x4+a.x3+b.x2+c.x+d
= 0 ( if the leading coefficient is not 1 , divide all the
equation by this coefficient ).
-First, the term in x3 is removed by a change of argument,
leading to x4+p.x2+q.x+r = 0 (E')
-Then, the resolvant z3+p.z2/2+(p2-4r).z/16-q2/64
= 0 is solved by "P3" and if we call z1 , z2
, z3 the 3 roots of this equation, the zeros of (E') are
x = z11/2 sign(-q) +/- ( z21/2+z31/2 ) ; x = -(z11/2) sign(-q) +/- ( z21/2-z31/2 )
Data Registers: R00 thru R04: temp
Flags: F01 F02
Subroutine: "P3"
| STACK | INPUTS | OUTPUTS |
| T | a | t |
| Z | b | z |
| Y | c | y |
| X | d | x |
-If CF 01 & CF 02 the 4 solutions are
x ; y ; z ; t
-If SF 01 & CF 02 ------------------
x+i.y ; x-i.y ; z ; t
-If SF 01 & SF 02 ------------------
x+i.y ; x-i.y ; z+i.t ; z-i.t
Example1: Solve x4-2.x3-35.x2+36.x+180 = 0
-2 ENTER^
-35 ENTER^
36 ENTER^
180 XEQ "P4" >>>>
6 RDN 3 RDN -2 RDN
-5 Flags F01 & F02 are clear , the 4 solutions
are 6 ; 3 ; -2 ; -5
Example2: Solve x4-5.x3+11.x2-189.x+522 = 0
-5 ENTER^
11 ENTER^
-189 ENTER^
522
R/S >>>> -2 RDN
5 RDN 3 RDN 6 Flag
F01 is set & F02 is clear , the 4 solutions are -2+5.i ;
-2-5.i ; 3 ; 6
Example3: Solve x4-8.x3+26.x2-168.x+1305 = 0
-8 ENTER^
26 ENTER^
-168 ENTER^
1305 R/S
>>>> -2 RDN 5 RDN
6 RDN 3 Flags F01 & F02
are set , the 4 solutions are -2+5.i ; -2-5.i ; 6+3.i ; 6-3.i
Evaluation of a Polynomial P(x)
-"PL" calculates p(x) = an.xn+an-1.xn-1+
... + a1.x+a0 for any given x-value.
Data Registers: •
Rbb = an • Rbb+1 = an-1 , ....... , •
Ree = a0 (
These n+1 registers are to be initialized before executing "PL" )
Flags: /
Subroutines: /
| STACK | INPUTS | OUTPUTS |
| Y | bbb.eee | / |
| X | x | p(x) |
| L | / | x |
Example: p(x) = 2.x3+3.x2-4.x+7
Find p(5)
-If we store these 4 coefficients into R01 to R04 ( 2 STO 01 3 STO 02 -4 STO 03 7 STO 04 )
1.004 ENTER^
5 XEQ
"PL" >>>> p(5) = 312
-Let a(x) = an.xn+an-1.xn-1+ ... + a1.x+a0 and b(x) = bm.xm+bm-1.xm-1+ ... + b1.x+b0 2 polynomials,
"PPRD" computes and stores the coefficients of p(x) = a(x).b(x)
into contiguous registers. ( you have to specify the first of these registers
)
Data Registers: R00 to R03: temp ( Registers Rbb thru Ree & Rbb' thru Ree' are to be initialized before executing "PPRD" )
• Rbb = an • Rbb+1 = an-1
, ....... , • Ree = a0
bb > 03
• Rbb' = bm • Rbb'+1 = bm-1 , .......
, • Ree' = b0
bb' > 03
Flags: /
Subroutines: /
| STACK | INPUTS | OUTPUTS |
| Z | (bbb.eee)a | / |
| Y | (bbb.eee)b | / |
| X | bbbp | (bbb.eee)p |
where (bbb.eee)a
= the control number of a(x)
and
(bbb.eee)b = the control number of b(x)
(bbb.eee)p = the control number of the product
( all bbb > 003 )
Example: a(x) = 2.x5+5.x4-21.x3+23.x2+3.x+5 b(x) = 2.x2-3.x+1 Find p(x) = a(x).b(x)
For instance: 2 STO 04
5 STO 05 -21 STO 06 23 STO 07 3 STO
08 5 STO 09 ( control number = 4.009
)
2 STO 11 -3 STO 12 1 STO 13
( control number = 11.013 )
and if we want to have p(x) in registers R21 R22 ... etc ...
4.009 ENTER^
11.013 ENTER^
21 XEQ
"PPRD" yields 21.028
-We have R21 = 4 R22 = 4 R23 = -55 R24 = 114 R25 = -84 R26 = 24 R27 = -12 R28 = 5
whence p(x) = 4.x7+4.x6-55.x5+114.x4-84.x3+24.x2-12.x+5
Note: Neither (bbb.eee)a
& (bbb.eee)p nor
(bbb.eee)b & (bbb.eee)p
can overlap.
Primality Test ( Non-negative integers )
-"PR?" takes a non-negative integer n in X-register,
and returns 1 if n is a prime or if n = 1 and 0
otherwise
-The smallest divisor d is in L-register.
Data Registers: R00 = 2
R01 = 4 R02 = 6
Flags: /
Subroutines: /
| STACK | INPUTS | OUTPUTS |
| Y | / | n |
| X | n | 1 or 0 |
| L | / | d |
Examples:
n = 999983
XEQ "PR?" >>>> 1 so
999983 is a prime
n = 9999865009
R/S >>>> 0
so n is not a prime and the smallest non-trivial
divisor = L = 10007
References:
[1] F. Scheid - "Numerical Analysis" - Mc Graw
Hill ISBN 07-055197-9
[2] J-M Ferrard - "Mastering your HP-48" -
D3I Diffusion ISBN 2-908791-04-8
[3] R. Thangadurai - Mathematics Newsletter - Vol
17 #2 - "Irreducibility of Polynomials whose Coefficients are Integers"
[4] Marian Trenkler - "An Algorithm for making Magic
Cubes"
[5] Marian Trenkler - "An Algorithm for Magic Tesseracts"
[6] W. S. Andrews - "Magic Squares and Cubes"
[7] Harm Derksen, Christian Eggermont, Arno van den
Essen - "Multimagic Squares"
[8] W. H. Benson & O. Jacoby - "Magic Cubes,
New Recreations" - Dover - ISBN 0-486-24140-8
[9] http://www.multimagie.com
[10] Pierre Tougne - "Pour la Science" - pages 121 to 127 - Issue#46
- Aout 1981 ( in French )