Index: Objects/longobject.c =================================================================== --- Objects/longobject.c (revision 67665) +++ Objects/longobject.c (working copy) @@ -6,6 +6,7 @@ #include "Python.h" #include "longintrepr.h" +#include "float.h" #include @@ -727,33 +728,123 @@ #undef NBITS_WANTED } -/* Get a C double from a long int object. */ +/* Get a C double from a long int object. Rounds to the nearest double, + using the round-half-to-even rule in the case of a tie. */ double PyLong_AsDouble(PyObject *vv) { - int e = -1; + PyLongObject *v = (PyLongObject *)vv; + Py_ssize_t rnd_digit, rnd_bit, m, n, i; + digit lsb, *d, parity_bit, rounding_bit; double x; if (vv == NULL || !PyLong_Check(vv)) { PyErr_BadInternalCall(); - return -1; + return -1.0; } - x = _PyLong_AsScaledDouble(vv, &e); - if (x == -1.0 && PyErr_Occurred()) - return -1.0; - /* 'e' initialized to -1 to silence gcc-4.0.x, but it should be - set correctly after a successful _PyLong_AsScaledDouble() call */ - assert(e >= 0); - if (e > INT_MAX / PyLong_SHIFT) + + n = (Py_ssize_t)_PyLong_NumBits(vv); + if (n > DBL_MAX_EXP || n < 0) goto overflow; + + /* Terminology: for the purposes of the comments and code below, the + *rounding bit* is the first bit of v that will be rounded away; in + other words, it's bit n - DBL_MANT_DIG - 1, where n is the total + number of bits in v. This is the bit that usually determines + whether we need to round up or down (for round-half-up it would be + the only bit we'd need to look at; but we want to do + round-half-even.). + + The *parity bit* is the bit immediately above (i.e. more + significant than) the rounding bit. For round-half-to-even, this + is the bit that needs to be examined in the case of a tie: if the + parity bit is set then we round up; otherwise, down. + + Strategy: write d[0] ... d[m-1] for the digits of v (least to most + significant). Let rnd_digit be the index of the digit containing + the rounding bit. Then we initialize a double x to be the integer + given by digits d[rnd_digit:m], but with the bottom + (1+rnd_bit%DBL_MANT_DIG) bits of d[rnd_digit] set to zero. + + In other words, x will have value given by the top DBL_MANT_DIG + bits of v, multiplied by lsb, where lsb is the power of 2 given by: + + lsb = 1 << (1 + rnd_bit % DBL_MANT_DIG) + + Now if we're rounding down, the value we want to return is + + x * 2**(PyLong_SHIFT * rnd_digit). + + and if we're rounding up, it's + + (x + lsb) * 2**(PyLong_SHIFT * rnd_digit). + + where lsb is 1 << (1+rnd_bit%DBL_MANT_DIG); that is, it's the value + of the least significant bit that will appear in the mantissa. + */ + + /* d[0], ..., d[m-1] are the digits of abs(v) */ + m = ABS(Py_SIZE(v)); + d = v->ob_digit; + + if (n <= DBL_MANT_DIG) { + /* no rounding to do; simply accumulate and return */ + x = 0.0; + for (i = m-1; i >= 0; i--) + x = x*PyLong_BASE + d[i]; + return Py_SIZE(v) < 0 ? -x : x; + } + + /* rounding necessary; identify rounding bit position */ + rnd_bit = n - DBL_MANT_DIG - 1; + rnd_digit = rnd_bit/PyLong_SHIFT; + lsb = 1 << (1 + rnd_bit % PyLong_SHIFT); + + /* accumulate top DBL_MANT_DIG bits of v. Note that + all floating-point operations here are *exact*; x never has + more than DBL_MANT_DIG significant digits. */ + x = 0.0; + for (i = m-1; i > rnd_digit; i--) + x = x*PyLong_BASE + d[i]; + x = x*PyLong_BASE + (d[rnd_digit] & (PyLong_BASE - lsb)); + + /* get rounding bit, parity bit */ + rounding_bit = d[rnd_digit] & lsb/2; + if (lsb == PyLong_BASE) + parity_bit = d[rnd_digit+1] & 1; + else + parity_bit = d[rnd_digit] & lsb; + + /* round up if necessary: to implement the round-half-to-even rule, we + round up if the rounding bit is set, and at least one of the + following is true: + + (1) any of the bits of d[rnd_digit] below the rounding + bit is set, or + (2) the bit above the rounding bit (the parity bit) is set, or + (3) any one of the digits d[i], i < rnd_digit is nonzero. + */ + if (rounding_bit) { + if (d[rnd_digit] & (lsb/2-1)) + x += lsb; + else if (parity_bit) + x += lsb; + else + for (i = rnd_digit-1; i >= 0; i--) + if (d[i]) { + x += lsb; + break; + } + } + /* now x has the correct bits; all that remains is to shift */ errno = 0; - x = ldexp(x, e * PyLong_SHIFT); + x = ldexp(x, rnd_digit * PyLong_SHIFT); if (Py_OVERFLOWED(x)) goto overflow; - return x; - -overflow: + /* put the correct sign in, and return */ + return Py_SIZE(v) < 0 ? -x : x; + overflow: PyErr_SetString(PyExc_OverflowError, "long int too large to convert to float"); return -1.0; Index: Lib/test/test_long.py =================================================================== --- Lib/test/test_long.py (revision 67665) +++ Lib/test/test_long.py (working copy) @@ -615,6 +615,45 @@ else: self.assertRaises(TypeError, pow,longx, longy, long(z)) + def test_float_conversion(self): + # the tests are IEEE-754 specific; skip them on non IEEE machines + if not float.__getformat__("double").startswith("IEEE"): + return + + exact_values = [-2L, -1L, 0L, 1L, 2L, + long(2**53-3), + long(2**53-2), + long(2**53-1), + long(2**53), + long(2**53+2), + long(2**54-4), + long(2**54-2), + long(2**54), + long(2**54+4)] + for x in exact_values: + self.assertEqual(long(float(x)), x) + self.assertEqual(long(float(-x)), -x) + + # test round-half-even + for x, y in [(1, 0), (2, 2), (3, 4), (4, 4), (5, 4), (6, 6), (7, 8)]: + for p in xrange(15): + self.assertEqual(long(float(2L**p*(2**53+x))), 2L**p*(2**53+y)) + + for x, y in [(0, 0), (1, 0), (2, 0), (3, 4), (4, 4), (5, 4), (6, 8), + (7, 8), (8, 8), (9, 8), (10, 8), (11, 12), (12, 12), + (13, 12), (14, 16), (15, 16)]: + for p in xrange(15): + self.assertEqual(long(float(2L**p*(2**54+x))), 2L**p*(2**54+y)) + + for p in xrange(100): + x = long(2**p * (2**53 + 1) + 1) + y = long(2**p * (2**53+ 2)) + self.assertEqual(long(float(x)), y) + + x = long(2**p * (2**53 + 1)) + y = long(2**p * 2**53) + self.assertEqual(long(float(x)), y) + def test_float_overflow(self): import math