HP-41 Complex Determinants w/ SandMatrix and 41Z

[Ángel Martin]

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  END

The 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=de...,+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

[Valentin Albillo]

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

[Ángel Martin]

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,
Á.

[Valentin Albillo]

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

[J-F Garnier]

(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

[Valentin Albillo]

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