Index: Objects/longobject.c =================================================================== --- Objects/longobject.c (revision 67706) +++ Objects/longobject.c (working copy) @@ -6,6 +6,7 @@ #include "Python.h" #include "longintrepr.h" +#include "float.h" #include @@ -727,33 +728,192 @@ #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; + const digit *d; + digit lsb, rounding_bit, msd; double x; + int sign; 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) + + /* Here's an overview of the rounding algorithm. For simplicity, + assume that v is positive, and is at least 2**DBL_MANT_DIG. (For + smaller positive v no rounding is necessary, while for negative v + we just ignore the sign, adding it back into the rounded result at + the end.) Write n for the number of bits in v; that is, the unique + integer such that 2**(n-1) <= v < 2**n. + + We introduce some terminology to streamline the explanation: for + the purposes of the comments and code below, the *rounding bit* of + v is the first bit of v that will be rounded away, and the *parity + bit* is the bit immediately above (i.e., more significant than) the + rounding bit. + + In other words, if we number the bits in the usual way with bit 0 + being the least significant, the rounding bit is bit n - + DBL_MANT_DIG - 1. This is the bit that usually determines whether + we need to round up or down (indeed, for round-half-up the rounding + bit would be the *only* bit we'd need to look at; but we want to do + round-half-to-even). The parity bit is bit n - DBL_MANT_DIG. This + is the bit that will need to be examined in the case that v is + exactly halfway between two representable floats: if the parity bit + is set then we round up; otherwise, down. + + Write d[0] ... d[m] for the digits of v (least to most + significant). Let rnd_bit be the index of the rounding bit, and + rnd_digit the index of the digit containing the rounding bit. Then + the bits of d[rnd_digit] look something like: + + rounding bit + | + v + msb -> sssssrttttttttt <- lsb + ^ + | + parity bit + + where each 'r', 's' or 't' is a placeholder for one of the bits, + 's' represents a 'significant bit' that will be included in the + mantissa of the result, 'r' is the rounding bit, and 't' represents + a 'trailing bit' following the rounding bit. Note that if the + rounding bit is at the top of d[rnd_digit] then the parity bit will + be in the next digit up, d[rnd_digit+1]. The number of trailing + bits in d[rnd_digit] is rnd_bit % PyLong_SHIFT. If we set + + lsb = 1 << (rnd_bit % PyLong_SHIFT) + + then d[rnd_digit] & (PyLong_BASE - 2*lsb) selects just the + significant bits of d[rnd_digit]; d[rnd_digit] & (lsb-1) + selects the trailing bits, and d[rnd_digit] & lsb gives the + rounding bit. + + The first step in the algorithm is to initialize a double x so that + its value is the integer given by digits d[rnd_digit:m-1], but with + the rounding bit and trailing bits of d[rnd_digit] set to zero. So + x will have value given by the top DBL_MANT_DIG bits of v, + multiplied by 2*lsb. + + Now if we're rounding down, the value we want to return is simply + + x * 2**(PyLong_SHIFT * rnd_digit). + + and if we're rounding up, it's + + (x + 2*lsb) * 2**(PyLong_SHIFT * rnd_digit). + + Under the round-half-to-even rule, we round up if, and only + if, the rounding bit is set *and* at least one of the + following three conditions is satisfied: + + (1) the parity bit is set, or + (2) at least one of the trailing bits of d[rnd_digit] is set, or + (3) at least one of the digits d[i], 0 <= i < rnd_digit + is nonzero. + + Here, the condition "(2) or (3)" is just the condition that some + bit below the rounding bit is set; in other words, it's the + condition that v does *not* lie exactly halfway between two + representable floats. + + Finally, we have to worry about overflow. If v >= 2**DBL_MAX_EXP, + or equivalently n > DBL_MAX_EXP, then overflow occurs. If v < + 2**DBL_MAX_EXP then we're usually safe, but there's a corner case + to consider: if v is very close to 2**DBL_MAX_EXP then it's + possible that v is rounded up to exactly 2**DBL_MAX_EXP, and then + again overflow occurs. + */ + + if (Py_SIZE(v) == 0) + return 0.0; + m = ABS(Py_SIZE(v)) - 1; + sign = Py_SIZE(v) < 0 ? -1 : 1; + d = v->ob_digit; + msd = d[m]; + + /* fast path for case where abs(v) < 2**DBL_MANT_DIG */ + if (m < DBL_MANT_DIG / PyLong_SHIFT || + (m == DBL_MANT_DIG / PyLong_SHIFT && + msd < (digit)1 << DBL_MANT_DIG%PyLong_SHIFT)) { + x = msd; + for (i = m-1; i >= 0; i--) + x = x*PyLong_BASE + d[i]; + return sign * x; + } + + /* if m is huge then we overflow immediately; if the condition below + holds then n >= m * PyLong_SHIFT + 1 > DBL_MAX_EXP. */ + if (m > (DBL_MAX_EXP-1)/PyLong_SHIFT) goto overflow; - errno = 0; - x = ldexp(x, e * PyLong_SHIFT); - if (Py_OVERFLOWED(x)) + + /* compute number of bits; no danger of overflow in the line below, + since we've already eliminated large m */ + n = m * PyLong_SHIFT; + assert(msd != 0); + do { + ++n; + msd >>= 1; + } while (msd != 0); + + /* catch almost all remaining overflow cases here */ + if (n > DBL_MAX_EXP) goto overflow; + + assert(n > DBL_MANT_DIG); /* dealt with |v| < 2**DBL_MANT_DIG above */ + rnd_bit = n - DBL_MANT_DIG - 1; + rnd_digit = rnd_bit/PyLong_SHIFT; + lsb = (digit)1 << (rnd_bit%PyLong_SHIFT); + + /* accumulate top DBL_MANT_DIG bits of v in x. */ + x = d[m]; + for (i = m-1; i > rnd_digit; i--) + x = x*PyLong_BASE + d[i]; + x = x*PyLong_BASE + (d[rnd_digit] & (PyLong_BASE-2*lsb)); + + /* apply round-half-to-even; round up if necessary */ + rounding_bit = d[rnd_digit] & lsb; + if (rounding_bit != 0) { + int round_up = 0; + digit parity_bit; + if (lsb == PyLong_BASE/2) + parity_bit = d[rnd_digit+1] & 1; + else + parity_bit = d[rnd_digit] & 2*lsb; + if (parity_bit != 0) + round_up = 1; + else if ((d[rnd_digit] & (lsb-1)) != 0) + round_up = 1; + else { + for (i = rnd_digit-1; i >= 0; i--) { + if (d[i] != 0) { + round_up = 1; + break; + } + } + } + if (round_up == 1) { + x += 2*lsb; + if (n == DBL_MAX_EXP && + x == ldexp((double)(2*lsb), DBL_MANT_DIG)) { + /* overflow corner case */ + goto overflow; + } + } + } + /* shift, adjust for sign, and return */ + x = sign * ldexp(x, rnd_digit*PyLong_SHIFT); return x; -overflow: + 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 67706) +++ Lib/test/test_long.py (working copy) @@ -615,6 +615,70 @@ else: self.assertRaises(TypeError, pow,longx, longy, long(z)) + def test_float_conversion(self): + # these tests are IEEE-754 specific; skip them on non IEEE + # machines (the code for converting from long to float is not + # IEEE-754 specific) + if not float.__getformat__("double").startswith("IEEE"): + return + + 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 = [-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)) + + # behaviour near extremes of floating-point range + long_dbl_max = long(DBL_MAX) + top_power = 2**DBL_MAX_EXP + halfway = (long_dbl_max + top_power)/2 + self.assertEqual(float(long_dbl_max), DBL_MAX) + self.assertEqual(float(long_dbl_max+1), DBL_MAX) + self.assertEqual(float(halfway-1), DBL_MAX) + self.assertRaises(OverflowError, float, halfway) + self.assertEqual(float(1-halfway), -DBL_MAX) + self.assertRaises(OverflowError, float, -halfway) + self.assertRaises(OverflowError, float, top_power-1) + self.assertRaises(OverflowError, float, top_power) + self.assertRaises(OverflowError, float, top_power+1) + self.assertRaises(OverflowError, float, 2*top_power-1) + self.assertRaises(OverflowError, float, 2*top_power) + self.assertRaises(OverflowError, float, top_power*top_power) + + + 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