HP-15C program to compute 3x3 complex matrix determinant



#10

Please refer to the following archived thread:

Complex Matrix Determinant on the 15C

The HP-15C represents an m x n complex-valued matrix as an m x 2n real-valued matrix. The HP-15C's built-in utility functions will transform an
m x m complex-valued square matrix represented in this 'complex' format into an 2m x 2m real-valued square matrix. The determinant of such a matrix is the squared magnitude of the complex-valued determinant of the m x m original matrix.

This procedure can be employed to obtain the magnitude of the determinant of a complex-valued matrix dimensioned up to 4 x 4. However, no built-in capability is available to return the complex-valued determinant itself.

Here's my completed program to calculate the complex-valued determinant of a 3 x 3 complex-valued matrix on the HP-15C. The matrix must be entered as a 3 x 6 real-valued matrix "A" in standard HP-15C 'complex' format. The program was tested with the following matrix A on the HP-15C, and verified using the built-in DET function on the HP-42S:

A =

[ 1+i3  2-i5   7+i1 ] 
[ 4-i2 6+i9 -8+i4 ]
[-3-i7 3+i2 -1+i6 ]

det A = -217 + i231, as computed exactly by this HP-15C program and by the HP-42S. The HP-15C's native capabilities return 100,450.0000 as the squared magnitude of the determinant, whose calculation time is 15 seconds.

If A(1,3) is edited to 7-i1, my HP-15C program returns the exact determinant of -75 + i289, and the HP-15C's native capabilities return 89,146.00002. (The HP-42S computes a slightly-inexact -74.9999999999 + i289, suggesting that the HP-42S does not calculate determinants directly.)

It may be that the exact answers produced by the HP-42S for the original matrix are attributable to the squared magnitude of the determinant exceeding 105. If so, that was pure luck. My only criterion for populating A was that each component of each element be a single-digit integer, with a variety of signs.

This HP-15C program consumes about 25 registers and takes 21 seconds to run. It utilizes indirect access of any matrix {A-E} by storing the descriptor of the matrix into the indirect register I.

Procedure:  Store 3 x 6 matrix in 'complex' form with identifier {A-E}; RCL {A-E}; GSB D.  

  1. Check if dimension of matrix is 3 x 6
  2. (If not, display blinking descriptor and stop)
  3. Compute and store first complex sub-determinant
  4. Compute second sub-determinant
  5. Combine with first sub-determinant and store
  6. Compute third sub-determinant
  7. Combine with first two sub-determinants
  8. Display complex result
  9. NOTE: "(u)" denotes user mode, which is necessary for advancing the matrix-element index.

001 LBL D 056 MATRIX 1 102 MATRIX 1
002 STO I 057 2 103 2
003 RCL DIM I 058 STO 0 104 STO 0
004 f I 059 RCL(i) (u) 105 RCL(i) (u)
005 3 060 RCL(i) 106 RCL(i)
006 ENTER 061 f I 107 f I
007 6 062 3 108 3
008 f I 063 STO 0 109 STO 0
009 - 064 5 110 STO 1
010 x=0? 065 STO 1 111 R_down
011 GTO 8 066 R_down 112 RCL(i) (u)
012 CF 8 067 R_down 113 RCL(i)
013 RCL I 068 RCL(i) (u) 114 f I
014 SF 9 069 RCL(i) 115 *
015 RTN 070 f I 116 2
016 LBL 8 071 * 117 STO 0
017 2 072 2 118 3
018 STO 0 073 STO 0 119 STO 1
019 3 074 5 120 R_down
020 STO 1 075 STO 1 121 R_down
021 RCL(i) (u) 076 R_down 122 RCL(i) (u)
022 RCL(i) (u) 077 R_down 123 RCL(i)
023 f I 078 RCL(i) (u) 124 f I
024 3 079 RCL(i) (u) 125 MATRIX 1
025 STO 0 080 f I 126 3
026 R_down 081 RCL(i) (u) 127 STO 0
027 RCL(i) (u) 082 RCL(i) 128 R_down
028 RCL(i) 083 f I 129 RCL(i) (u)
029 f I 084 * 130 RCL(i)
030 * 085 - 131 f I
031 2 086 MATRIX 1 132 *
032 STO 0 087 3 133 -
033 CLx 088 STO 1 134 MATRIX 1
034 5 089 R_down 135 5
035 STO 1 090 RCL(i) (u) 136 STO 1
036 R_down 091 RCL(i) 137 R_down
037 RCL(i) (u) 092 f I 138 RCL(i) (u)
038 RCL(i) (u) 093 * 139 RCL(i)
039 f I 094 RCL 2 140 f I
040 3 095 RCL 3 141 *
041 STO 1 096 f I 142 RCL 2
042 R_down 097 x<>y 143 RCL 3
043 RCL(i) (u) 098 - 144 f I
044 RCL(i) 099 STO 2 145 +
045 f I 100 Re<>Im 146 PSE
046 * 101 STO 3 147 Re<>Im
047 - 148 PSE
048 MATRIX 1 149 Re<>Im
049 RCL(i) (u) 150 RTN
050 RCL(i)
051 f I
052 *
053 STO 2
054 Re<>Im
055 STO 3


