inverse of a matrix with Fraction entries

[Daniel Fetchinson]

>>>> 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'>

This sounds like a good idea!

Cheers,
Daniel


-- 
Psss, psss, put it down! - http://www.cafepress.com/putitdown