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
 202
 203
 204 def rational_from_exact_float(x):
 205     """Returns the exact Rational equivalent of x."""
 206     mantissa, exponent = _math.frexp(x)
 207     mantissa = int(mantissa * 2 ** 53)
 208     exponent -= 53
 209     if exponent < 0:
 210         return Rational(mantissa, 2 ** (-exponent))
 211     else:
 212         return Rational(mantissa * 2 ** exponent)
 213
 214 def rational_from_exact_decimal(x):
 215     """Returns the exact Rational equivalent of x."""
 216     sign, mantissa, exponent = x.as_tuple()
 217     sign = (1, -1)[sign]
 218     mantissa = sign * reduce(lambda a, b: 10 * a + b, mantissa)
 219     if exponent < 0:
 220         return Rational(mantissa, 10 ** (-exponent))
 221     else:
 222         return Rational(mantissa * 10 ** exponent)
 223
 224
 225 def rational_approx_smallest_denominator(x, tolerance):
 226     """
 227     Returns a Rational approximation of x.
 228     Minimizes the denominator given a constraint on the error.
 229
 230     x = the float or Decimal value to convert
 231     tolerance = maximum absolute error allowed,
 232                 must be of the same type as x
 233     """
 234     tolerance = abs(tolerance)
 235     n = 1
 236     while True:
 237         m = int(round(x * n))
 238         result = Rational(m, n)
 239         if abs(result - x) < tolerance:
 240             return result
 241         n += 1
 242
 243
 244 def rational_approx_smallest_error(x, maxDenominator):
 245     """
 246     Returns a Rational approximation of x.
 247     Minimizes the error given a constraint on the denominator.
 248
 249     x = the float or Decimal value to convert
 250     maxDenominator = maximum denominator allowed
 251     """
 252     result = None
 253     minError = x
 254     for n in xrange(1, maxDenominator + 1):
 255         m = int(round(x * n))
 256         r = Rational(m, n)
 257         error = abs(r - x)
 258         if error == 0:
 259             return r
 260         elif error < minError:
 261             result = r
 262             minError = error
 263     return result
 264
 265 def divide(x, y):
 266     """Same as x/y, but returns a Rational if both are ints."""
 267     if isinstance(x, (int, long)) and isinstance(y, (int, long)):
 268         return Rational(x, y)
 269     else:
 270         return x / y
 271