Edited: 19 June 2010, 2:52 p.m. after one or more responses were posted


#11

i wonder if the 42s does its complex matrices through LU decomposition, like the 15c does real ones. This would then give it a working complex matrix math which the 15c lacks.

it might explain the slight inaccuracy in this case.


#12

Quote:
I wonder if the 42s does its complex matrices through LU decomposition, like the 15c does real ones.

It's certainly a possibility, although that is not documented, in my recollection. If so, however, I would assume that the HP-42S internally transforms the complex matrices to the 'quadrupled' real-valued matrices.

The HP-48 and descendants provide an "LU" function that returns the L, U, and P matrices. The Crout algorithm -- instead of the Doolittle algorithm of the HP-15C -- is employed.

-- KS


Edited: 12 June 2010, 11:58 p.m.


#13

Regarding the algorithms used, I tried DET on the matrix

   [ 1+i3  2-i5   7-i1 ] 
A= [ 4-i2 6+i9 -8+i4 ]
[-3-i7 3+i2 -1+i6 ]

I get these results:

HP-42S: -74.9999999999 + i289

Free42: -74.99999999999999999999982 + i288.9999999999999999999998

Free42 uses a straightforward LU decomposition algorithm straight from Numerical Recipes in C. The similarity of the errors makes me suspect the 42S does, too.

Trying this same test case on the 48G, the result is -75 + i289 exactly. This is interesting because the 48G uses the same floating-point representation as the 42S. I wonder: is the 48G more accurate because it uses a different algorithm, or does it use the same algorithm but with 15-digit intermediaries? Given that the 48 series has a lot more RAM than the 42S, it would make sense to use some of that to gain some extra precision.


#14

- The HP-42S uses LU-decomp for complex matrices, not the transformation to real-valued ones. That can be verified: with a matrix in X, executing DET will need to have spare memory for a copy of the matrix (since the original one must be kept in L).
If the transformation to a real-valued matrix would be used, the amount of spare memory needed would be twice the size of the original matrix.
In my case, I have about 5900 bytes free.
With a Cplx 13x13 matrix in level X, that is reduced to about 3200 bytes, leaving enough room for a second 13x13 Complex matrix (a bit over 2700 bytes), but not enough for a 26x26 Real matrix, which is twice the size. And I can run DET without problems.

- The HP48 and successors have a refinement built-in: since the determinant of a matrix with integer elements is itself an integer (or a complex number with integer real and imag components), the calculated result is *rounded* to the nearest integer, which of course occasionally produces a wrong result.

#15

Inspired by Andrew's request I had a look at
Wikipedia
where I found the following formula to calculate the determinant of a 3x3 matrix A:

Calculating the powers of matrix A would be easy but what about the trace function?

|R * *|      |R * * I * *|      |R * * I * * *|      |R R R|          |R|
|* R *| -> |* R * * I *| -> |R * * I * * *| -> |* * *| |*|
|* * R| |* * R * * I| |R * * I 0 0 0| |* * *| |1| |*|
|I * *| |I I I| x |1| = |I| -> R + iI
|* I *| |* * 0| |1| |*|
|* * I| |* * 0| |*|
|* * 0| |*|

With the picture above in my mind I wrote the following program:

                                               X:              Y:              Z:              T:

001 - 42,21,11 LBL A
002 - 45,16,11 RCL MATRIX A A 6 3
003 - 42,16, 4 f MATRIX 4 A 3 6
004 - 45,23,11 RCL DIM A 6 3 A 3 6
005 - 1 1 1 6 3 A 3 6
006 - 40 + 7 3 A 3 6 A 3 6
007 - 42,23,11 f DIM A 7 3 A 3 7 A 3 7
008 - 34 x<>y 3 7 A 3 7 A 3 7
009 - 43 36 g LSTx 1 3 7 A 3 7
010 - 42,23,12 f DIM B 1 3 7 A 3 7
011 - 44,16,12 STO MATRIX B 1 3 7 A 3 7
012 - 43 33 g R-^ A 3 7 1 3 7
013 - 45,16,12 RCL MATRIX B b 3 1 A 3 7 1 3
014 - 42,26,13 f RESULT C b 3 1 A 3 7 1 3
015 - 42,16, 5 f MATRIX 5 C 7 1 1 3 3
016 - 42,16, 1 f MATRIX 1 C 7 1 1 3 3
017 - 45 13 RCL C R C 7 1 1 3
018 - 43 33 g R-^ 3 R C 7 1 1
019 - 44,40, 0 STO + 0 3 R C 7 1 1
020 - 34 x<>y R 3 C 7 1 1
021 - 45 13 RCL C I R 3 C 7 1
022 - 42 25 f I R + iI 3 C 7 1 C 7 1
023 - 43 32 RTN R + iI 3 C 7 1 C 7 1

