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         else:
  99             return NotImplemented
 100     __radd__ = __add__
 101     def __sub__(self, other):
 102         if isinstance(other, Rational):
 103             return Rational(self._n * other._d - self._d * other._n,
 104                             self._d * other._d)
 105         elif isinstance(other, (int, long)):
 106             return Rational(self._n - self._d * other, self._d)
 107         elif isinstance(other, (float, complex)):
 108             return float(self) - other
 109         else:
 110             return NotImplemented
 111     def __rsub__(self, other):
 112         if isinstance(other, (int, long)):
 113             return Rational(other * self._d - self._n, self._d)
 114         elif isinstance(other, (float, complex)):
 115             return other - float(self)
 116         else:
 117             return NotImplemented
 118     def __mul__(self, other):
 119         if isinstance(other, Rational):
 120             return Rational(self._n * other._n, self._d * other._d)
 121         elif isinstance(other, (int, long)):
 122             return Rational(self._n * other, self._d)
 123         elif isinstance(other, (float, complex)):
 124             return float(self) * other
 125         else:
 126             return NotImplemented
 127     __rmul__ = __mul__
 128     def __truediv__(self, other):
 129         if isinstance(other, Rational):
 130             return Rational(self._n * other._d, self._d * other._n)
 131         elif isinstance(other, (int, long)):
 132             return Rational(self._n, self._d * other)
 133         elif isinstance(other, (float, complex)):
 134             return float(self) / other
 135         else:
 136             return NotImplemented
 137     __div__ = __truediv__
 138     def __rtruediv__(self, other):
 139         if isinstance(other, (int, long)):
 140             return Rational(other * self._d, self._n)
 141         elif isinstance(other, (float, complex)):
 142             return other / float(self)
 143         else:
 144             return NotImplemented
 145     __rdiv__ = __rtruediv__
 146     def __floordiv__(self, other):
 147         truediv = self / other
 148         if isinstance(truediv, Rational):
 149             return truediv._n // truediv._d
 150         else:
 151             return truediv // 1
 152     def __rfloordiv__(self, other):
 153         return (other / self) // 1
 154     def __mod__(self, other):
 155         return self - self // other * other
 156     def __rmod__(self, other):
 157         return other - other // self * self
 158     def __divmod__(self, other):
 159         return self // other, self % other
 160     def __cmp__(self, other):
 161         if other == 0:
 162             return cmp(self._n, 0)
 163         else:
 164             return cmp(self - other, 0)
 165     def __pow__(self, other):
 166         if isinstance(other, (int, long)):
 167             if other < 0:
 168                 return Rational(self._d ** -other, self._n ** -other)
 169             else:
 170                 return Rational(self._n ** other, self._d ** other)
 171         else:
 172                 return float(self) ** other
 173     def __rpow__(self, other):
 174         return other ** float(self)
 175     def round(self, denominator):
 176         """Return self rounded to nearest multiple of 1/denominator."""
 177         int_part, frac_part = divmod(self * denominator, 1)
 178         round_direction = cmp(frac_part * 2, 1)
 179         if round_direction == 0:
 180            numerator = int_part + (int_part & 1) # round to even
 181         elif round_direction < 0:
 182            numerator = int_part
 183         else:
 184            numerator = int_part + 1
 185         return Rational(numerator, denominator)
 186
 187
 188
 189 def rational_from_exact_float(x):
 190     """Returns the exact Rational equivalent of x."""
 191     mantissa, exponent = _math.frexp(x)
 192     mantissa = int(mantissa * 2 ** 53)
 193     exponent -= 53
 194     if exponent < 0:
 195         return Rational(mantissa, 2 ** (-exponent))
 196     else:
 197         return Rational(mantissa * 2 ** exponent)
 198
 199
 200
 201 def rational_approx_smallest_denominator(x, tolerance):
 202     """
 203     Returns a Rational approximation of x.
 204     Minimizes the denominator given a constraint on the error.
 205
 206     x = the float or Decimal value to convert
 207     tolerance = maximum absolute error allowed,
 208                 must be of the same type as x
 209     """
 210     tolerance = abs(tolerance)
 211     n = 1
 212     while True:
 213         m = int(round(x * n))
 214         result = Rational(m, n)
 215         if abs(result - x) < tolerance:
 216             return result
 217         n += 1
 218
 219
 220 def rational_approx_smallest_error(x, maxDenominator):
 221     """
 222     Returns a Rational approximation of x.
 223     Minimizes the error given a constraint on the denominator.
 224
 225     x = the float or Decimal value to convert
 226     maxDenominator = maximum denominator allowed
 227     """
 228     result = None
 229     minError = x
 230     for n in xrange(1, maxDenominator + 1):
 231         m = int(round(x * n))
 232         r = Rational(m, n)
 233         error = abs(r - x)
 234         if error == 0:
 235             return r
 236         elif error < minError:
 237             result = r
 238             minError = error
 239     return result
 240
 241 def divide(x, y):
 242     """Same as x/y, but returns a Rational if both are ints."""
 243     if isinstance(x, (int, long)) and isinstance(y, (int, long)):
 244         return Rational(x, y)
 245     else:
 246         return x / y
 247