+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: Not HP Calculators (/forum-7.html)
+--- Forum: Not quite HP Calculators - but related (/forum-8.html)
+--- Thread: Inverse of 3x3 Matrix with Complex Numbers (/thread-13286.html)
Inverse of 3x3 Matrix with Complex Numbers - deetee - 07-15-2019 11:59 AM
Hello!
My next calculator should be able to work with 3x3 matrices with complex numbers (determinant, invert, add, substract, multiply, ...). That for I want to impement fast explicit formulae (found on Wikipedia) rather than some algorithms.
The calculator itself works with two double stacks for real and imaginary parts of numbers (stack[4] and stacki[4]). Basic stack operations are done with _add(), _sub() and _mult().
For matrix operations I use matrix A (ma[3][3] and mai[3][3]) and B (mb[3][3] and mbi[3][3]) and a temporary matrix C (mc[3][3] and mci[3][3]).
So far I implemented the following code which works fine:
Code:
static void _mdet(void) { // MATRIX DET A
mdet = mdeti = 0.0;
stack[2] = ma[0][0]; stacki[2] = mai[0][0];
stack[1] = ma[1][1]; stacki[1] = mai[1][1];
stack[0] = ma[2][2]; stacki[0] = mai[2][2];
_mult(); _mult(); mdet += stack[0]; mdeti += stacki[0];
stack[2] = ma[0][1]; stacki[2] = mai[0][1];
stack[1] = ma[1][2]; stacki[1] = mai[1][2];
stack[0] = ma[2][0]; stacki[0] = mai[2][0];
_mult(); _mult(); mdet += stack[0]; mdeti += stacki[0];
stack[2] = ma[0][2]; stacki[2] = mai[0][2];
stack[1] = ma[1][0]; stacki[1] = mai[1][0];
stack[0] = ma[2][1]; stacki[0] = mai[2][1];
_mult(); _mult(); mdet += stack[0]; mdeti += stacki[0];
stack[2] = ma[2][0]; stacki[2] = mai[2][0];
stack[1] = ma[1][1]; stacki[1] = mai[1][1];
stack[0] = ma[0][2]; stacki[0] = mai[0][2];
_mult(); _mult(); mdet -= stack[0]; mdeti -= stacki[0];
stack[2] = ma[2][1]; stacki[2] = mai[2][1];
stack[1] = ma[1][2]; stacki[1] = mai[1][2];
stack[0] = ma[0][0]; stacki[0] = mai[0][0];
_mult(); _mult(); mdet -= stack[0]; mdeti -= stacki[0];
stack[2] = ma[2][2]; stacki[2] = mai[2][2];
stack[1] = ma[1][0]; stacki[1] = mai[1][0];
stack[0] = ma[0][1]; stacki[0] = mai[0][1];
_mult(); _mult(); mdet -= stack[0]; mdeti -= stacki[0];
stack[0] = mdet; stacki[0] = mdeti;
}
static void _minv(void) { // MATRIX INVERSE A
_mdet();
stack[3] = ma[1][1]; stacki[3] = mai[1][1]; // C00
stack[2] = ma[2][2]; stacki[2] = mai[2][2];
stack[1] = ma[2][1]; stacki[2] = mai[2][1];
stack[0] = ma[1][2]; stacki[3] = mai[1][2];
_mult(); _rot(); _mult(); _rotup(); _sub();
mc[0][0] = stack[0]; mci[0][0] = stacki[0];
stack[3] = ma[0][2]; stacki[3] = mai[0][2]; // C01
stack[2] = ma[2][1]; stacki[2] = mai[2][1];
stack[1] = ma[0][1]; stacki[2] = mai[0][1];
stack[0] = ma[2][2]; stacki[3] = mai[2][2];
_mult(); _rot(); _mult(); _rotup(); _sub();
mc[0][1] = stack[0]; mci[0][1] = stacki[0];
stack[3] = ma[0][1]; stacki[3] = mai[0][1]; // C02
stack[2] = ma[1][2]; stacki[2] = mai[1][2];
stack[1] = ma[0][2]; stacki[2] = mai[0][1];
stack[0] = ma[1][1]; stacki[3] = mai[1][2];
_mult(); _rot(); _mult(); _rotup(); _sub();
mc[0][2] = stack[0]; mci[0][2] = stacki[0];
stack[3] = ma[1][2]; stacki[3] = mai[1][2]; // C10
stack[2] = ma[2][0]; stacki[2] = mai[2][0];
stack[1] = ma[1][0]; stacki[2] = mai[1][2];
stack[0] = ma[2][2]; stacki[3] = mai[2][0];
_mult(); _rot(); _mult(); _rotup(); _sub();
mc[1][0] = stack[0]; mci[1][0] = stacki[0];
stack[3] = ma[0][0]; stacki[3] = mai[0][0]; // C11
stack[2] = ma[2][2]; stacki[2] = mai[2][2];
stack[1] = ma[0][2]; stacki[2] = mai[0][0];
stack[0] = ma[2][0]; stacki[3] = mai[2][2];
_mult(); _rot(); _mult(); _rotup(); _sub();
mc[1][1] = stack[0]; mci[1][1] = stacki[0];
stack[3] = ma[0][2]; stacki[3] = mai[0][2]; // C12
stack[2] = ma[1][0]; stacki[2] = mai[1][0];
stack[1] = ma[0][0]; stacki[2] = mai[0][0];
stack[0] = ma[1][2]; stacki[3] = mai[1][2];
_mult(); _rot(); _mult(); _rotup(); _sub();
mc[1][2] = stack[0]; mci[1][2] = stacki[0];
stack[3] = ma[1][0]; stacki[3] = mai[1][0]; // C20
stack[2] = ma[2][1]; stacki[2] = mai[2][1];
stack[1] = ma[1][1]; stacki[2] = mai[1][1];
stack[0] = ma[2][0]; stacki[3] = mai[2][0];
_mult(); _rot(); _mult(); _rotup(); _sub();
mc[2][0] = stack[0]; mci[2][0] = stacki[0];
stack[3] = ma[0][1]; stacki[3] = mai[0][1]; // C21
stack[2] = ma[2][0]; stacki[2] = mai[2][0];
stack[1] = ma[0][0]; stacki[2] = mai[0][0];
stack[0] = ma[2][1]; stacki[3] = mai[2][1];
_mult(); _rot(); _mult(); _rotup(); _sub();
mc[2][1] = stack[0]; mci[2][1] = stacki[0];
stack[3] = ma[0][0]; stacki[3] = mai[0][0]; // C22
stack[2] = ma[1][1]; stacki[2] = mai[1][1];
stack[1] = ma[0][1]; stacki[2] = mai[0][1];
stack[0] = ma[1][0]; stacki[3] = mai[1][0];
_mult(); _rot(); _mult(); _rotup(); _sub();
mc[2][2] = stack[0]; mci[2][2] = stacki[0];
for (byte j = 0; j < 3; j++) for (byte i = 0; i < 3; i++) { // Aij = 1/det * Cij
stack[0] = mdet; stacki[0] = mdeti;
stack[1] = mc[j][i]; stacki[1] = mci[j][i];
_div();
ma[j][i] = stack[0]; mai[j][i] = stacki[0];
}
}Unfortunately this code "eats up" (too) many valuable bytes. I'm sure it can be optimized and shrinked - at least with some loops. But I must resign to manage such loop indices.
Does anyone of you (very smart guys) see how to automate this code in an efficient way?
Any idea and hint is highly welcome.
Regards deetee
RE: Inverse of 3x3 Matrix with Complex Numbers - deetee - 07-16-2019 05:47 AM
Hi!
At least I found some improvement by myself to reduce the code size by summarizing similar code in subroutines - and call them with (a lot of) index parameters:
Code:
static void _mdetrc(byte r1, byte r2, byte r3, byte c1, byte c2, byte c3, int8_t mult) {
stack[2] = ma[r1][c1]; stacki[2] = mai[r1][c1];
stack[1] = ma[r2][c2]; stacki[1] = mai[r2][c2];
stack[0] = ma[r3][c3]; stacki[0] = mai[r3][c3];
_mult(); _mult();
mdet += mult * stack[0]; mdeti += mult * stacki[0];
}
static void _mdet(void) { // MATRIX DET A
mdet = mdeti = 0.0;
_mdetrc(0, 1, 2, 0, 1, 2, 1);
_mdetrc(0, 1, 2, 1, 2, 0, 1);
_mdetrc(0, 1, 2, 2, 0, 1, 1);
_mdetrc(2, 1, 0, 0, 1, 2, -1);
_mdetrc(2, 1, 0, 1, 2, 0, -1);
_mdetrc(2, 1, 0, 2, 0, 1, -1);
stack[0] = mdet; stacki[0] = mdeti;
}
static void _minvrc(byte r1, byte r2, byte r3, byte r4, byte c1, byte c2, byte c3, byte c4, byte r, byte c) {
stack[3] = ma[r1][c1]; stacki[3] = mai[r1][c1]; // Calculates inverse matrix element
stack[2] = ma[r2][c2]; stacki[2] = mai[r2][c2];
stack[1] = ma[r3][c3]; stacki[1] = mai[r3][c3];
stack[0] = ma[r4][c4]; stacki[0] = mai[r4][c4];
_mult(); _rot(); _mult(); _rotup(); _sub(); _push();
stack[0] = mdet; stacki[0] = mdeti;
_div();
mc[r][c] = stack[0]; mci[r][c] = stacki[0];
}
static void _minv(void) { // MATRIX INVERSE A
_mdet();
_minvrc(1, 2, 2, 1, 1, 2, 1, 2, 0, 0); // C00
_minvrc(0, 2, 0, 2, 2, 1, 1, 2, 0, 1); // C01
_minvrc(0, 1, 0, 1, 1, 2, 2, 1, 0, 2); // C02
_minvrc(1, 2, 1, 2, 2, 0, 0, 2, 1, 0); // C10
_minvrc(0, 2, 0, 2, 0, 2, 2, 0, 1, 1); // C11
_minvrc(0, 1, 0, 1, 2, 0, 0, 2, 1, 2); // C12
_minvrc(1, 2, 1, 2, 0, 1, 1, 0, 2, 0); // C20
_minvrc(0, 2, 0, 2, 1, 0, 0, 1, 2, 1); // C21
_minvrc(0, 1, 0, 1, 0, 1, 1, 0, 2, 2); // C22
memcpy(ma, mc, 3 * 3 * sizeof(double)); // A = C
memcpy(mai, mci, 3 * 3 * sizeof(double));
}Thanks in advance for any hint.
Regards deetee
RE: Inverse of 3x3 Matrix with Complex Numbers - toml_12953 - 07-16-2019 10:34 AM
(07-16-2019 05:47 AM)deetee Wrote: Hi!
At least I found some improvement by myself to reduce the code size by summarizing similar code in subroutines - and call them with (a lot of) index parameters:
Code:
static void _mdetrc(byte r1, byte r2, byte r3, byte c1, byte c2, byte c3, int8_t mult) {
stack[2] = ma[r1][c1]; stacki[2] = mai[r1][c1];
stack[1] = ma[r2][c2]; stacki[1] = mai[r2][c2];
stack[0] = ma[r3][c3]; stacki[0] = mai[r3][c3];
_mult(); _mult();
mdet += mult * stack[0]; mdeti += mult * stacki[0];
}
static void _mdet(void) { // MATRIX DET A
mdet = mdeti = 0.0;
_mdetrc(0, 1, 2, 0, 1, 2, 1);
_mdetrc(0, 1, 2, 1, 2, 0, 1);
_mdetrc(0, 1, 2, 2, 0, 1, 1);
_mdetrc(2, 1, 0, 0, 1, 2, -1);
_mdetrc(2, 1, 0, 1, 2, 0, -1);
_mdetrc(2, 1, 0, 2, 0, 1, -1);
stack[0] = mdet; stacki[0] = mdeti;
}
static void _minvrc(byte r1, byte r2, byte r3, byte r4, byte c1, byte c2, byte c3, byte c4, byte r, byte c) {
stack[3] = ma[r1][c1]; stacki[3] = mai[r1][c1]; // Calculates inverse matrix element
stack[2] = ma[r2][c2]; stacki[2] = mai[r2][c2];
stack[1] = ma[r3][c3]; stacki[1] = mai[r3][c3];
stack[0] = ma[r4][c4]; stacki[0] = mai[r4][c4];
_mult(); _rot(); _mult(); _rotup(); _sub(); _push();
stack[0] = mdet; stacki[0] = mdeti;
_div();
mc[r][c] = stack[0]; mci[r][c] = stacki[0];
}
static void _minv(void) { // MATRIX INVERSE A
_mdet();
_minvrc(1, 2, 2, 1, 1, 2, 1, 2, 0, 0); // C00
_minvrc(0, 2, 0, 2, 2, 1, 1, 2, 0, 1); // C01
_minvrc(0, 1, 0, 1, 1, 2, 2, 1, 0, 2); // C02
_minvrc(1, 2, 1, 2, 2, 0, 0, 2, 1, 0); // C10
_minvrc(0, 2, 0, 2, 0, 2, 2, 0, 1, 1); // C11
_minvrc(0, 1, 0, 1, 2, 0, 0, 2, 1, 2); // C12
_minvrc(1, 2, 1, 2, 0, 1, 1, 0, 2, 0); // C20
_minvrc(0, 2, 0, 2, 1, 0, 0, 1, 2, 1); // C21
_minvrc(0, 1, 0, 1, 0, 1, 1, 0, 2, 2); // C22
memcpy(ma, mc, 3 * 3 * sizeof(double)); // A = C
memcpy(mai, mci, 3 * 3 * sizeof(double));
}
Thanks in advance for any hint.
Regards deetee
Have you benchmarked your routine vs. the built-in routines? I'd be curious to see if there's a significant difference one way or the other.
RE: Inverse of 3x3 Matrix with Complex Numbers - deetee - 07-16-2019 12:07 PM
Hi Tom!
No I didn't benchmark the explicit method and an algorithm based method, but I think the explicit method is much faster. On my Arduino (ATMEGA328P) the result was shown immediately.
On the other hand algorithms are more flexible regarding the (matrix) dimensions.
My intention was to solve matrices with complex numbers up to 3x3. The input of 9 complex numbers is quite painful. For bigger dimensions one should rather use a PC.
Regards
deetee
PS: BTW I'm not using an existing calculator. I try to build my own hardware.