[pypy-svn] pypy default: Implement _AsDouble() with the "maximum precision" algorithm from

[arigo]

Author: Armin Rigo <arigo at tunes.org>
Branch: 
Changeset: r40888:35de438ccb66
Date: 2011-01-18 17:18 +0100
http://bitbucket.org/pypy/pypy/changeset/35de438ccb66/

Log:	Implement _AsDouble() with the "maximum precision" algorithm from
	CPython 2.7.

diff --git a/pypy/rlib/test/test_rbigint.py b/pypy/rlib/test/test_rbigint.py
--- a/pypy/rlib/test/test_rbigint.py
+++ b/pypy/rlib/test/test_rbigint.py
@@ -176,6 +176,24 @@
         d = f2.tofloat()
         assert d == float(2097152 << SHIFT)
 
+    def test_tofloat_precision(self):
+        assert rbigint.fromlong(0).tofloat() == 0.0
+        for sign in [1, -1]:
+            for p in xrange(100):
+                x = long(2**p * (2**53 + 1) + 1) * sign
+                y = long(2**p * (2**53+ 2)) * sign
+                rx = rbigint.fromlong(x)
+                rxf = rx.tofloat()
+                assert rxf == float(y)
+                assert rbigint.fromfloat(rxf).tolong() == y
+                #
+                x = long(2**p * (2**53 + 1)) * sign
+                y = long(2**p * 2**53) * sign
+                rx = rbigint.fromlong(x)
+                rxf = rx.tofloat()
+                assert rxf == float(y)
+                assert rbigint.fromfloat(rxf).tolong() == y
+
     def test_fromfloat(self):
         x = 1234567890.1234567890
         f1 = rbigint.fromfloat(x)

diff --git a/pypy/rlib/rbigint.py b/pypy/rlib/rbigint.py
--- a/pypy/rlib/rbigint.py
+++ b/pypy/rlib/rbigint.py
@@ -1,6 +1,7 @@
 from pypy.rlib.rarithmetic import LONG_BIT, intmask, r_uint, r_ulonglong
 from pypy.rlib.rarithmetic import ovfcheck, r_longlong, widen, isinf, isnan
 from pypy.rlib.debug import make_sure_not_resized
+from pypy.rlib.objectmodel import we_are_translated
 
 import math, sys
 
@@ -1292,15 +1293,54 @@
 # XXX make sure that we don't ignore this!
 # YYY no, we decided to do ignore this!
 
-def _AsDouble(v):
+def _AsDouble(n):
     """ Get a C double from a bigint object. """
-    x, e = _AsScaledDouble(v)
-    if e <= sys.maxint / SHIFT:
-        x = math.ldexp(x, e * SHIFT)
-        #if not isinf(x):
-        # this is checked by math.ldexp
-        return x
-    raise OverflowError # can't say "long int too large to convert to float"
+    # This is a "correctly-rounded" version from Python 2.7.
+    #
+    from pypy.rlib import rfloat
+    DBL_MANT_DIG = rfloat.DBL_MANT_DIG  # 53 for IEEE 754 binary64
+    DBL_MAX_EXP = rfloat.DBL_MAX_EXP    # 1024 for IEEE 754 binary64
+    assert DBL_MANT_DIG < r_ulonglong.BITS
+
+    # Reduce to case n positive.
+    sign = n.sign
+    if sign == 0:
+        return 0.0
+    elif sign < 0:
+        n = n.neg()
+
+    # Find exponent: 2**(exp - 1) <= n < 2**exp
+    exp = n.bit_length()
+
+    # Get top DBL_MANT_DIG + 2 significant bits of n, with a 'sticky'
+    # last bit: that is, the least significant bit of the result is 1
+    # iff any of the shifted-out bits is set.
+    shift = DBL_MANT_DIG + 2 - exp
+    if shift >= 0:
+        q = _AsULonglong_mask(n) << shift
+        if not we_are_translated():
+            assert q == n.tolong() << shift   # no masking actually done
+    else:
+        shift = -shift
+        n2 = n.rshift(shift)
+        q = _AsULonglong_mask(n2)
+        if not we_are_translated():
+            assert q == n2.tolong()           # no masking actually done
+        if not n.eq(n2.lshift(shift)):
+            q |= 1
+
+    # Now remove the excess 2 bits, rounding to nearest integer (with
+    # ties rounded to even).
+    q = (q >> 2) + (bool(q & 2) and bool(q & 5))
+
+    if (exp > DBL_MAX_EXP or (exp == DBL_MAX_EXP and
+                              q == r_ulonglong(2) ** DBL_MANT_DIG)):
+        raise OverflowError("integer too large to convert to float")
+
+    ad = math.ldexp(float(q), exp - DBL_MANT_DIG)
+    if sign < 0:
+        ad = -ad
+    return ad
 
 def _loghelper(func, arg):
     """