Merged revisions 71772 via svnmerge from

svn+ssh://pythondev@svn.python.org/python/trunk

........
  r71772 | mark.dickinson | 2009-04-20 22:13:33 +0100 (Mon, 20 Apr 2009) | 5 lines

  Issue #3166: Make long -> float (and int -> float) conversions
  correctly rounded, using round-half-to-even.  This ensures that the
  value of float(n) doesn't depend on whether we're using 15-bit digits
  or 30-bit digits for Python longs.
........

Lib/test/test_long.py

@@ -620,6 +620,65 @@ class LongTest(unittest.TestCase):
                             else:
                                 self.assertRaises(TypeError, pow,longx, longy, int(z))
 
+    @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,
+                         2**53-2,
+                         2**53-1,
+                         2**53,
+                         2**53+2,
+                         2**54-4,
+                         2**54-2,
+                         2**54,
+                         2**54+4]
+        for x in exact_values:
+            self.assertEqual(float(x), x)
+            self.assertEqual(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 range(15):
+                self.assertEqual(int(float(2**p*(2**53+x))), 2**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 range(15):
+                self.assertEqual(int(float(2**p*(2**54+x))), 2**p*(2**54+y))
+
+        # behaviour near extremes of floating-point range
+        int_dbl_max = int(DBL_MAX)
+        top_power = 2**DBL_MAX_EXP
+        halfway = (int_dbl_max + top_power)//2
+        self.assertEqual(float(int_dbl_max), DBL_MAX)
+        self.assertEqual(float(int_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 range(100):
+            x = 2**p * (2**53 + 1) + 1
+            y = 2**p * (2**53 + 2)
+            self.assertEqual(int(float(x)), y)
+
+            x = 2**p * (2**53 + 1)
+            y = 2**p * 2**53
+            self.assertEqual(int(float(x)), y)
+
     def test_float_overflow(self):
         import math
 

Misc/NEWS

@@ -12,6 +12,9 @@ What's New in Python 3.1 beta 1?
 Core and Builtins
 -----------------
 
+- Issue #3166: Make long -> float (and int -> float) conversions
+  correctly rounded.
+
 - Issue #1869 (and many duplicates): make round(x, n) correctly
   rounded for a float x, by using the decimal <-> binary conversions
   from Python/dtoa.c.  As a consequence, (e.g.) round(x, 2) now

Objects/longobject.c

@@ -6,6 +6,7 @@
 #include "longintrepr.h"
 #include "structseq.h"
 
+#include <float.h>
 #include <ctype.h>
 #include <stddef.h>
 
@@ -100,6 +101,9 @@ maybe_small_long(PyLongObject *v)
 		if (PyErr_CheckSignals()) PyTryBlock \
 	}
 
+/* forward declaration */
+static int bits_in_digit(digit d);
+
 /* Normalize (remove leading zeros from) a long int object.
    Doesn't attempt to free the storage--in most cases, due to the nature
    of the algorithms used, this could save at most be one word anyway. */
@@ -962,33 +966,166 @@ _PyLong_AsScaledDouble(PyObject *vv, int *exponent)
 #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;
+	digit lsb, *d;
+	int round_up = 0;
 	double x;
 
 	if (vv == NULL || !PyLong_Check(vv)) {
 		PyErr_BadInternalCall();
-		return -1;
-	}
-	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)
+	}
+
+	/* Notes on the method: for simplicity, assume v is positive and >=
+	   2**DBL_MANT_DIG. (For negative v we just ignore the sign until the
+	   end; for small v no rounding is necessary.)  Write n for the number
+	   of bits in v, so that 2**(n-1) <= v < 2**n, and n > DBL_MANT_DIG.
+
+	   Some terminology: the *rounding bit* of v is the 1st bit of v that
+	   will be rounded away (bit n - DBL_MANT_DIG - 1); the *parity bit*
+	   is the bit immediately above.  The round-half-to-even rule says
+	   that we round up if the rounding bit is set, unless v is exactly
+	   halfway between two floats and the parity bit is zero.
+
+	   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 PyLong digit containing the rounding bit.  Then the
+	   bits of the digit d[rnd_digit] look something like:
+
+	              rounding bit
+	                  |
+	                  v
+	      msb -> sssssrttttttttt <- lsb
+	                 ^
+	                 |
+	              parity bit
+
+	   where '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 the lsb of d[rnd_digit+1].  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) gets the
+	   trailing bits, and d[rnd_digit] & lsb gives the rounding bit.
+
+	   We initialize the double x to the integer given by digits
+	   d[rnd_digit:m-1], but with the rounding bit and trailing bits of
+	   d[rnd_digit] masked out.  So the value of x comes from the top
+	   DBL_MANT_DIG bits of v, multiplied by 2*lsb.  Note that in the loop
+	   that produces x, all floating-point operations are exact (assuming
+	   that FLT_RADIX==2).  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.
+
+	   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;
+	d = v->ob_digit;
+	assert(d[m]);  /* v should be normalized */
+
+	/* fast path for case where 0 < abs(v) < 2**DBL_MANT_DIG */
+	if (m < DBL_MANT_DIG / PyLong_SHIFT ||
+	    (m == DBL_MANT_DIG / PyLong_SHIFT &&
+	     d[m] < (digit)1 << DBL_MANT_DIG%PyLong_SHIFT)) {
+		x = d[m];
+		while (--m >= 0)
+			x = x*PyLong_BASE + d[m];
+		return Py_SIZE(v) < 0 ? -x : x;
+	}
+
+	/* if m is huge then overflow immediately; otherwise, compute the
+	   number of bits n in v.  The condition below implies n (= #bits) >=
+	   m * PyLong_SHIFT + 1 > DBL_MAX_EXP, hence v >= 2**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))
+	n = m * PyLong_SHIFT + bits_in_digit(d[m]);
+	if (n > DBL_MAX_EXP)
 		goto overflow;
-	return x;
 
-overflow:
+	/* find location of rounding bit */
+	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);
+
+	/* Get top DBL_MANT_DIG bits of v.  Assumes PyLong_SHIFT <
+	   DBL_MANT_DIG, so we'll need bits from at least 2 digits of v. */
+	x = d[m];
+	assert(m > rnd_digit);
+	while (--m > rnd_digit)
+		x = x*PyLong_BASE + d[m];
+	x = x*PyLong_BASE + (d[m] & (PyLong_BASE-2*lsb));
+
+	/* decide whether to round up, using round-half-to-even */
+	assert(m == rnd_digit);
+	if (d[m] & lsb) { /* if (rounding bit is set) */
+		digit parity_bit;
+		if (lsb == PyLong_BASE/2)
+			parity_bit = d[m+1] & 1;
+		else
+			parity_bit = d[m] & 2*lsb;
+		if (parity_bit)
+			round_up = 1;
+		else if (d[m] & (lsb-1))
+			round_up = 1;
+		else {
+			while (--m >= 0) {
+				if (d[m]) {
+					round_up = 1;
+					break;
+				}
+			}
+		}
+	}
+
+	/* and round up if necessary */
+	if (round_up) {
+		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 = ldexp(x, rnd_digit*PyLong_SHIFT);
+	return Py_SIZE(v) < 0 ? -x : x;
+
+  overflow:
 	PyErr_SetString(PyExc_OverflowError,
 		"Python int too large to convert to C double");
 	return -1.0;