Index: Objects/longobject.c =================================================================== --- Objects/longobject.c (revision 63454) +++ Objects/longobject.c (working copy) @@ -8,6 +8,7 @@ #include "longintrepr.h" #include "formatter_string.h" +#include #include /* For long multiplication, use the O(N**2) school algorithm unless @@ -33,6 +34,7 @@ /* Forward */ static PyLongObject *long_normalize(PyLongObject *); +static PyObject *long_lshift(PyObject *, PyObject *); static PyLongObject *mul1(PyLongObject *, wdigit); static PyLongObject *muladd1(PyLongObject *, wdigit, wdigit); static PyLongObject *divrem1(PyLongObject *, digit, digit *); @@ -2656,47 +2658,149 @@ static PyObject * long_true_divide(PyObject *v, PyObject *w) { - PyLongObject *a, *b; + PyLongObject *a, *b, *q, *r, *temp, *lshift; double ad, bd; - int failed, aexp = -1, bexp = -1; + const double multiplier = (double)(1L << PyLong_SHIFT); + int failed, positive, shift, increment, i; + digit qlow; + Py_ssize_t e; + /* Method: compute the value q = int(a/b*2**-e) for a suitable integer + * e; adjust q and e to deal with correct rounding; then return + * q*2**e. Roughly speaking, we choose e so that q has the right + * number of bits for a double, plus one or two extra bits to + * faciliate rounding decisions. More precisely, the choices for e + * below ensure that either: + * + * q = int(2**-e * a/b) has DBL_MANT_DIG+1 or DBL_MANT_DIG+2 bits, + * and e >= DBL_MIN_EXP-DBL_MANT_DIG-1, or + * + * q = int(2**-e * a/b) has <= DBL_MANT_DIG+1 bits, and e == + * DBL_MIN_EXP-DBL_MANT_DIG-1. + * + * The restriction e >= DBL_MIN_EXP-DBL_MANT_DIG-1 is used to avoid + * double rounding, leading to incorrect results, when the result is + * an IEEE 754 subnormal. It should be harmless on non IEEE + * platforms. + */ + CONVERT_BINOP(v, w, &a, &b); - ad = _PyLong_AsScaledDouble((PyObject *)a, &aexp); - bd = _PyLong_AsScaledDouble((PyObject *)b, &bexp); - failed = (ad == -1.0 || bd == -1.0) && PyErr_Occurred(); + + if (Py_SIZE(b) == 0) { + Py_DECREF(a); + Py_DECREF(b); + PyErr_SetString(PyExc_ZeroDivisionError, + "long division or modulo by zero"); + return NULL; + } + + /* fast path for a and b small */ + if (ABS(Py_SIZE(a)) <= DBL_MANT_DIG/PyLong_SHIFT && + ABS(Py_SIZE(b)) <= DBL_MANT_DIG/PyLong_SHIFT) { + ad = PyLong_AsDouble((PyObject *)a); + bd = PyLong_AsDouble((PyObject *)b); + failed = (ad == -1.0 || bd == -1.0) && PyErr_Occurred(); + Py_DECREF(a); + Py_DECREF(b); + if (failed) + return NULL; + return PyFloat_FromDouble(ad/bd); + } + + /* determine whether result should be positively signed */ + positive = (Py_SIZE(a) >= 0) == (Py_SIZE(b) >= 0); + + /* a == 0 -> +-0.0 */ + if (Py_SIZE(a) == 0) { + Py_DECREF(a); + Py_DECREF(b); + return PyFloat_FromDouble(positive ? 0.0 : -0.0); + } + + /* choose value for e */ + e = _PyLong_NumBits((PyObject *)a) - _PyLong_NumBits((PyObject *)b); + if (e < DBL_MIN_EXP) + e = DBL_MIN_EXP - DBL_MANT_DIG - 1; + else + e -= DBL_MANT_DIG + 1; + + /* find q = int(a/b*2**-e) using integer arithmetic */ + if (e <= 0) { + /* a <<= -e */ + lshift = (PyLongObject *)PyLong_FromSsize_t(-e); + if (lshift == NULL) { + Py_DECREF(a); + Py_DECREF(b); + return NULL; + } + temp = (PyLongObject *) + long_lshift((PyObject *)a, (PyObject *)lshift); + Py_DECREF(lshift); + Py_DECREF(a); + if (temp == NULL) { + Py_DECREF(b); + return NULL; + } + a = temp; + } + else { + /* b <<= e */ + lshift = (PyLongObject *)PyLong_FromSsize_t(e); + if (lshift == NULL) { + Py_DECREF(a); + Py_DECREF(b); + return NULL; + } + temp = (PyLongObject *) + long_lshift((PyObject *)b, (PyObject *)lshift); + Py_DECREF(lshift); + Py_DECREF(b); + if (temp == NULL) { + Py_DECREF(a); + return NULL; + } + b = temp; + } + failed = long_divrem(a, b, &q, &r); Py_DECREF(a); Py_DECREF(b); if (failed) return NULL; - /* 'aexp' and 'bexp' were initialized to -1 to silence gcc-4.0.x, - but should really be set correctly after sucessful calls to - _PyLong_AsScaledDouble() */ - assert(aexp >= 0 && bexp >= 0); - if (bd == 0.0) { - PyErr_SetString(PyExc_ZeroDivisionError, - "long division or modulo by zero"); - return NULL; + /* q = 0 -> underflow to +-0.0 */ + if (Py_SIZE(q) == 0) { + Py_DECREF(q); + Py_DECREF(r); + return PyFloat_FromDouble(positive ? 0.0 : -0.0); } - /* True value is very close to ad/bd * 2**(PyLong_SHIFT*(aexp-bexp)) */ - ad /= bd; /* overflow/underflow impossible here */ - aexp -= bexp; - if (aexp > INT_MAX / PyLong_SHIFT) - goto overflow; - else if (aexp < -(INT_MAX / PyLong_SHIFT)) - return PyFloat_FromDouble(0.0); /* underflow to 0 */ + /* shift and round q, putting result in ad */ + qlow = q->ob_digit[0]; + if (_PyLong_NumBits((PyObject *)q) == DBL_MANT_DIG+2) { + shift = 2; + increment = (qlow & 2) && ((qlow & 5) || Py_SIZE(r)); + } + else { + shift = 1; + increment = (qlow & 1) && ((qlow & 2) || Py_SIZE(r)); + } + ad = (double)((qlow >> shift) + increment); + bd = (double)(1L << (PyLong_SHIFT-shift)); + for (i=1; iob_digit[i]; + bd *= multiplier; + } + Py_DECREF(q); + Py_DECREF(r); + e += shift; errno = 0; - ad = ldexp(ad, aexp * PyLong_SHIFT); - if (Py_OVERFLOWED(ad)) /* ignore underflow to 0.0 */ - goto overflow; - return PyFloat_FromDouble(ad); - -overflow: - PyErr_SetString(PyExc_OverflowError, - "long/long too large for a float"); - return NULL; - + assert(e >= DBL_MIN_EXP - DBL_MANT_DIG); + if (Py_OVERFLOWED(ad)) { + PyErr_SetString(PyExc_OverflowError, + "long/long too large for a float"); + return NULL; + } + return PyFloat_FromDouble(positive ? ad : -ad); } static PyObject * Index: Lib/test/test_long_future.py =================================================================== --- Lib/test/test_long_future.py (revision 63454) +++ Lib/test/test_long_future.py (working copy) @@ -5,6 +5,7 @@ import unittest from test.test_support import run_unittest +from random import randrange class TrueDivisionTests(unittest.TestCase): def test(self): @@ -25,6 +26,73 @@ self.assertEqual(huge / (huge << 1), 0.5) self.assertEqual((1000000 * huge) / huge, 1000000) + # issue #1811: check result is correctly rounded + for i in xrange(100): + self.assertEqual(10**(i+1) / 10**i, 10.0); + + # the next few tests assume IEEE 754 floats + if float.__getformat__("double").startswith("IEEE"): + # 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 xrange(1000): + a = randrange(1, M) + b = randrange(a, 2*a+1) + if randrange(2): + a = -a + if randrange(2): + b = -b + ans = a/b + q = int(ans*2**53) + err = abs(a*2**54-b*q*2) + self.assert_(err <= abs(b) - q%2) + + # one particularly bad case: in Python 2.5 the result of + # the following division is 0.99999999697597697, which is + # out by almost 3.5 ulps. + m = 295147931372582273023 + n = 295147932265116303360 + self.assertEqual(m/n, 0.9999999969759773) + + # test round-half-to-even behaviour, normal result + for m in [1, 2, 7, 17, 12345, 7**100, + -1, -2, -5, -23, -67891, -41**50]: + self.assertEqual(2**53*m/m, 9007199254740992.0) + self.assertEqual((2**53+1)*m/m, 9007199254740992.0) + self.assertEqual((2**53+2)*m/m, 9007199254740994.0) + self.assertEqual((2**53+3)*m/m, 9007199254740996.0) + self.assertEqual((2**53+4)*m/m, 9007199254740996.0) + self.assertEqual((2**53+5)*m/m, 9007199254740996.0) + self.assertEqual(-2**53*m/m, -9007199254740992.0) + self.assertEqual(-(2**53+1)*m/m, -9007199254740992.0) + self.assertEqual(-(2**53+2)*m/m, -9007199254740994.0) + self.assertEqual(-(2**53+3)*m/m, -9007199254740996.0) + self.assertEqual(-(2**53+4)*m/m, -9007199254740996.0) + self.assertEqual(-(2**53+5)*m/m, -9007199254740996.0) + + # round-half-to-even, subnormal result + self.assertEqual((-6)/2**1076, -9.8813129168249309e-324) + self.assertEqual((-5)/2**1076, -4.9406564584124654e-324) + self.assertEqual((-4)/2**1076, -4.9406564584124654e-324) + self.assertEqual((-3)/2**1076, -4.9406564584124654e-324) + self.assertEqual((-2)/2**1076, -0.0) + self.assertEqual((-1)/2**1076, -0.0) + self.assertEqual(0/2**1076, 0.0) + self.assertEqual(1/2**1076, 0.0) + self.assertEqual(2/2**1076, 0.0) + self.assertEqual(3/2**1076, 4.9406564584124654e-324) + self.assertEqual(4/2**1076, 4.9406564584124654e-324) + self.assertEqual(5/2**1076, 4.9406564584124654e-324) + self.assertEqual(6/2**1076, 9.8813129168249309e-324) + self.assertEqual(7/2**1076, 9.8813129168249309e-324) + self.assertEqual(8/2**1076, 9.8813129168249309e-324) + self.assertEqual(9/2**1076, 9.8813129168249309e-324) + self.assertEqual(10/2**1076, 9.8813129168249309e-324) + + + namespace = {'huge': huge, 'mhuge': mhuge} for overflow in ["float(huge)", "float(mhuge)",