Merged revisions 77062 via svnmerge from
svn+ssh://pythondev@svn.python.org/python/trunk
........
r77062 | mark.dickinson | 2009-12-27 14:55:57 +0000 (Sun, 27 Dec 2009) | 2 lines
Issue #1811: Improve accuracy and consistency of true division for integers.
........
diff --git a/Lib/test/test_long.py b/Lib/test/test_long.py
index a7dbb5c..98221e6 100644
--- a/Lib/test/test_long.py
+++ b/Lib/test/test_long.py
@@ -14,6 +14,11 @@
def __str__(self):
return self.format % self.args
+# decorator for skipping tests on non-IEEE 754 platforms
+requires_IEEE_754 = unittest.skipUnless(
+ float.__getformat__("double").startswith("IEEE"),
+ "test requires IEEE 754 doubles")
+
# SHIFT should match the value in longintrepr.h for best testing.
SHIFT = sys.int_info.bits_per_digit
BASE = 2 ** SHIFT
@@ -35,6 +40,43 @@
# add complements & negations
special += [~x for x in special] + [-x for x in special]
+DBL_MAX = sys.float_info.max
+DBL_MAX_EXP = sys.float_info.max_exp
+DBL_MIN_EXP = sys.float_info.min_exp
+DBL_MANT_DIG = sys.float_info.mant_dig
+DBL_MIN_OVERFLOW = 2**DBL_MAX_EXP - 2**(DBL_MAX_EXP - DBL_MANT_DIG - 1)
+
+# pure Python version of correctly-rounded true division
+def truediv(a, b):
+ """Correctly-rounded true division for integers."""
+ negative = a^b < 0
+ a, b = abs(a), abs(b)
+
+ # exceptions: division by zero, overflow
+ if not b:
+ raise ZeroDivisionError("division by zero")
+ if a >= DBL_MIN_OVERFLOW * b:
+ raise OverflowError("int/int too large to represent as a float")
+
+ # find integer d satisfying 2**(d - 1) <= a/b < 2**d
+ d = a.bit_length() - b.bit_length()
+ if d >= 0 and a >= 2**d * b or d < 0 and a * 2**-d >= b:
+ d += 1
+
+ # compute 2**-exp * a / b for suitable exp
+ exp = max(d, DBL_MIN_EXP) - DBL_MANT_DIG
+ a, b = a << max(-exp, 0), b << max(exp, 0)
+ q, r = divmod(a, b)
+
+ # round-half-to-even: fractional part is r/b, which is > 0.5 iff
+ # 2*r > b, and == 0.5 iff 2*r == b.
+ if 2*r > b or 2*r == b and q % 2 == 1:
+ q += 1
+
+ result = float(q) * 2.**exp
+ return -result if negative else result
+
+
class LongTest(unittest.TestCase):
# Get quasi-random long consisting of ndigits digits (in base BASE).
@@ -306,10 +348,6 @@
@unittest.skipUnless(float.__getformat__("double").startswith("IEEE"),
"test requires IEEE 754 doubles")
def test_float_conversion(self):
- import sys
- DBL_MAX = sys.float_info.max
- DBL_MAX_EXP = sys.float_info.max_exp
- DBL_MANT_DIG = sys.float_info.mant_dig
exact_values = [0, 1, 2,
2**53-3,
@@ -614,6 +652,128 @@
for zero in ["huge / 0", "mhuge / 0"]:
self.assertRaises(ZeroDivisionError, eval, zero, namespace)
+ def check_truediv(self, a, b, skip_small=True):
+ """Verify that the result of a/b is correctly rounded, by
+ comparing it with a pure Python implementation of correctly
+ rounded division. b should be nonzero."""
+
+ # skip check for small a and b: in this case, the current
+ # implementation converts the arguments to float directly and
+ # then applies a float division. This can give doubly-rounded
+ # results on x87-using machines (particularly 32-bit Linux).
+ if skip_small and max(abs(a), abs(b)) < 2**DBL_MANT_DIG:
+ return
+
+ try:
+ # use repr so that we can distinguish between -0.0 and 0.0
+ expected = repr(truediv(a, b))
+ except OverflowError:
+ expected = 'overflow'
+ except ZeroDivisionError:
+ expected = 'zerodivision'
+
+ try:
+ got = repr(a / b)
+ except OverflowError:
+ got = 'overflow'
+ except ZeroDivisionError:
+ got = 'zerodivision'
+
+ if expected != got:
+ self.fail("Incorrectly rounded division {}/{}: expected {!r}, "
+ "got {!r}.".format(a, b, expected, got))
+
+ @requires_IEEE_754
+ def test_correctly_rounded_true_division(self):
+ # more stringent tests than those above, checking that the
+ # result of true division of ints is always correctly rounded.
+ # This test should probably be considered CPython-specific.
+
+ # Exercise all the code paths not involving Gb-sized ints.
+ # ... divisions involving zero
+ self.check_truediv(123, 0)
+ self.check_truediv(-456, 0)
+ self.check_truediv(0, 3)
+ self.check_truediv(0, -3)
+ self.check_truediv(0, 0)
+ # ... overflow or underflow by large margin
+ self.check_truediv(671 * 12345 * 2**DBL_MAX_EXP, 12345)
+ self.check_truediv(12345, 345678 * 2**(DBL_MANT_DIG - DBL_MIN_EXP))
+ # ... a much larger or smaller than b
+ self.check_truediv(12345*2**100, 98765)
+ self.check_truediv(12345*2**30, 98765*7**81)
+ # ... a / b near a boundary: one of 1, 2**DBL_MANT_DIG, 2**DBL_MIN_EXP,
+ # 2**DBL_MAX_EXP, 2**(DBL_MIN_EXP-DBL_MANT_DIG)
+ bases = (0, DBL_MANT_DIG, DBL_MIN_EXP,
+ DBL_MAX_EXP, DBL_MIN_EXP - DBL_MANT_DIG)
+ for base in bases:
+ for exp in range(base - 15, base + 15):
+ self.check_truediv(75312*2**max(exp, 0), 69187*2**max(-exp, 0))
+ self.check_truediv(69187*2**max(exp, 0), 75312*2**max(-exp, 0))
+
+ # overflow corner case
+ for m in [1, 2, 7, 17, 12345, 7**100,
+ -1, -2, -5, -23, -67891, -41**50]:
+ for n in range(-10, 10):
+ self.check_truediv(m*DBL_MIN_OVERFLOW + n, m)
+ self.check_truediv(m*DBL_MIN_OVERFLOW + n, -m)
+
+ # check detection of inexactness in shifting stage
+ for n in range(250):
+ # (2**DBL_MANT_DIG+1)/(2**DBL_MANT_DIG) lies halfway
+ # between two representable floats, and would usually be
+ # rounded down under round-half-to-even. The tiniest of
+ # additions to the numerator should cause it to be rounded
+ # up instead.
+ self.check_truediv((2**DBL_MANT_DIG + 1)*12345*2**200 + 2**n,
+ 2**DBL_MANT_DIG*12345)
+
+ # 1/2731 is one of the smallest division cases that's subject
+ # to double rounding on IEEE 754 machines working internally with
+ # 64-bit precision. On such machines, the next check would fail,
+ # were it not explicitly skipped in check_truediv.
+ self.check_truediv(1, 2731)
+
+ # a particularly bad case for the old algorithm: gives an
+ # error of close to 3.5 ulps.
+ self.check_truediv(295147931372582273023, 295147932265116303360)
+ for i in range(1000):
+ self.check_truediv(10**(i+1), 10**i)
+ self.check_truediv(10**i, 10**(i+1))
+
+ # test round-half-to-even behaviour, normal result
+ for m in [1, 2, 4, 7, 8, 16, 17, 32, 12345, 7**100,
+ -1, -2, -5, -23, -67891, -41**50]:
+ for n in range(-10, 10):
+ self.check_truediv(2**DBL_MANT_DIG*m + n, m)
+
+ # test round-half-to-even, subnormal result
+ for n in range(-20, 20):
+ self.check_truediv(n, 2**1076)
+
+ # largeish random divisions: a/b where |a| <= |b| <=
+ # 2*|a|; |ans| is between 0.5 and 1.0, so error should
+ # always be bounded by 2**-54 with equality possible only
+ # if the least significant bit of q=ans*2**53 is zero.
+ for M in [10**10, 10**100, 10**1000]:
+ for i in range(1000):
+ a = random.randrange(1, M)
+ b = random.randrange(a, 2*a+1)
+ self.check_truediv(a, b)
+ self.check_truediv(-a, b)
+ self.check_truediv(a, -b)
+ self.check_truediv(-a, -b)
+
+ # and some (genuinely) random tests
+ for _ in range(10000):
+ a_bits = random.randrange(1000)
+ b_bits = random.randrange(1, 1000)
+ x = random.randrange(2**a_bits)
+ y = random.randrange(1, 2**b_bits)
+ self.check_truediv(x, y)
+ self.check_truediv(x, -y)
+ self.check_truediv(-x, y)
+ self.check_truediv(-x, -y)
def test_small_ints(self):
for i in range(-5, 257):