views:

286

answers:

2

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).

A: 

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 )
S.Lott
Hmmm, better have _some_ base case to terminate the recursion eventually!-)
Alex Martelli
Thanks!......15
S.Lott
+5  A: 

"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!-)

Alex Martelli
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!
Noah
@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).
Alex Martelli
The error from your numerical approximation of exp will be usually bigger then the error from approximating the approximation result with a float.
quant_dev
@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;-).
Alex Martelli
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.
J.F. Sebastian
@Alex Martelli, but the question is: can you compute exp(x) with correct 18 digits?
quant_dev
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 :)
Noah
Alex Martelli
@Noah, clnum is rich in transcendental functions (trig 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).
Alex Martelli
@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 `**`).
Alex Martelli