[GiNaC-list] Inverse of complex matrix

[Richard B. Kreckel]

On Thu, 4 Nov 2004, Jonathan Brandmeyer wrote:
> Sometimes, inverting a matrix with complex values does not produce
> something that, when multiplied by the original, returns the identity.

In a way, it does.  It just isn't apparent that it is the identity.  :-(

> I have attached a program that demonstrates this for one case, which was
> found by altering checks/check_matricies.cpp to use complex values from
> sparse_tree().
>
> Is this a known problem?

Yes, it is a known issue.  The problem is that normal(), as it is
implemented, cannot deal with polynomials containing complex coefficients.
It's easy to construct an even simpler failure.  In ginsh:

> normal((x^2+1)/((x+I)*(x-I)));
(I+x)^(-1)*(1+x^2)*(-I+x)^(-1)

So, normal() returned the input unchanged although it is equal to 1.  The
problem you are describing has nothing to do with the matrix code proper.

Algebraic extension fields are something that hasn't been touched by GiNaC
till now.  Note also, that in GiNaC, the imaginary unit `I' is not a
symbol, as in some other systems but a number.  (I do not know what `I' as
a symbol could buy you anyways, because eventually you *will* transform
x^2-I^2 to x^2+1 or even have x^2+1 as input and then the knowledge that
the 1 came from -I^2 is completely gone.  So, some algorithm capable of
dealing with extension fields has to be employed in anyway.)

Regards
     -richy.
-- 
Richard B. Kreckel
<http://www.ginac.de/~kreckel/>