Issue #1811: Improve accuracy and consistency of true division for integers.

Lib/test/test_long_future.py

@@ -3,9 +3,52 @@ from __future__ import division
 # test_long.py instead.  In the meantime, it's too obscure to try to
 # trick just part of test_long into using future division.
 
+import sys
+import random
 import unittest
 from test.test_support import run_unittest
 
+# decorator for skipping tests on non-IEEE 754 platforms
+requires_IEEE_754 = unittest.skipUnless(
+    float.__getformat__("double").startswith("IEEE"),
+    "test requires IEEE 754 doubles")
+
+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 TrueDivisionTests(unittest.TestCase):
     def test(self):
         huge = 1L << 40000
@@ -45,6 +88,131 @@ class TrueDivisionTests(unittest.TestCase):
             with 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."""
+
+        a, b = long(a), long(b)
+
+        # 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_main():
     run_unittest(TrueDivisionTests)

Misc/NEWS

@@ -12,6 +12,11 @@ What's New in Python 2.7 alpha 2?
 Core and Builtins
 -----------------
 
+- Issue #1811: improve accuracy and cross-platform consistency for
+  true division of integers: the result of a/b is now correctly
+  rounded for ints a and b (at least on IEEE 754 platforms), and in
+  particular does not depend on the internal representation of a long.
+
 - Issue #6108: unicode(exception) and str(exception) should return the same
   message when only __str__ (and not __unicode__) is overridden in the subclass.
 

Objects/intobject.c

@@ -645,16 +645,35 @@ int_classic_div(PyIntObject *x, PyIntObject *y)
 }
 
 static PyObject *
-int_true_divide(PyObject *v, PyObject *w)
+int_true_divide(PyIntObject *x, PyIntObject *y)
 {
+	long xi, yi;
 	/* If they aren't both ints, give someone else a chance.  In
 	   particular, this lets int/long get handled by longs, which
 	   underflows to 0 gracefully if the long is too big to convert
 	   to float. */
-	if (PyInt_Check(v) && PyInt_Check(w))
-		return PyFloat_Type.tp_as_number->nb_true_divide(v, w);
-	Py_INCREF(Py_NotImplemented);
-	return Py_NotImplemented;
+	CONVERT_TO_LONG(x, xi);
+	CONVERT_TO_LONG(y, yi);
+	if (yi == 0) {
+		PyErr_SetString(PyExc_ZeroDivisionError,
+				"division by zero");
+		return NULL;
+	}
+	if (xi == 0)
+		return PyFloat_FromDouble(yi < 0 ? -0.0 : 0.0);
+
+#define WIDTH_OF_ULONG (CHAR_BIT*SIZEOF_LONG)
+#if DBL_MANT_DIG < WIDTH_OF_ULONG
+	if ((xi >= 0 ? 0UL + xi : 0UL - xi) >> DBL_MANT_DIG ||
+	    (yi >= 0 ? 0UL + yi : 0UL - yi) >> DBL_MANT_DIG)
+		/* Large x or y.  Use long integer arithmetic. */
+		return PyLong_Type.tp_as_number->nb_true_divide(
+			(PyObject *)x, (PyObject *)y);
+	else
+#endif
+		/* Both ints can be exactly represented as doubles.  Do a
+		   floating-point division. */
+		return PyFloat_FromDouble((double)xi / (double)yi);
 }
 
 static PyObject *
@@ -1355,7 +1374,7 @@ static PyNumberMethods int_as_number = {
 	0,			/*nb_inplace_xor*/
 	0,			/*nb_inplace_or*/
 	(binaryfunc)int_div,	/* nb_floor_divide */
-	int_true_divide,	/* nb_true_divide */
+	(binaryfunc)int_true_divide, /* nb_true_divide */
 	0,			/* nb_inplace_floor_divide */
 	0,			/* nb_inplace_true_divide */
 	(unaryfunc)int_int,	/* nb_index */