468,733 Members | 1,474 Online

# precision problems in base conversion of rational numbers

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'm pretty close but I know I am not handling the precision quite
right. Nothing other than understanding hangs on it, but, that's
enough :-)

To do all this I'm using the decimal module (for the first time) and
I've been setting the precision with

decimal.getcontext().prec = max(getcontext().prec,
x * self.precision )

This is in my class's __init__ method before I convert every number to
Decimal type and self.precision is at this point an int passed.

The first term is to be sure not to reduce precision anywhere else. In
the last term I'd started off with x = 1, and that works well enough
for "small" cases (i.e. cases where I demanded a relatively low degree
of precision in the output).

I've no idea how to work out what x should be in general. (I realize
the answer may be a function of my choice of algorithm. If it is
needed, I'm happy to extract the relevant chunks of code in a
follow-up -- this is long already.)

For getcontext.prec = 80 (i.e x = 1) when I work out 345 / 756 in base
17 to 80 places (i.e. self.precision = 80) I get:
print Rational_in_base_n(345, 756, 17, 80)

0.7CF0CA7CF0CA7CF0CA7CF0CA7CF0CA7CF0CA7CF0CA7CF0CA 7CF0CA7CF0CA7CF0D5C1G603999179EB

(Rational_in_base_n(numerator, denominator, base, precision) is the
rational specific subclass. When I convert the result back to decimal
notation by hand it agrees with the correct answer to as many places
as I cared to check.)

I've discovered empirically that I have to set getcontext.prec = 99 or
greater (i.e. x >= 1.2375) to get
0.7CF0CA7CF0CA7CF0CA7CF0CA7CF0CA7CF0CA7CF0CA7CF0CA 7CF0CA7CF0CA7CF0CA7CF0CA7CF0CA7C

(I do feel safe in assuming that is the right answer. :-)

How would I go about working out what degree of precision for the
decimal module's context is needed to get n 'digits' of precision in
an expression of base m in general? (I've no formal training in Comp
Sci, nor in the relevant branches of mathematics.)

Thanks and best,

Brian vdB

Jul 21 '05 #1
2 2589

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'm pretty close but I know I am not handling the precision quite
right. Nothing other than understanding hangs on it, but, that's
enough :-)

To do all this I'm using the decimal module (for the first time) and
I've been setting the precision with

decimal.getcontext().prec = max(getcontext().prec,
x * self.precision )

This is in my class's __init__ method before I convert every number to
Decimal type and self.precision is at this point an int passed.

The first term is to be sure not to reduce precision anywhere else. In
the last term I'd started off with x = 1, and that works well enough
for "small" cases (i.e. cases where I demanded a relatively low degree
of precision in the output).

I've no idea how to work out what x should be in general. (I realize
the answer may be a function of my choice of algorithm. If it is
needed, I'm happy to extract the relevant chunks of code in a
follow-up -- this is long already.)

For getcontext.prec = 80 (i.e x = 1) when I work out 345 / 756 in base
17 to 80 places (i.e. self.precision = 80) I get:
>>> print Rational_in_base_n(345, 756, 17, 80) 0.7CF0CA7CF0CA7CF0CA7CF0CA7CF0CA7CF0CA7CF0CA7CF0CA 7CF0CA7CF0CA7CF0D5C1G603999179EB

(Rational_in_base_n(numerator, denominator, base, precision) is the
rational specific subclass. When I convert the result back to decimal
notation by hand it agrees with the correct answer to as many places
as I cared to check.)

I've discovered empirically that I have to set getcontext.prec = 99 or
greater (i.e. x >= 1.2375) to get
The value you want for x is log(17)/log(10) =
1.2304489213782739285401698943283

Multiplied by 80 gives you 98.435913710261914283213591546267, but you
can't
have a fraction of a digit, so round up to 99 (giving x=1.235).
0.7CF0CA7CF0CA7CF0CA7CF0CA7CF0CA7CF0CA7CF0CA7CF0CA 7CF0CA7CF0CA7CF0CA7CF0CA7CF0CA7C

(I do feel safe in assuming that is the right answer. :-)
It's not. The last digit should be "D".

This is how the GMPY module handles it:
from gmpy import * # Create an unlimited precision rational.
r = mpq(345,756) # Convert rational to floating point with 300 bits precision.
# GMPY default precision is based on the nature of the rational
# but can be specified. Note that prcision is in bits, not digits.
f = mpf(r,300)
print f 0.456349206349206349206349206349206349206349206349 20634920634920634920634920634920634920634920634921
# Convert the floating point to base 17.
# Note that it has been correctly rounded up, so last digit is 'd'.
print fdigits(f,17)

7.cf0ca7cf0ca7cf0ca7cf0ca7cf0ca7cf0ca7cf0ca7cf0ca7 cf0ca7cf0ca7cf0ca7cf0ca7cf0ca7d@-1

How would I go about working out what degree of precision for the
decimal module's context is needed to get n 'digits' of precision in
an expression of base m in general? (I've no formal training in Comp
Sci, nor in the relevant branches of mathematics.)

Thanks and best,

Brian vdB

Jul 21 '05 #2
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.7CF0CA7CF0CA7CF0CA7CF0CA7CF0CA7CF0CA7CF0CA7CF0CA 7CF0CA7CF0CA7CF0CA7CF0CA7CF0CA7C

(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
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
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 __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

Jul 21 '05 #3

### This discussion thread is closed

Replies have been disabled for this discussion.