+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Software Libraries (/forum-10.html)
+--- Forum: HP-41C Software Library (/forum-11.html)
+--- Thread: HP-41 Complex Determinants w/ SandMatrix and 41Z (/thread-10897.html)
HP-41 Complex Determinants w/ SandMatrix and 41Z - Ángel Martin - 06-11-2018 08:51 AM
Complex Determinants with the SandMatrix and 41Z Modules
The programs below are a first-pass successful attempt at calculating Complex Matrix determinants up to order 4.
The Complex Matrix is to be stored using the SandMatrix convention - which is identical to the HP-41 Advantage's. With this convention each complex number is represented by four elements in the complex matrix - refer to the manuals for details.
Code:
[[ Re(aij) -Im(aij) ]
[ Im(aij) Re(aij) ]]The SandMatrix comes well-equipped with routines to calculate the TRACE and integer powers of a matrix (M^2 and MPWR), therefore it lends itself rather nicely to the direct formulas using those elements, as described at:
https://en.wikipedia.org/wiki/Determinant
The program includes 4 sections:
1. A Complex Matrix trace routine, "CMTR"
2. Determinant of a Complex 2x2 Matrix, "CMDT2"
3. Determinant of a Complex 3x3 Matrix, "CMDT3"
4. Determinant of a Complex 4x4 Matrix, "CMDT4"
The last three can be grouped into a single one, using the matrix DIMension as criteria to branch to the appropriate section - but I left it as-is for easier readability. Also the code can be optimized putting repeated lines in independent subroutines; yet ditto as before.
Only CMDT4 is commented, as it's a superset of the other two and can be used as a reference. Here's the program listing:
Code:
01 LBL "CMTR" 63 2
02 XEQ 01 64 /
03 ZAVIEW 65 MNAME?
04 RTN 66 RTN
05 LBL "CMDT2" 67 LBL "CMDT4"
06 ASTO 02 68 ASTO 02
07 XEQ 01 69 XEQ 01 tr(A) = 3+i
08 Z^2 70 ZRPL^ fills Z-stack
09 ZENTER^ 71 "|-,#"
10 "|-,#" 72 MAT= copy to scratch
11 MAT= 73 "#" scratch matrix
12 "#" 74 3
13 M^2 75 MPWR cube power
14 XEQ 01 76 XEQ 01 tr(A^3) = -75+50i
15 Z- 77 Z* tr(A).tr(A^3)
16 2 78 8
17 GTO 02 79 ST* Z
18 LBL "CMDT3" 80 * 8.tr(A).tr(A^3)
19 ASTO 02 81 Z<>W tr(A)
20 XEQ 01 82 Z^2 tr(A)^2
21 Z^3 83 6
22 LASTZ 84 ST* Z
23 ZENTER^ 85 * 6. tr(A)^2
24 "|-,#" 86 ZENTER^
25 MAT= 87 CLA
26 "#" 88 ARCL 02
27 M^2 89 "|-,#"
28 XEQ 01 90 MAT=
29 Z* 91 "#" scratch matrix
30 3 92 M^2 squared
31 ST* Z 93 XEQ 01 tr(A^2) = 2 -4i
32 * 94 Z* 6. tr(A)^2 . tr(A^2)
33 Z- 95 LASTZ TR(A^2)
34 ZENTER^ 96 ZRDN puts it in level "V"
35 CLA 97 Z- 8.tr(A).tr(A^3) - 6. tr(A)^2 . tr(A^2)
36 ARCL 02 98 Z<>W tr(A)
37 "|-,#" 99 Z^2
38 MAT= 100 Z^2 tr(A)^4
39 "#" 101 Z+ tr(A)^4 - 6. tr(A)^2 . tr(A^2) + 8.tr(A).tr(A^3)
40 3 102 ZRUP tr(A^2)
41 MPWR 103 Z^2 tr(A^2)^2
42 XEQ 01 104 3
43 2 105 ST* Z
44 ST* Z 106 *
45 * 107 Z+ 3.tr(A^2)^2 + tr(A)^4 - 6. tr(A)^2 . tr(A^2) + 8.tr(A).tr(A^3)
46 Z+ 108 ZENTER^
47 6 109 M^2 fourth power
48 GTO 02 110 XEQ 01 tr(A^4) = -160(1-i)
49 LBL 01 111 6
50 , 112 ST* Z
51 MSIJA 113 * 6.tr(A^4)
52 I+ 114 Z- 3.tr(A^2)^2 + tr(A)^4 - 6. tr(A)^2 . tr(A^2) + 8.tr(A).tr(A^3) - 6.r(A^4)
53 MRC+ 115 24
54 LBL 00 116 LBL 02
55 I+ 117 ST/ Z
56 J+ 118 /
57 J+ 119 ZAVIEW
58 MRC+ 120 "#"
59 + 121 PURFL
60 FC? 09 122 CLA
61 GTO 00 123 ARCL 02
62 MTRACE 124 ENDThe complex matrix won't be altered in any way, as all operations are made on a scratch copy. It can be stored in X-Mem, CL_Y-Mem, or standard Data registers area. The easiest way to enter the matrix is by using the CMEDIT routine - which expects the matrix name in ALPHA. It expects the matrix already created, using 2n x 2n as dimension - with "n" being the order.
If you place it in the standard registers area, be aware that data registers R00, R01 are used by the routine MPWR for scratch. Additionally, data register R02 is used to store the Matrix Name (thus can't exceed 6 characters).
As you can see there are numerous 41Z functions - used for the complex arithmetic using the Complex Stack. This has the additional advantage that doesn't require additional data registers, be that standard or CL Y-RAM.
Examples.-
Calculate the determinant of the 4x4 Complex Matrix:
Code:
[[ 1+i 2+2i 3+3i 4+4i ]
[ 0 1 -3-3i -4-4i ]
[ -1+i 1-i 1 i ]
[ -i -1+1 1 0 ]]Solution: det = -62-8i
https://www.wolframalpha.com/input/?i=det(%7B%7B1%2Bi,+2%2B2i,+3%2B3i,+4%2B4i%7D,+%7B0,+1,+-3-3i,+-4-4i%7D,+%7B-1%2Bi,+1-i,+1,+i%7D,+%7B-i,+-1%2Bi,+1,+0%7D%7D)
The program is slow in non-turbo settings- there are lots of moving pars behind the scenes, despite the straight-forward program listing.
Using TURBO_50 the 4x4 determinant is obtained in 5 seconds approx.
The accuracy for integer matrices holds up nicely, giving exact integer real and imaginary parts in the solution.
Feel free to add your comments and suggestions for improvement.
Cheers,
ÁM
RE: HP-41 Complex Determinants w/ SandMatrix and 41Z - Valentin Albillo - 06-11-2018 11:09 PM
.
Hi, Ángel:
As I told you in my previous post on this subject (which was in another thread), for low dimensionality (say up to 5x5 or 6x6) I prefer the simple, general expansion by minors procedure over complicated, specific formulae for each dimension, which entail traces (cheap) and powers of matrices (expensive).
The expansion is very simple to code and works for any dimensions from 2x2 upwards though it's only reasonably efficient for low orders because the computation time grows like n!. This is my non-optimized, recursive implementation as a 4-line subprogram in HP-71B code valid for both real and complex matrices of any size:
8 SUB MDET(A(,),D) @ N=UBND(A,1) @ IF N=2 THEN D=A(1,1)*A(2,2)-A(1,2)*A(2,1) @ END
9 D=0 @ IF TYPE(A)<6 THEN REAL E,B(N-1,N-1) ELSE COMPLEX E,B(N-1,N-1)
10 FOR K=1 TO N @ FOR I=2 TO N @ C=1 @ FOR J=1 TO N @ IF J#K THEN B(I-1,C)=A(I,J) @ C=C+1
11 NEXT J @ NEXT I @ CALL MDET(B,E) @ D=D-(-1)^K*A(1,K)*E @ NEXT K
To use it, you simply call it from the keyboard or your own program, passing as parameters (both by reference) the matrix itself (real or complex) and a variable where the determinant will be returned (same type real or complex as the matrix), i.e.:
CALL MDET(M,D)
Upon return from the subprogram, D holds the value of the determinant. The subprogram finds out the type (real or complex) of the matrix and its dimensions so there's no need to pass this information as additional parameters.
The code could be greatly optimized. For instance, it could attempt to find the row or column which has the most zero values and do the minor expansion by that row or column, which would save significant time, and it could also implement the exact formula for 3x3 determinants as well for extra speed, but this is intended just for proof-of-concept, not production.
One way to check the MDET subprogram is by running this sample driver program:
1 DESTROY ALL @ OPTION BASE 1
2 DATA 58,71,67,36,35,19,60,50,71,71,56,45,20,52,64,40,84,50,51,43,69,31,28
3 DATA 41,54,31,18,33,45,23,46,38,50,43,50,41,10,28,17,33,41,46,66,72,71,38
4 DATA 40,27,69,(1,2),(2,3),(3,1),(-1,2),(2,-1),(-1,-1),(0,3),(-2,0),(2,2)
5 REAL A(7,7),D @ READ A @ DISP DET(A) @ DISP @ GOSUB 7
6 COMPLEX A(3,3),D @ READ A @ DISP @ GOSUB 7 @ END
7 MAT DISP A; @ DISP @ CALL MDET(A,D) @ DISP "D =";D @ RETURN
which defines both a 7x7 real matrix and a 3x3 complex one and upon running displays them and calls MDET to compute and subsequently display their determinants, like this:
>RUN
.97095056196
58 71 67 36 35 19 60
50 71 71 56 45 20 52
64 40 84 50 51 43 69
31 28 41 54 31 18 33
45 23 46 38 50 43 50
41 10 28 17 33 41 46
66 72 71 38 40 27 69
D = 1
(1,2) (2,3) (3,1)
(-1,2) (2,-1) (-1,-1)
(0,3) (-2,0) (2,2)
D = (44,-6)
The very first number displayed, .97095056196, is the result returned by the Math ROM function DET when applied to the 7x7 real matrix (DET doesn't work for complex matrices). As you can see from the output produced by MDET, the actual, correct value for the 7x7 determinant is exactly 1, not .97095056196, so this is another important advantage of MDET: it produces exact integer values for integer matrices (as long as no intermediate products or sums exceed 12 digits), which DET might not.
The 7x7 particular matrix used as a test is my 7x7 integer-element "Albillo Matrix #1", as discussed in my article "Mean Matrices". Its determinant is exactly 1, but the DET function in the Matrix ROM evaluates it as .97095056196 despite its small size, simple 2-digit integer elements and 15-digit internal accuracy. See the article for more details on this and other hugely more difficult Albillo Matrices, some of which can put to shame the WP-34S famous 34-digit accuracy. In particular perhaps you would consider:
- implementing expansion by minors for complex matrices' determinants in 41C code using your awesome microcoded complex and matrix functions, to see how it fares in terms of speed versus your current 4x4 implementation.
- check your matrix functions or program using my 7x7 matrix as a test case, to see if they return the correct determinant (1). Even though the matrix is real it can be used with your algorithms expecting a complex matrix, right ?
For dimensions over 6x6 or 7x7, I'd simply perform the classical LU decomposition for speed but the accuracy for integer matrices would certainly not be as good.
Regards.
V.
.
RE: HP-41 Complex Determinants w/ SandMatrix and 41Z - Ángel Martin - 06-12-2018 04:57 AM
Hola Valentín,
Thanks for adding your comments and providing such a detailed explanation. When I see those deceivingly simple four-liners of your code I really feel machine envy!!
At the end you're going to make me switch to the 71B ! :-D
I need some time to absorb your suggestions; what you're pointing out of course makes all the sense in the world and I too would prefer a general solution.
Saludos,
Á.
RE: HP-41 Complex Determinants w/ SandMatrix and 41Z - Valentin Albillo - 06-12-2018 10:41 PM
.
Hola, Ángel:
(06-12-2018 04:57 AM)Ángel Martin Wrote: Thanks for adding your comments and providing such a detailed explanation. When I see those deceivingly simple four-liners of your code I really feel machine envy!!
You're welcome. As for the 4-liner, if HP had done their work and included subarray handling in the 71B Math ROM (as featured in the excellent Series 80 Matrix ROM) we would be talking about a 3-liner (and a much faster one, two nested FOR loops rendered unnecessary).
But alas, they didn't, they'd rather fill up 5 Kb of precious ROM space with zeros than provide that simple and most useful functionality. Bunch of <utterly unpublishable expletive> !!
Quote:At the end you're going to make me switch to the 71B ! :-D
One can only hope. Perhaps some day you'll finally see the light, I pray you do ... :-D
Quote:I need some time to absorb your suggestions; what you're pointing out of course makes all the sense in the world and I too would prefer a general solution.
Certainly. Low-dimensionality is borderly acceptable for lesser machines with limited RAM but if I'm not mistaken the whole point of your latest work with the Advantage ROM was to be able to access up to 55\(^2\) = 3,025 data registers, thus to then go and produce a 4x4 program is kinda lame ... Just kidding ! ;-)
NxN or bust !
Nasnoches.
V.
.
RE: HP-41 Complex Determinants w/ SandMatrix and 41Z - J-F Garnier - 06-13-2018 07:19 AM
(06-12-2018 04:57 AM)Ángel Martin Wrote: Thanks for adding your comments and providing such a detailed explanation. When I see those deceivingly simple four-liners of your code I really feel machine envy!!
At the end you're going to make me switch to the 71B ! :-D
Interesting how the 41 and 71 numbers appear in the "Albillo Matrix #1" :
58 71 67 36 35 19 60
50 71 71 56 45 20 52
64 40 84 50 51 43 69
31 28 41 54 31 18 33
45 23 46 38 50 43 50
41 10 28 17 33 41 46
66 72 71 38 40 27 69
No 75 however ... :-)
J-F
RE: HP-41 Complex Determinants w/ SandMatrix and 41Z - Valentin Albillo - 06-13-2018 10:15 AM
(06-13-2018 07:19 AM)J-F Garnier Wrote: Interesting how the 41 and 71 numbers appear in the "Albillo Matrix #1" :
58 71 67 36 35 19 60
50 71 71 56 45 20 52
64 40 84 50 51 43 69
31 28 41 54 31 18 33
45 23 46 38 50 43 50
41 10 28 17 33 41 46
66 72 71 38 40 27 69
That was intentional, there's also a 67 and a 35 ... :-D
Quote:No 75 however ... :-)
Right, once the vastly superior 71 is there there's simply no place for the clunker. ;-)
Best regards.
V.
.