Double precision problem

[TheDestroyer]

I'm writing a class for matrices. In this class, I have a function for getting the inverse of the matrix. I use gaussian elimination (with same elementary row operation over an identity matrix) for obtaining the inverse.

The problem is, when I invert the following matrix:


2 2 5 4
5 1 7 8
9 13 11 12
13 14 15 16

I get the following result:


-0.2 -0.2 -0.6 0.6
1.11022e-016 -0.0909091 0.181818 -0.0909091
0.6 -0.309091 -0.381818 0.290909
-0.4 0.531818 0.686364 -0.618182

The element (2,1) is supposed to have the value 0 exactly, which I found on Mathematica. But there's a precision problem here apparently. My class is templatised over the variable type of each element. If I use long double, the element (2,1) becomes 0.

This is a real problem for me. I can't set any tolerance, because the tolerance may hinder some really small values to go, in case if some really small values are input in the matrix.

Any ideas about solutions for such a problem?

Thank you for any efforts :)

[EricDu]

As C/C++ native data types cannot support any all small numbers, you have to think about tolerance, or use third party library to do high precision calculation.

[helios]

These kinds of errors are inevitable when using floating point numbers by themselves. What you could do is use rationals instead, which don't have these sorts of problems because they can exactly represent more numbers than floating point values can.

[TheDestroyer]

Thank you for your answers.

@ErocDu: correction to what you said: double precision supports small numbers, since the exponent of the binary form could go to something like 2^-1024 (I'm not sure), but the problem is with the mantissa, that has a very limited precision.


Aren't there some common tricks to get over some situations where this could happen?