precision problems in base conversion of rational numbers

[Dan Bishop]

Brian van den Broek wrote:
> Hi all,
>
> I guess it is more of a maths question than a programming one, but it
> involves use of the decimal module, so here goes:
>
> As a self-directed learning exercise I've been working on a script to
> convert numbers to arbitrary bases. It aims to take any of whole
> numbers (python ints, longs, or Decimals), rational numbers (n / m n,
> m whole) and floating points (the best I can do for reals), and
> convert them to any base between 2 and 36, to desired precision.
...
> I've discovered empirically that I have to set getcontext.prec = 99 or
> greater (i.e. x >= 1.2375) to get
> 0.7CF0CA7CF0CA7CF0CA7CF0CA7CF0CA7CF0CA7CF0CA7CF0CA7CF0CA7CF0CA7CF0CA7CF0CA7CF0CA7C
>
> (I do feel safe in assuming that is the right answer. :-)

And you are correct.

~$ repdigits 345/756 17
0.(7CF0CA)

(That's from a script I wrote a couple of years ago.  For your next
self-directed learning exercise, you could get your program to produce
that output :-)

Anyhow, if you need such precise precision, perhaps you should consider
doing an exact calculation, with rational numbers.





"""Implementation of rational arithmetic."""

from __future__ import division

import decimal as _decimal
import math as _math

def _gcf(a, b):
    """Returns the greatest common factor of a and b."""
    a = abs(a)
    b = abs(b)
    while b:
        a, b = b, a % b
    return a

class Rational(object):
    """
    This class provides an exact representation of rational numbers.

    All of the standard arithmetic operators are provided.  In
mixed-type
    expressions, an int or a long can be converted to a Rational
without
    loss of precision, and will be done as such.

    Rationals can be implicity (using binary operators) or explicity
    (using float(x) or x.decimal()) converted to floats or Decimals;
    this may cause a loss of precision.  The reverse conversions can be
    done without loss of precision, and are performed with the
    from_exact_float and from_exact_decimal static methods.  However,
    because of rounding error in the original values, this tends to
    produce "ugly" fractions.  "Nicer" conversions to Rational can be
made
    with approx_smallest_denominator or approx_smallest_error.
    """
    def __init__(self, numerator, denominator=1):
       """Contructs the Rational object for numerator/denominator."""
       if not isinstance(numerator, (int, long)):
           raise TypeError('numerator must have integer type')
       if not isinstance(denominator, (int, long)):
           raise TypeError('denominator must have integer type')
       if not denominator:
           raise ZeroDivisionError('rational construction')
       # Store the fraction in reduced form as _n/_d
       factor = _gcf(numerator, denominator)
       self._n = numerator // factor
       self._d = denominator // factor
       if self._d < 0:
           self._n = -self._n
           self._d = -self._d
    def __repr__(self):
        if self._d == 1:
            return "Rational(%d)" % self._n
        else:
            return "Rational(%d, %d)" % (self._n, self._d)
    def __str__(self):
        if self._d == 1:
            return str(self._n)
        else:
            return "%d/%d" % (self._n, self._d)
    def __hash__(self):
        try:
            return hash(float(self))
        except OverflowError:
            return hash(long(self))
    def __float__(self):
        return self._n / self._d
    def __int__(self):
        if self._n < 0:
            return -int(-self._n // self._d)
        else:
            return int(self._n // self._d)
    def __long__(self):
        return long(int(self))
    def __nonzero__(self):
        return bool(self._n)
    def __pos__(self):
        return self
    def __neg__(self):
        return Rational(-self._n, self._d)
    def __abs__(self):
        if self._n < 0:
            return -self
        else:
            return self
    def __add__(self, other):
        if isinstance(other, Rational):
            return Rational(self._n * other._d + self._d * other._n,
                            self._d * other._d)
        elif isinstance(other, (int, long)):
            return Rational(self._n + self._d * other, self._d)
        elif isinstance(other, (float, complex)):
            return float(self) + other
        elif isinstance(other, _decimal.Decimal):
            return self.decimal() + other
        else:
            return NotImplemented
    __radd__ = __add__
    def __sub__(self, other):
        if isinstance(other, Rational):
            return Rational(self._n * other._d - self._d * other._n,
                            self._d * other._d)
        elif isinstance(other, (int, long)):
            return Rational(self._n - self._d * other, self._d)
        elif isinstance(other, (float, complex)):
            return float(self) - other
        elif isinstance(other, _decimal.Decimal):
            return self.decimal() - other
        else:
            return NotImplemented
    def __rsub__(self, other):
        if isinstance(other, (int, long)):
            return Rational(other * self._d - self._n, self._d)
        elif isinstance(other, (float, complex)):
            return other - float(self)
        elif isinstance(other, _decimal.Decimal):
            return other - self.decimal()
        else:
            return NotImplemented
    def __mul__(self, other):
        if isinstance(other, Rational):
            return Rational(self._n * other._n, self._d * other._d)
        elif isinstance(other, (int, long)):
            return Rational(self._n * other, self._d)
        elif isinstance(other, (float, complex)):
            return float(self) * other
        elif isinstance(other, _decimal.Decimal):
            return self.decimal() * other
        else:
            return NotImplemented
    __rmul__ = __mul__
    def __truediv__(self, other):
        if isinstance(other, Rational):
            return Rational(self._n * other._d, self._d * other._n)
        elif isinstance(other, (int, long)):
            return Rational(self._n, self._d * other)
        elif isinstance(other, (float, complex)):
            return float(self) / other
        elif isinstance(other, _decimal.Decimal):
            return self.decimal() / other
        else:
            return NotImplemented
    __div__ = __truediv__
    def __rtruediv__(self, other):
        if isinstance(other, (int, long)):
            return Rational(other * self._d, self._n)
        elif isinstance(other, (float, complex)):
            return other / float(self)
        elif isinstance(other, _decimal.Decimal):
            return other / self.decimal()
        else:
            return NotImplemented
    __rdiv__ = __rtruediv__
    def __floordiv__(self, other):
        truediv = self / other
        if isinstance(truediv, Rational):
            return truediv._n // truediv._d
        else:
            return truediv // 1
    def __rfloordiv__(self, other):
        return (other / self) // 1
    def __mod__(self, other):
        return self - self // other * other
    def __rmod__(self, other):
        return other - other // self * self
    def __divmod__(self, other):
        return self // other, self % other
    def __cmp__(self, other):
        if other == 0:
            return cmp(self._n, 0)
        else:
            return cmp(self - other, 0)
    def __pow__(self, other):
        if isinstance(other, (int, long)):
            if other < 0:
                return Rational(self._d ** -other, self._n ** -other)
            else:
                return Rational(self._n ** other, self._d ** other)
        else:
                return float(self) ** other
    def __rpow__(self, other):
        return other ** float(self)
    def decimal(self):
        """Return a Decimal approximation of self in the current
context."""
        return _decimal.Decimal(self._n) / _decimal.Decimal(self._d)
    def round(self, denominator):
        """Return self rounded to nearest multiple of 1/denominator."""
        int_part, frac_part = divmod(self * denominator, 1)
        round_direction = cmp(frac_part * 2, 1)
        if round_direction == 0:
           numerator = int_part + (int_part & 1) # round to even
        elif round_direction < 0:
           numerator = int_part
        else:
           numerator = int_part + 1
        return Rational(numerator, denominator)
    @staticmethod
    def from_exact_float(x):
        """Returns the exact Rational equivalent of x."""
        mantissa, exponent = _math.frexp(x)
        mantissa = int(mantissa * 2 ** 53)
        exponent -= 53
        if exponent < 0:
            return Rational(mantissa, 2 ** (-exponent))
        else:
            return Rational(mantissa * 2 ** exponent)
    @staticmethod
    def from_exact_decimal(x):
        """Returns the exact Rational equivalent of x."""
        sign, mantissa, exponent = x.as_tuple()
        sign = (1, -1)[sign]
        mantissa = sign * reduce(lambda a, b: 10 * a + b, mantissa)
        if exponent < 0:
            return Rational(mantissa, 10 ** (-exponent))
        else:
            return Rational(mantissa * 10 ** exponent)
    @staticmethod
    def approx_smallest_denominator(x, tolerance):
        """
        Returns a Rational approximation of x.
        Minimizes the denominator given a constraint on the error.

        x = the float or Decimal value to convert
        tolerance = maximum absolute error allowed,
                    must be of the same type as x
        """
        tolerance = abs(tolerance)
        n = 1
        while True:
            m = int(round(x * n))
            result = Rational(m, n)
            if abs(result - x) < tolerance:
                return result
            n += 1
    @staticmethod
    def approx_smallest_error(x, maxDenominator):
        """
        Returns a Rational approximation of x.
        Minimizes the error given a constraint on the denominator.

        x = the float or Decimal value to convert
        maxDenominator = maximum denominator allowed
        """
        result = None
        minError = x
        for n in xrange(1, maxDenominator + 1):
            m = int(round(x * n))
            r = Rational(m, n)
            error = abs(r - x)
            if error == 0:
                return r
            elif error < minError:
                result = r
                minError = error
        return result

def divide(x, y):
    """Same as x/y, but returns a Rational if both are ints."""
    if isinstance(x, (int, long)) and isinstance(y, (int, long)):
        return Rational(x, y)
    else:
        return x / y