python/rational.py

   1 """Implementation of rational arithmetic."""
   2
   3 from __future__ import division
   4
   5 import math as _math
   6
   7 def _gcf(a, b):
   8     """Returns the greatest common factor of a and b."""
   9     a = abs(a)
  10     b = abs(b)
  11     while b:
  12         a, b = b, a % b
  13     return a
  14
  15 class Rational(object):
  16     """
  17     This class provides an exact representation of rational numbers.
  18
  19     All of the standard arithmetic operators are provided.  In mixed-type
  20     expressions, an int or a long can be converted to a Rational without
  21     loss of precision, and will be done as such.
  22
  23     Rationals can be implicity (using binary operators) or explicity
  24     (using float(x) or x.decimal()) converted to floats or Decimals;
  25     this may cause a loss of precision.  The reverse conversions can be
  26     done without loss of precision, and are performed with the
  27     from_exact_float and from_exact_decimal static methods.  However,
  28     because of rounding error in the original values, this tends to
  29     produce "ugly" fractions.  "Nicer" conversions to Rational can be made
  30     with approx_smallest_denominator or approx_smallest_error.
  31     """
  32
  33     def __init__(self, numerator, denominator=1):
  34        """Contructs the Rational object for numerator/denominator."""
  35        if not isinstance(numerator, (int, long)):
  36            raise TypeError('numerator must have integer type')
  37        if not isinstance(denominator, (int, long)):
  38            raise TypeError('denominator must have integer type')
  39        if not denominator:
  40            raise ZeroDivisionError('rational construction')
  41        # Store the fraction in reduced form as _n/_d
  42        factor = _gcf(numerator, denominator)
  43        self._n = numerator // factor
  44        self._d = denominator // factor
  45        if self._d < 0:
  46            self._n = -self._n
  47            self._d = -self._d
  48     def numerator(self):
  49         return self._n
  50
  51     def denominator(self):
  52         return self._d
  53
  54     def __repr__(self):
  55         if self._d == 1:
  56             return "Rational(%d)" % self._n
  57         else:
  58             return "Rational(%d, %d)" % (self._n, self._d)
  59     def __str__(self):
  60         if self._d == 1:
  61             return str(self._n)
  62         else:
  63             return "%d/%d" % (self._n, self._d)
  64     def __hash__(self):
  65         try:
  66             return hash(float(self))
  67         except OverflowError:
  68             return hash(long(self))
  69     def __float__(self):
  70         return self._n / self._d
  71     def __int__(self):
  72         if self._n < 0:
  73             return -int(-self._n // self._d)
  74         else:
  75             return int(self._n // self._d)
  76     def __long__(self):
  77         return long(int(self))
  78     def __nonzero__(self):
  79         return bool(self._n)
  80     def __pos__(self):
  81         return self
  82     def __neg__(self):
  83         return Rational(-self._n, self._d)
  84     def __abs__(self):
  85         if self._n < 0:
  86             return -self
  87         else:
  88             return self
  89     def __add__(self, other):
  90         if isinstance(other, Rational):
  91             return Rational(self._n * other._d + self._d * other._n,
  92                             self._d * other._d)
  93         elif isinstance(other, (int, long)):
  94             return Rational(self._n + self._d * other, self._d)
  95         elif isinstance(other, (float, complex)):
  96             return float(self) + other
  97         else:
  98             return NotImplemented
  99     __radd__ = __add__
 100     def __sub__(self, other):
 101         if isinstance(other, Rational):
 102             return Rational(self._n * other._d - self._d * other._n,
 103                             self._d * other._d)
 104         elif isinstance(other, (int, long)):
 105             return Rational(self._n - self._d * other, self._d)
 106         elif isinstance(other, (float, complex)):
 107             return float(self) - other
 108         else:
 109             return NotImplemented
 110     def __rsub__(self, other):
 111         if isinstance(other, (int, long)):
 112             return Rational(other * self._d - self._n, self._d)
 113         elif isinstance(other, (float, complex)):
 114             return other - float(self)
 115         else:
 116             return NotImplemented
 117     def __mul__(self, other):
 118         if isinstance(other, Rational):
 119             return Rational(self._n * other._n, self._d * other._d)
 120         elif isinstance(other, (int, long)):
 121             return Rational(self._n * other, self._d)
 122         elif isinstance(other, (float, complex)):
 123             return float(self) * other
 124         else:
 125             return NotImplemented
 126     __rmul__ = __mul__
 127     def __truediv__(self, other):
 128         if isinstance(other, Rational):
 129             return Rational(self._n * other._d, self._d * other._n)
 130         elif isinstance(other, (int, long)):
 131             return Rational(self._n, self._d * other)
 132         elif isinstance(other, (float, complex)):
 133             return float(self) / other
 134         else:
 135             return NotImplemented
 136     __div__ = __truediv__
 137     def __rtruediv__(self, other):
 138         if isinstance(other, (int, long)):
 139             return Rational(other * self._d, self._n)
 140         elif isinstance(other, (float, complex)):
 141             return other / float(self)
 142         else:
 143             return NotImplemented
 144     __rdiv__ = __rtruediv__
 145     def __floordiv__(self, other):
 146         truediv = self / other
 147         if isinstance(truediv, Rational):
 148             return truediv._n // truediv._d
 149         else:
 150             return truediv // 1
 151     def __rfloordiv__(self, other):
 152         return (other / self) // 1
 153     def __mod__(self, other):
 154         return self - self // other * other
 155     def __rmod__(self, other):
 156         return other - other // self * self
 157     def __divmod__(self, other):
 158         return self // other, self % other
 159     def __cmp__(self, other):
 160         if other == 0:
 161             return cmp(self._n, 0)
 162         else:
 163             return cmp(self - other, 0)
 164     def __pow__(self, other):
 165         if isinstance(other, (int, long)):
 166             if other < 0:
 167                 return Rational(self._d ** -other, self._n ** -other)
 168             else:
 169                 return Rational(self._n ** other, self._d ** other)
 170         else:
 171                 return float(self) ** other
 172     def __rpow__(self, other):
 173         return other ** float(self)
 174     def round(self, denominator):
 175         """Return self rounded to nearest multiple of 1/denominator."""
 176         int_part, frac_part = divmod(self * denominator, 1)
 177         round_direction = cmp(frac_part * 2, 1)
 178         if round_direction == 0:
 179            numerator = int_part + (int_part & 1) # round to even
 180         elif round_direction < 0:
 181            numerator = int_part
 182         else:
 183            numerator = int_part + 1
 184         return Rational(numerator, denominator)
 185
 186
 187
 188 def rational_from_exact_float(x):
 189     """Returns the exact Rational equivalent of x."""
 190     mantissa, exponent = _math.frexp(x)
 191     mantissa = int(mantissa * 2 ** 53)
 192     exponent -= 53
 193     if exponent < 0:
 194         return Rational(mantissa, 2 ** (-exponent))
 195     else:
 196         return Rational(mantissa * 2 ** exponent)
 197
 198
 199
 200 def rational_approx_smallest_denominator(x, tolerance):
 201     """
 202     Returns a Rational approximation of x.
 203     Minimizes the denominator given a constraint on the error.
 204
 205     x = the float or Decimal value to convert
 206     tolerance = maximum absolute error allowed,
 207                 must be of the same type as x
 208     """
 209     tolerance = abs(tolerance)
 210     n = 1
 211     while True:
 212         m = int(round(x * n))
 213         result = Rational(m, n)
 214         if abs(result - x) < tolerance:
 215             return result
 216         n += 1
 217
 218
 219 def rational_approx_smallest_error(x, maxDenominator):
 220     """
 221     Returns a Rational approximation of x.
 222     Minimizes the error given a constraint on the denominator.
 223
 224     x = the float or Decimal value to convert
 225     maxDenominator = maximum denominator allowed
 226     """
 227     result = None
 228     minError = x
 229     for n in xrange(1, maxDenominator + 1):
 230         m = int(round(x * n))
 231         r = Rational(m, n)
 232         error = abs(r - x)
 233         if error == 0:
 234             return r
 235         elif error < minError:
 236             result = r
 237             minError = error
 238     return result
 239
 240 def divide(x, y):
 241     """Same as x/y, but returns a Rational if both are ints."""
 242     if isinstance(x, (int, long)) and isinstance(y, (int, long)):
 243         return Rational(x, y)
 244     else:
 245         return x / y
 246