But unfortunately there's not enough memory to do the calculation of the powers of a complex 3x3 Matrix.

Still the program might be useful.

Best regards

Thomas


Edited: 15 June 2010, 2:18 a.m.


#16

Thomas --

Thank you for the contribution. The HP-15C, due to modest data-storage capability and computational speed, is certainly not the ideal tool for heavy-duty matrix calculations. However, it's good to demonstrate that certain things can be done within the limited resources available.

-- KS

Edited: 17 June 2010, 1:54 a.m.

#17

Karl,

Thanks very much for all the work you have done on this.

I originally thought it might be possible to achieve this without having to resort to 'brute force', e.g. using the trace formula, but it appears that the 15C just doesn't have the capacity to do so.

However, I do like your original solution: get a 42S!

Thanks again,

Andrew


#18

Andrew H --

Quote:
Thanks very much for all the work you have done on this.

You're welcome. It's a topic that came up a few years ago, but which I never tackled.

My program is now improved, in that any matrix {A-E} can be used, not just A. The HP-15C allows matrix descriptors to be stored in any register; the "I" register allows indirect access to an element in any matrix. Note also that the program checks matrix dimensions in a concise manner for correctness by using built-in logical comparisons in complex mode. The thoroughness and attention to detail exhibited by the HP-15C's developers has always impressed me -- it is a case study of first-rate software engineering.

Quote:
I originally thought it might be possible to achieve this without having to resort to 'brute force', e.g. using the trace formula, but it appears that the 15C just doesn't have the capacity to do so.

'Brute force' is a bit pejorative; I'd prefer to call it the 'direct' or 'straightforward' method. There's a certain amount of computational grinding involved with multiplying two complex-valued matrices, as well -- the traces are the easy part.

In fact, my modest program could be utilized four times in succession to calculate the complex-valued determinant of a 4 x 4 matrix:

  1. Store three columns of the last three rows of a 4 x 4 complex matrix in two parts: a 3 x 6 real matrix, and the remaining 'short column' as a separate 3 x 2 matrix.
  2. Use the program to calculate one sub-determinant, multiply by the corresponding complex element of the top row, and accumulate in registers R4 and R5.
  3. Swap one 'short column' with the separate column (maybe write a short utility program for that), and repeat step 2.
  4. Repeat step 3 twice more.

This procedure will fit within the HP-15C's capabilities: 18 registers for the 3 x 3 complex submatrix, 6 more for the 3 x 1 complex submatrix, 25 registers for the program, 4 registers for temporary storage, 5 registers for complex mode, with space left over for a column-swapping utility program.

Or, use an HP-42S or RPL-based calc...

-- KS


Edited: 20 June 2010, 10:32 p.m.


Possibly Related Threads...
Thread Author Replies Views Last Post
  HP Prime: complex numbers in CAS. Alberto Candel 1 396 12-06-2013, 02:36 PM
Last Post: parisse
  AFTER HP-Prime update, Shift+Matrix CRASHES Joseph Ec 3 509 12-06-2013, 11:06 AM
Last Post: Joseph Ec
  [HP Prime] Plots containing complex numbers bug? Chris Pem10 7 745 12-05-2013, 07:40 AM
Last Post: cyrille de Brébisson
  HP Prime Matrix TERRIBLE bug and question uklo 19 1,432 11-25-2013, 12:10 PM
Last Post: Mic
  HP Prime: editing a matrix Alberto Candel 6 540 11-20-2013, 06:26 PM
Last Post: Helge Gabert
  Complex Number Entry on Prime Jeff O. 19 1,522 11-16-2013, 12:34 PM
Last Post: Jeff O.
  Absolute Value and Matrix BruceTTT 5 549 11-11-2013, 11:52 PM
Last Post: Walter B
  HP Prime: run a program in another program Davi Ribeiro de Oliveira 6 607 11-11-2013, 08:28 PM
Last Post: Davi Ribeiro de Oliveira
  HP Prime: Rounding error in determinant Stephan Matthys 3 402 10-25-2013, 09:29 PM
Last Post: Walter B
  HP Prime complex results Javier Goizueta 0 248 10-06-2013, 12:59 PM
Last Post: Javier Goizueta

Forum Jump: