diff options
| -rw-r--r-- | Lib/test/test_long.py | 168 | ||||
| -rw-r--r-- | Misc/NEWS | 5 | ||||
| -rw-r--r-- | Objects/longobject.c | 280 |
3 files changed, 419 insertions, 34 deletions
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 @@ class Frm(object): 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 @@ del p2 # 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 @@ class LongTest(unittest.TestCase): @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 @@ class LongTest(unittest.TestCase): 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): @@ -12,6 +12,11 @@ What's New in Python 3.2 Alpha 1? 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 an int. + - Issue #6834: replace the implementation for the 'python' and 'pythonw' executables on OSX. diff --git a/Objects/longobject.c b/Objects/longobject.c index 8e4093c..429b635 100644 --- a/Objects/longobject.c +++ b/Objects/longobject.c @@ -3213,47 +3213,267 @@ long_div(PyObject *a, PyObject *b) return (PyObject *)div; } +/* PyLong/PyLong -> float, with correctly rounded result. */ + +#define MANT_DIG_DIGITS (DBL_MANT_DIG / PyLong_SHIFT) +#define MANT_DIG_BITS (DBL_MANT_DIG % PyLong_SHIFT) + static PyObject * -long_true_divide(PyObject *a, PyObject *b) +long_true_divide(PyObject *v, PyObject *w) { - double ad, bd; - int failed, aexp = -1, bexp = -1; + PyLongObject *a, *b, *x; + Py_ssize_t a_size, b_size, shift, extra_bits, diff, x_size, x_bits; + digit mask, low; + int inexact, negate, a_is_small, b_is_small; + double dx, result; - CHECK_BINOP(a, b); - ad = _PyLong_AsScaledDouble((PyObject *)a, &aexp); - bd = _PyLong_AsScaledDouble((PyObject *)b, &bexp); - failed = (ad == -1.0 || bd == -1.0) && PyErr_Occurred(); - 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); + CHECK_BINOP(v, w); + a = (PyLongObject *)v; + b = (PyLongObject *)w; + + /* + Method in a nutshell: + + 0. reduce to case a, b > 0; filter out obvious underflow/overflow + 1. choose a suitable integer 'shift' + 2. use integer arithmetic to compute x = floor(2**-shift*a/b) + 3. adjust x for correct rounding + 4. convert x to a double dx with the same value + 5. return ldexp(dx, shift). + + In more detail: + + 0. For any a, a/0 raises ZeroDivisionError; for nonzero b, 0/b + returns either 0.0 or -0.0, depending on the sign of b. For a and + b both nonzero, ignore signs of a and b, and add the sign back in + at the end. Now write a_bits and b_bits for the bit lengths of a + and b respectively (that is, a_bits = 1 + floor(log_2(a)); likewise + for b). Then + + 2**(a_bits - b_bits - 1) < a/b < 2**(a_bits - b_bits + 1). + + So if a_bits - b_bits > DBL_MAX_EXP then a/b > 2**DBL_MAX_EXP and + so overflows. Similarly, if a_bits - b_bits < DBL_MIN_EXP - + DBL_MANT_DIG - 1 then a/b underflows to 0. With these cases out of + the way, we can assume that + + DBL_MIN_EXP - DBL_MANT_DIG - 1 <= a_bits - b_bits <= DBL_MAX_EXP. + + 1. The integer 'shift' is chosen so that x has the right number of + bits for a double, plus two or three extra bits that will be used + in the rounding decisions. Writing a_bits and b_bits for the + number of significant bits in a and b respectively, a + straightforward formula for shift is: + + shift = a_bits - b_bits - DBL_MANT_DIG - 2 + + This is fine in the usual case, but if a/b is smaller than the + smallest normal float then it can lead to double rounding on an + IEEE 754 platform, giving incorrectly rounded results. So we + adjust the formula slightly. The actual formula used is: + + shift = MAX(a_bits - b_bits, DBL_MIN_EXP) - DBL_MANT_DIG - 2 + + 2. The quantity x is computed by first shifting a (left -shift bits + if shift <= 0, right shift bits if shift > 0) and then dividing by + b. For both the shift and the division, we keep track of whether + the result is inexact, in a flag 'inexact'; this information is + needed at the rounding stage. + + With the choice of shift above, together with our assumption that + a_bits - b_bits >= DBL_MIN_EXP - DBL_MANT_DIG - 1, it follows + that x >= 1. + + 3. Now x * 2**shift <= a/b < (x+1) * 2**shift. We want to replace + this with an exactly representable float of the form + + round(x/2**extra_bits) * 2**(extra_bits+shift). + + For float representability, we need x/2**extra_bits < + 2**DBL_MANT_DIG and extra_bits + shift >= DBL_MIN_EXP - + DBL_MANT_DIG. This translates to the condition: + + extra_bits >= MAX(x_bits, DBL_MIN_EXP - shift) - DBL_MANT_DIG + + To round, we just modify the bottom digit of x in-place; this can + end up giving a digit with value > PyLONG_MASK, but that's not a + problem since digits can hold values up to 2*PyLONG_MASK+1. + + With the original choices for shift above, extra_bits will always + be 2 or 3. Then rounding under the round-half-to-even rule, we + round up iff the most significant of the extra bits is 1, and + either: (a) the computation of x in step 2 had an inexact result, + or (b) at least one other of the extra bits is 1, or (c) the least + significant bit of x (above those to be rounded) is 1. + + 4. Conversion to a double is straightforward; all floating-point + operations involved in the conversion are exact, so there's no + danger of rounding errors. + + 5. Use ldexp(x, shift) to compute x*2**shift, the final result. + The result will always be exactly representable as a double, except + in the case that it overflows. To avoid dependence on the exact + behaviour of ldexp on overflow, we check for overflow before + applying ldexp. The result of ldexp is adjusted for sign before + returning. + */ - if (bd == 0.0) { + /* Reduce to case where a and b are both positive. */ + a_size = ABS(Py_SIZE(a)); + b_size = ABS(Py_SIZE(b)); + negate = (Py_SIZE(a) < 0) ^ (Py_SIZE(b) < 0); + if (b_size == 0) { PyErr_SetString(PyExc_ZeroDivisionError, - "integer division or modulo by zero"); - return NULL; + "division by zero"); + goto error; } - - /* 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) + if (a_size == 0) + goto underflow_or_zero; + + /* Fast path for a and b small (exactly representable in a double). + Relies on floating-point division being correctly rounded; results + may be subject to double rounding on x86 machines that operate with + the x87 FPU set to 64-bit precision. */ + a_is_small = a_size <= MANT_DIG_DIGITS || + (a_size == MANT_DIG_DIGITS+1 && + a->ob_digit[MANT_DIG_DIGITS] >> MANT_DIG_BITS == 0); + b_is_small = b_size <= MANT_DIG_DIGITS || + (b_size == MANT_DIG_DIGITS+1 && + b->ob_digit[MANT_DIG_DIGITS] >> MANT_DIG_BITS == 0); + if (a_is_small && b_is_small) { + double da, db; + da = a->ob_digit[--a_size]; + while (a_size > 0) + da = da * PyLong_BASE + a->ob_digit[--a_size]; + db = b->ob_digit[--b_size]; + while (b_size > 0) + db = db * PyLong_BASE + b->ob_digit[--b_size]; + result = da / db; + goto success; + } + + /* Catch obvious cases of underflow and overflow */ + diff = a_size - b_size; + if (diff > PY_SSIZE_T_MAX/PyLong_SHIFT - 1) + /* Extreme overflow */ goto overflow; - else if (aexp < -(INT_MAX / PyLong_SHIFT)) - return PyFloat_FromDouble(0.0); /* underflow to 0 */ - errno = 0; - ad = ldexp(ad, aexp * PyLong_SHIFT); - if (Py_OVERFLOWED(ad)) /* ignore underflow to 0.0 */ + else if (diff < 1 - PY_SSIZE_T_MAX/PyLong_SHIFT) + /* Extreme underflow */ + goto underflow_or_zero; + /* Next line is now safe from overflowing a Py_ssize_t */ + diff = diff * PyLong_SHIFT + bits_in_digit(a->ob_digit[a_size - 1]) - + bits_in_digit(b->ob_digit[b_size - 1]); + /* Now diff = a_bits - b_bits. */ + if (diff > DBL_MAX_EXP) goto overflow; - return PyFloat_FromDouble(ad); + else if (diff < DBL_MIN_EXP - DBL_MANT_DIG - 1) + goto underflow_or_zero; + + /* Choose value for shift; see comments for step 1 above. */ + shift = MAX(diff, DBL_MIN_EXP) - DBL_MANT_DIG - 2; + + inexact = 0; + + /* x = abs(a * 2**-shift) */ + if (shift <= 0) { + Py_ssize_t i, shift_digits = -shift / PyLong_SHIFT; + digit rem; + /* x = a << -shift */ + if (a_size >= PY_SSIZE_T_MAX - 1 - shift_digits) { + /* In practice, it's probably impossible to end up + here. Both a and b would have to be enormous, + using close to SIZE_T_MAX bytes of memory each. */ + PyErr_SetString(PyExc_OverflowError, + "intermediate overflow during division"); + goto error; + } + x = _PyLong_New(a_size + shift_digits + 1); + if (x == NULL) + goto error; + for (i = 0; i < shift_digits; i++) + x->ob_digit[i] = 0; + rem = v_lshift(x->ob_digit + shift_digits, a->ob_digit, + a_size, -shift % PyLong_SHIFT); + x->ob_digit[a_size + shift_digits] = rem; + } + else { + Py_ssize_t shift_digits = shift / PyLong_SHIFT; + digit rem; + /* x = a >> shift */ + assert(a_size >= shift_digits); + x = _PyLong_New(a_size - shift_digits); + if (x == NULL) + goto error; + rem = v_rshift(x->ob_digit, a->ob_digit + shift_digits, + a_size - shift_digits, shift % PyLong_SHIFT); + /* set inexact if any of the bits shifted out is nonzero */ + if (rem) + inexact = 1; + while (!inexact && shift_digits > 0) + if (a->ob_digit[--shift_digits]) + inexact = 1; + } + long_normalize(x); + x_size = Py_SIZE(x); + + /* x //= b. If the remainder is nonzero, set inexact. We own the only + reference to x, so it's safe to modify it in-place. */ + if (b_size == 1) { + digit rem = inplace_divrem1(x->ob_digit, x->ob_digit, x_size, + b->ob_digit[0]); + long_normalize(x); + if (rem) + inexact = 1; + } + else { + PyLongObject *div, *rem; + div = x_divrem(x, b, &rem); + Py_DECREF(x); + x = div; + if (x == NULL) + goto error; + if (Py_SIZE(rem)) + inexact = 1; + Py_DECREF(rem); + } + x_size = ABS(Py_SIZE(x)); + assert(x_size > 0); /* result of division is never zero */ + x_bits = (x_size-1)*PyLong_SHIFT+bits_in_digit(x->ob_digit[x_size-1]); + + /* The number of extra bits that have to be rounded away. */ + extra_bits = MAX(x_bits, DBL_MIN_EXP - shift) - DBL_MANT_DIG; + assert(extra_bits == 2 || extra_bits == 3); + + /* Round by directly modifying the low digit of x. */ + mask = (digit)1 << (extra_bits - 1); + low = x->ob_digit[0] | inexact; + if (low & mask && low & (3*mask-1)) + low += mask; + x->ob_digit[0] = low & ~(mask-1U); + + /* Convert x to a double dx; the conversion is exact. */ + dx = x->ob_digit[--x_size]; + while (x_size > 0) + dx = dx * PyLong_BASE + x->ob_digit[--x_size]; + Py_DECREF(x); -overflow: + /* Check whether ldexp result will overflow a double. */ + if (shift + x_bits >= DBL_MAX_EXP && + (shift + x_bits > DBL_MAX_EXP || dx == ldexp(1.0, x_bits))) + goto overflow; + result = ldexp(dx, shift); + + success: + return PyFloat_FromDouble(negate ? -result : result); + + underflow_or_zero: + return PyFloat_FromDouble(negate ? -0.0 : 0.0); + + overflow: PyErr_SetString(PyExc_OverflowError, - "int/int too large for a float"); + "integer division result too large for a float"); + error: return NULL; - } static PyObject * |