Matrix inversion for 4x4

[frenchjam]

Hi,

Someone referred me to your library, which is a big help to me. I don't know if you are maintaining this library or if you intended it for public consumption, but I do have one question / issue. Have you tested it with 4x4 or higher order matrices? I seem to be getting anomalous results when I try to invert a 4x4 identity matrix. Code is pretty simple:

matrix A = identity_matrix(4);
matrix B = A.invert();
cout << B << endl

Output is:

[ 1, 0, 0, 0; 0, -1, 0, 0; 0, 0, 1, 0; 0, 0, 0, -1 ]

Perhaps I am doing something stupid. Any insights for me?

I am using Visual C++ 2017 Community.

[frenchjam]

By the way, I had to modify det() to get it to compile under Visual C++. It is entirely possible that I introduced a bug in doing so. I am going to try and track it down. But if you could tell me if the above works under gcc, for instance, I can be more confident that I am the one who introduced the problem in my code.

[frenchjam]

The problem appears to be for square matrices with even numbers of rows and columns, except 2x2. To see the effect:

for (i = 3; i < 10; i++ ) cout << identity_matrix(i).invert(); << endl;

Multiplying the result of invert() by diag(1 -1 1 -1 ... ) for n x n matrices where n is even and > 2 seems to provide the correct result, although I have not extensively tested it with more than very simple matrices. This suggests a bug in the cofactor calculation. I am looking at it, but if someone has some insight, I would be happy to have it.

[frenchjam]

By the way, I had to modify det() to get it to compile under Visual C++. It is entirely possible that I introduced a bug in doing so. I am going to try and track it down. But if you could tell me if the above works under gcc, for instance, I can be more confident that I am the one who introduced the problem in my code.

The modifications that I had to make to compile under VC++ 2017 had to do with the allocation of the permutation arrays p and v. The original code

uint p[n_] and v[n_]

would not compile, so I do:

uint *p = new uint[n_];
uint *v = new uint[n_];

...

delete p;
delete v;

I don't think that this could be the source of the issue with the inverse of higher order matrices, but who knows?