Fractions with decimal precision

Question

Is there a pure python implementation of fractions.Fraction that supports longs as numerator and denominator? Unfortunately, exponentiation appears to be coded in to return a float (ack!!!), which should at least support using decimal.Decimal.

If there isn't, I suppose I can probably make a copy of the library and try to replace occurrences of float() with something appropriate from Decimal but I'd rather something that's been tested by others before.

Here's a code example:

base = Fraction.from_decimal(Decimal(1).exp())
a = Fraction(69885L, 53L)
x = Fraction(9L, 10L)

print base**(-a*x), type(base**(-a*x))

results in 0.0 <type 'float'> where the answer should be a really small decimal.

Update: I've got the following work-around for now (assuming, for a**b, that both are fractions; of course, I'll need another function when exp_ is a float or is itself a Decimal):

def fracpow(base, exp_):
    base = Decimal(base.numerator)/Decimal(base.denominator)
    exp_ = Decimal(exp_.numerator)/Decimal(exp_.denominator)

    return base**exp_

which gives the answer 4.08569925773896097019795484811E-516.

I'd still be interested if there's a better way of doing this without the extra functions (I'm guessing if I work with the Fraction class enough, I'll find other floats working their way into my results).

Ovine

When I do fractions.Fraction(11,17) ** 2 I get fractions.Fraction(121, 289). What do you get? – Dagoba

try type(fractions.Fraction(11,17)**2.1) – Ovine

IMO, the behavior you see is correct, although not what you want. Since the exponent is a float, it is reasonable that things are converted to float. Now, if the exponent was Fraction(21,10), that would be more interesting... – Dagoba

See the updated code above, where the exponent is in fact a Fraction. – Ovine

Answer 1

"Raise to a power" is not a closed operation over the rationals (differently from the usual four arithmetic operations): there is no rational number r such that r == 2 ** 0.5. Legend has it that Pythagoras (from whose theorem this fact so simply follows) had his disciple Hippasus killed for the horrible crime of proving this; looks like you sympathize wit Pythagoras' alleged reaction;-), given your weird use of "should".

Python's fractions are meant to be exact, so inevitably there are case in which raising a fraction to another fraction's power will be absolutely unable to return a fraction as its result; and "should" just cannot be sensibly applied to a mathematical impossibility.

So the best you can do is to approximate your desired result, e.g. by getting a result that's not an exact fraction (floats are generally considered sufficient for the purpose) and then further approximating it back with a fraction. Most existing pure-Python implementations (there are many rationals.py files found around the net;-) prefer not to implement a ** operator at all, but of course there's nothing stopping you from making a different design decision in your own implementation!-)

Supersede

This is of course correct mathematically, but the reason I'm using the Fraction class is because I want high numerical accuracy, and I'm trying to update code someone else wrote using the clnum library. Since Fraction's happy to take a long as a numerator or denominator, it seems only fair that it should be able to return a Decimal upon exponentiation. Now, if only I could figure out if fracpow(0,0) should be defined as 1, or if that was an error in the code! – Ovine

@Noah, can't really help you with pure Python code here, since for such tasks I always use gmpy, of course; but gmpy.mpq does not allow inexact-roots with fractional exponents either (my design decision in this case), so I'd go through high-precision floats (gmpy.mpf in my case) just like you're doing with Decimal (which are floats too, just decimal rather than binary ones, with SW rather than HW implementations, and with precision settable as high as desired -- gmpy.mpf are similar, but binary rather than decimal). – Supersede

The error from your numerical approximation of exp will be usually bigger then the error from approximating the approximation result with a float. – Tindall

@quant_dev, I think you're wrong: a plain float has about 18 digits of precision, while by default a Decimal has 28 (and decimal.Decimal(1).exp() is indeed precise to 28 digits) -- plus of course you can always use a decimal.Context(prec=99) to force the precision up to (say) 99 digits (not an option for plain floats;-). – Supersede

If 2 ** 0.5 justifies using float for all cases instead of fractions.Fraction then why (-1) ** 0.5 doesn't justify using complex for all cases. – Breastbone

@Alex Martelli, but the question is: can you compute exp(x) with correct 18 digits? – Tindall

If I were, hypothetically, going to have our sysadmin put one of these packages on our cluster, how does gmpy compare to clnum? Your unbiased opinion, of course :) – Ovine

@quant_dev, of course you can compute Decimal(n).exp(context) to whatever prec the context specifies -- what makes you think otherwise? Check out the algorithm used in svn.python.org/view/python/trunk/Lib/… (lines 2819-2892) -- seems correct to me. – Supersede

@Noah, clnum is rich in transcendental functions (trig &c), and complex numbers; gmpy does not offer those but has more support for integers (fibonacci, gcd, lcm, primality testing, bit-level functions such as hamming distance and population count, jacobi, kronecker, legendre ...) and fractions ("Stern-Brocot tree" approximation). – Supersede

@J.F., the answer to your "why" is: fractions.Fraction is about exact computation, float (and decimal.Decimal) are not (so they do regularly supply approximate answers, e.g. as results of **). – Supersede

Answer 2

You can write your own "pow" function for fractions that doesn't use floating-point exponentiation. Is that what you're trying to do?

This will raise a fraction to an integer power with falling back to float.

def pow( fract, exp ):
    if exp == 0: 
        return fract
    elif exp % 2 == 0:
        t = pow( fract, exp//2 )
        return t*t
    else:
        return fract*pos( fract, exp-1 )

Stasiastasis

Hmmm, better have some base case to terminate the recursion eventually!-) – Supersede