rational.py

"""Implementation of rational arithmetic."""

from __future__ import division

import math as _math

def _gcf(a, b):
    """Returns the greatest common factor of a and b."""
    a = abs(a)
    b = abs(b)
    while b:
        a, b = b, a % b
    return a

class Rational(object):
    """
    This class provides an exact representation of rational numbers.
 
    All of the standard arithmetic operators are provided.  In mixed-type
    expressions, an int or a long can be converted to a Rational without
    loss of precision, and will be done as such.

    Rationals can be implicity (using binary operators) or explicity
    (using float(x) or x.decimal()) converted to floats or Decimals;
    this may cause a loss of precision.  The reverse conversions can be
    done without loss of precision, and are performed with the
    from_exact_float and from_exact_decimal static methods.  However,
    because of rounding error in the original values, this tends to
    produce "ugly" fractions.  "Nicer" conversions to Rational can be made
    with approx_smallest_denominator or approx_smallest_error.
    """

    def __init__(self, numerator, denominator=1):
       """Contructs the Rational object for numerator/denominator."""
       if not isinstance(numerator, (int, long)):
           raise TypeError('numerator must have integer type')
       if not isinstance(denominator, (int, long)):
           raise TypeError('denominator must have integer type')
       if not denominator:
           raise ZeroDivisionError('rational construction')
       self._d = denominator
       self._n = numerator
       self.normalize_self()
    # Cancel the fraction to reduced form
    def normalize_self(self):
       factor = _gcf(self._n, self._d)
       self._n = self._n // factor
       self._d = self._d // factor
       if self._d < 0:
           self._n = -self._n
           self._d = -self._d

    def numerator(self):
        return self._n

    def denominator(self):
        return self._d

    def __repr__(self):
        if self._d == 1:
            return "Rational(%d)" % self._n
        else:
            return "Rational(%d, %d)" % (self._n, self._d)
    def __str__(self):
        if self._d == 1:
            return str(self._n)
        else:
            return "%d/%d" % (self._n, self._d)
    def __hash__(self):
        try:
            return hash(float(self))
        except OverflowError:
            return hash(long(self))
    def __float__(self):
        return self._n / self._d
    def __int__(self):
        if self._n < 0:
            return -int(-self._n // self._d)
        else:
            return int(self._n // self._d)
    def __long__(self):
        return long(int(self))
    def __nonzero__(self):
        return bool(self._n)
    def __pos__(self):
        return self
    def __neg__(self):
        return Rational(-self._n, self._d)
    def __abs__(self):
        if self._n < 0:
            return -self
        else:
            return self
    def __add__(self, other):
        if isinstance(other, Rational):
            return Rational(self._n * other._d + self._d * other._n,
                            self._d * other._d)
        elif isinstance(other, (int, long)):
            return Rational(self._n + self._d * other, self._d)
        elif isinstance(other, (float, complex)):
            return float(self) + other
        else:
            return NotImplemented
    __radd__ = __add__
    def __sub__(self, other):
        if isinstance(other, Rational):
            return Rational(self._n * other._d - self._d * other._n,
                            self._d * other._d)
        elif isinstance(other, (int, long)):
            return Rational(self._n - self._d * other, self._d)
        elif isinstance(other, (float, complex)):
            return float(self) - other
        else:
            return NotImplemented
    def __rsub__(self, other):
        if isinstance(other, (int, long)):
            return Rational(other * self._d - self._n, self._d)
        elif isinstance(other, (float, complex)):
            return other - float(self)
        else:
            return NotImplemented
    def __mul__(self, other):
        if isinstance(other, Rational):
            return Rational(self._n * other._n, self._d * other._d)
        elif isinstance(other, (int, long)):
            return Rational(self._n * other, self._d)
        elif isinstance(other, (float, complex)):
            return float(self) * other
        else:
            return NotImplemented
    __rmul__ = __mul__
    def __truediv__(self, other):
        if isinstance(other, Rational):
            return Rational(self._n * other._d, self._d * other._n)
        elif isinstance(other, (int, long)):
            return Rational(self._n, self._d * other)
        elif isinstance(other, (float, complex)):
            return float(self) / other
        else:
            return NotImplemented
    __div__ = __truediv__
    def __rtruediv__(self, other):
        if isinstance(other, (int, long)):
            return Rational(other * self._d, self._n)
        elif isinstance(other, (float, complex)):
            return other / float(self)
        else:
            return NotImplemented
    __rdiv__ = __rtruediv__
    def __floordiv__(self, other):
        truediv = self / other
        if isinstance(truediv, Rational):
            return truediv._n // truediv._d
        else:
            return truediv // 1
    def __rfloordiv__(self, other):
        return (other / self) // 1
    def __mod__(self, other):
        return self - self // other * other
    def __rmod__(self, other):
        return other - other // self * self
    def __divmod__(self, other):
        return self // other, self % other
    def __cmp__(self, other):
        if other == 0:
            return cmp(self._n, 0)
        else:
            return cmp(self - other, 0)
    def __pow__(self, other):
        if isinstance(other, (int, long)):
            if other < 0:
                return Rational(self._d ** -other, self._n ** -other)
            else:
                return Rational(self._n ** other, self._d ** other)
        else:
                return float(self) ** other
    def __rpow__(self, other):
        return other ** float(self)
    def round(self, denominator):
        """Return self rounded to nearest multiple of 1/denominator."""
        int_part, frac_part = divmod(self * denominator, 1)
        round_direction = cmp(frac_part * 2, 1)
        if round_direction == 0:
           numerator = int_part + (int_part & 1) # round to even
        elif round_direction < 0:
           numerator = int_part
        else:
           numerator = int_part + 1
        return Rational(numerator, denominator)



def rational_from_exact_float(x):
    """Returns the exact Rational equivalent of x."""
    mantissa, exponent = _math.frexp(x)
    mantissa = int(mantissa * 2 ** 53)
    exponent -= 53
    if exponent < 0:
        return Rational(mantissa, 2 ** (-exponent))
    else:
        return Rational(mantissa * 2 ** exponent)



def rational_approx_smallest_denominator(x, tolerance):
    """
    Returns a Rational approximation of x.
    Minimizes the denominator given a constraint on the error.

    x = the float or Decimal value to convert
    tolerance = maximum absolute error allowed,
                must be of the same type as x
    """
    tolerance = abs(tolerance)
    n = 1
    while True:
        m = int(round(x * n))
        result = Rational(m, n)
        if abs(result - x) < tolerance:
            return result
        n += 1


def rational_approx_smallest_error(x, maxDenominator):
    """
    Returns a Rational approximation of x.
    Minimizes the error given a constraint on the denominator.

    x = the float or Decimal value to convert
    maxDenominator = maximum denominator allowed
    """
    result = None
    minError = x
    for n in xrange(1, maxDenominator + 1):
        m = int(round(x * n))
        r = Rational(m, n)
        error = abs(r - x)
        if error == 0:
            return r
        elif error < minError:
            result = r
            minError = error
    return result

def divide(x, y):
    """Same as x/y, but returns a Rational if both are ints."""
    if isinstance(x, (int, long)) and isinstance(y, (int, long)):
        return Rational(x, y)
    else:
        return x / y