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        self._d = denominator
  42        self._n = numerator
  43        self.normalize_self()
  44     # Cancel the fraction to reduced form
  45     def normalize_self(self):
  46        factor = _gcf(self._n, self._d)
  47        self._n = self._n // factor
  48        self._d = self._d // factor
  49        if self._d < 0:
  50            self._n = -self._n
  51            self._d = -self._d
  52
  53     def numerator(self):
  54         return self._n
  55
  56     def denominator(self):
  57         return self._d
  58
  59     def __repr__(self):
  60         if self._d == 1:
  61             return "Rational(%d)" % self._n
  62         else:
  63             return "Rational(%d, %d)" % (self._n, self._d)
  64     def __str__(self):
  65         if self._d == 1:
  66             return str(self._n)
  67         else:
  68             return "%d/%d" % (self._n, self._d)
  69     def __hash__(self):
  70         try:
  71             return hash(float(self))
  72         except OverflowError:
  73             return hash(long(self))
  74     def __float__(self):
  75         return self._n / self._d
  76     def __int__(self):
  77         if self._n < 0:
  78             return -int(-self._n // self._d)
  79         else:
  80             return int(self._n // self._d)
  81     def __long__(self):
  82         return long(int(self))
  83     def __nonzero__(self):
  84         return bool(self._n)
  85     def __pos__(self):
  86         return self
  87     def __neg__(self):
  88         return Rational(-self._n, self._d)
  89     def __abs__(self):
  90         if self._n < 0:
  91             return -self
  92         else:
  93             return self
  94     def __add__(self, other):
  95         if isinstance(other, Rational):
  96             return Rational(self._n * other._d + self._d * other._n,
  97                             self._d * other._d)
  98         elif isinstance(other, (int, long)):
  99             return Rational(self._n + self._d * other, self._d)
 100         elif isinstance(other, (float, complex)):
 101             return float(self) + other
 102         else:
 103             return NotImplemented
 104     __radd__ = __add__
 105     def __sub__(self, other):
 106         if isinstance(other, Rational):
 107             return Rational(self._n * other._d - self._d * other._n,
 108                             self._d * other._d)
 109         elif isinstance(other, (int, long)):
 110             return Rational(self._n - self._d * other, self._d)
 111         elif isinstance(other, (float, complex)):
 112             return float(self) - other
 113         else:
 114             return NotImplemented
 115     def __rsub__(self, other):
 116         if isinstance(other, (int, long)):
 117             return Rational(other * self._d - self._n, self._d)
 118         elif isinstance(other, (float, complex)):
 119             return other - float(self)
 120         else:
 121             return NotImplemented
 122     def __mul__(self, other):
 123         if isinstance(other, Rational):
 124             return Rational(self._n * other._n, self._d * other._d)
 125         elif isinstance(other, (int, long)):
 126             return Rational(self._n * other, self._d)
 127         elif isinstance(other, (float, complex)):
 128             return float(self) * other
 129         else:
 130             return NotImplemented
 131     __rmul__ = __mul__
 132     def __truediv__(self, other):
 133         if isinstance(other, Rational):
 134             return Rational(self._n * other._d, self._d * other._n)
 135         elif isinstance(other, (int, long)):
 136             return Rational(self._n, self._d * other)
 137         elif isinstance(other, (float, complex)):
 138             return float(self) / other
 139         else:
 140             return NotImplemented
 141     __div__ = __truediv__
 142     def __rtruediv__(self, other):
 143         if isinstance(other, (int, long)):
 144             return Rational(other * self._d, self._n)
 145         elif isinstance(other, (float, complex)):
 146             return other / float(self)
 147         else:
 148             return NotImplemented
 149     __rdiv__ = __rtruediv__
 150     def __floordiv__(self, other):
 151         truediv = self / other
 152         if isinstance(truediv, Rational):
 153             return truediv._n // truediv._d
 154         else:
 155             return truediv // 1
 156     def __rfloordiv__(self, other):
 157         return (other / self) // 1
 158     def __mod__(self, other):
 159         return self - self // other * other
 160     def __rmod__(self, other):
 161         return other - other // self * self
 162     def __divmod__(self, other):
 163         return self // other, self % other
 164     def __cmp__(self, other):
 165         if other == 0:
 166             return cmp(self._n, 0)
 167         else:
 168             return cmp(self - other, 0)
 169     def __pow__(self, other):
 170         if isinstance(other, (int, long)):
 171             if other < 0:
 172                 return Rational(self._d ** -other, self._n ** -other)
 173             else:
 174                 return Rational(self._n ** other, self._d ** other)
 175         else:
 176                 return float(self) ** other
 177     def __rpow__(self, other):
 178         return other ** float(self)
 179     def round(self, denominator):
 180         """Return self rounded to nearest multiple of 1/denominator."""
 181         int_part, frac_part = divmod(self * denominator, 1)
 182         round_direction = cmp(frac_part * 2, 1)
 183         if round_direction == 0:
 184            numerator = int_part + (int_part & 1) # round to even
 185         elif round_direction < 0:
 186            numerator = int_part
 187         else:
 188            numerator = int_part + 1
 189         return Rational(numerator, denominator)
 190
 191
 192
 193 def rational_from_exact_float(x):
 194     """Returns the exact Rational equivalent of x."""
 195     mantissa, exponent = _math.frexp(x)
 196     mantissa = int(mantissa * 2 ** 53)
 197     exponent -= 53
 198     if exponent < 0:
 199         return Rational(mantissa, 2 ** (-exponent))
 200     else:
 201         return Rational(mantissa * 2 ** exponent)
 202
 203
 204
 205 def rational_approx_smallest_denominator(x, tolerance):
 206     """
 207     Returns a Rational approximation of x.
 208     Minimizes the denominator given a constraint on the error.
 209
 210     x = the float or Decimal value to convert
 211     tolerance = maximum absolute error allowed,
 212                 must be of the same type as x
 213     """
 214     tolerance = abs(tolerance)
 215     n = 1
 216     while True:
 217         m = int(round(x * n))
 218         result = Rational(m, n)
 219         if abs(result - x) < tolerance:
 220             return result
 221         n += 1
 222
 223
 224 def rational_approx_smallest_error(x, maxDenominator):
 225     """
 226     Returns a Rational approximation of x.
 227     Minimizes the error given a constraint on the denominator.
 228
 229     x = the float or Decimal value to convert
 230     maxDenominator = maximum denominator allowed
 231     """
 232     result = None
 233     minError = x
 234     for n in xrange(1, maxDenominator + 1):
 235         m = int(round(x * n))
 236         r = Rational(m, n)
 237         error = abs(r - x)
 238         if error == 0:
 239             return r
 240         elif error < minError:
 241             result = r
 242             minError = error
 243     return result
 244
 245 def divide(x, y):
 246     """Same as x/y, but returns a Rational if both are ints."""
 247     if isinstance(x, (int, long)) and isinstance(y, (int, long)):
 248         return Rational(x, y)
 249     else:
 250         return x / y
 251