python/rational.py

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