inverse of a matrix with Fraction entries

[Peter Otten]

Daniel Fetchinson wrote:

>>> I guess this is a question to folks with some numpy background (but
>>> not necessarily).
>>>
>>> I'm using fractions.Fraction as entries in a matrix because I need to
>>> have very high precision and fractions.Fraction provides infinite
>>> precision (as I've learned from advice from this list).
>>
>> "Infinite" precision only for values that can be expressed as a fraction:
>>
>>>>> Fraction(2)**Fraction(1, 2)
>> 1.4142135623730951
>>>>> type(_)
>> <type 'float'>
> 
> True, but I only need to add, multiply and divide my fractions in
> order to end up with the entries in the matrix.
> 
>>> Now I need to
>>> calculate its inverse. Can numpy help in this regard? Can I tell numpy
>>> that the inverse matrix should also have entries in
>>> fractions.Fraction? Or numpy can only do floating point calculations?
>>
>> I tried it, and numpy.linalg.inv() converted the Fraction values to
>> float. If you aren't concerned about efficiency it should be easy to do
>> it yourself, in pure python.
> 
> If there is no other way, that's what I'll try to do.
> 
>> You could also ask on the numpy mailing list.
>>
>>> Probably it doesn't matter but the matrix has all components non-zero
>>> and is about a thousand by thousand in size.
>>
>> Hmm, where did you get that million of exact fractions from?
> 
> The million of exact fractions are coming from a recursion relation.
> This recursion relation only adds, multiplies and divides numbers so
> the end result is always a rational number.
> 
>> If some real
>> world data is involved the error introduced by the floating point
>> calculation may be negligable in comparison with the initial measurement
>> uncertainty.
> 
> It's a mathematical problem so no uncertainty is present in the
> initial values. And even if there was, if there are many orders of
> magnitude differences between the entries in the matrix floating point
> does not suffice for various things like eigenvalue calculation and
> stuff like that.

It may be worthwhile to have a look at http://www.sagemath.org/

sage: Matrix([[1,2],[3,4]])**-1

[  -2    1]
[ 3/2 -1/2]
sage: a = Matrix([[1,2],[3,4]])
sage: b = Matrix([[1,2],[3,4]])**-1
sage: a*b

[1 0]
[0 1]
sage: type(b[1,1])
<type 'sage.rings.rational.Rational'>

Peter