[pypy-commit] pypy numpypy-longdouble: implement uint64-tuple representation of long double, break inf, nan handling

[mattip]

Author: mattip <matti.picus at gmail.com>
Branch: numpypy-longdouble
Changeset: r59550:b504a5a69cf9
Date: 2012-12-25 03:14 +0200
http://bitbucket.org/pypy/pypy/changeset/b504a5a69cf9/

Log:	implement uint64-tuple representation of long double, break inf,nan
	handling

diff --git a/pypy/rlib/rstruct/ieee.py b/pypy/rlib/rstruct/ieee.py
--- a/pypy/rlib/rstruct/ieee.py
+++ b/pypy/rlib/rstruct/ieee.py
@@ -26,6 +26,8 @@
     return int_part
 
 def float_unpack80(QQ, size):
+    '''Unpack a (mant, exp) tuple of r_ulonglong into a long double float
+    '''
     if size == 16 or size == 12:
         #Implement a x86-hardware extended 80 bit format as two 64 bit uints
         # QQ[0] is the sign and exp, QQ[1] is the mant
@@ -42,13 +44,13 @@
         # This tests generates wrong code when translated:
         # with gcc, shifting a 64bit int by 64 bits does
         # not change the value.
-        if QQ[0] >> TOPBITS:
-            raise ValueError("input '%r' out of range '%r'" % (QQ, QQ[0]>>TOPBITS))
+        if QQ[1] >> TOPBITS:
+            raise ValueError("input '%r' out of range '%r'" % (QQ, QQ[1]>>TOPBITS))
 
     # extract pieces with explicit one in MANT_DIG
-    sign = rarithmetic.intmask(QQ[0] >> TOPBITS - 1)
-    exp = rarithmetic.intmask((QQ[0] & ((one << TOPBITS - 1) - 1)))
-    mant = QQ[1]
+    sign = rarithmetic.intmask(QQ[1] >> TOPBITS - 1)
+    exp = rarithmetic.intmask((QQ[1] & ((one << TOPBITS - 1) - 1)))
+    mant = QQ[0]
 
     if exp == MAX_EXP - MIN_EXP + 2:
         # nan or infinity
@@ -111,9 +113,6 @@
 def float_pack(x, size):
     """Convert a Python float x into a 64-bit unsigned integer
     with the same byte representation."""
-    return float_pack_helper(x, size, r_ulonglong)
-
-def float_pack_helper(x, size, r_ulonglong):
     if size == 8:
         MIN_EXP = -1021  # = sys.float_info.min_exp
         MAX_EXP = 1024   # = sys.float_info.max_exp
@@ -189,6 +188,65 @@
     sign = r_ulonglong(sign)
     return ((sign << BITS - 1) | (exp << MANT_DIG - 1)) | mant
 
+def float_pack80(x, size):
+    """Convert a Python float x into two 64-bit unsigned integers
+    with the same byte representation."""
+    if size == 16 or size == 12:
+        #Implement a x86-hardware extended 80 bit format
+        # with explicit 1 in bit 64
+        MIN_EXP = -16381
+        MAX_EXP = 16384
+        MANT_DIG = 64 
+        BITS = 80
+    else:
+        raise ValueError("invalid size value")
+
+    sign = rfloat.copysign(1.0, x) < 0.0
+    if not rfloat.isfinite(x):
+        if rfloat.isinf(x):
+            mant = r_ulonglong(0)
+            exp = MAX_EXP - MIN_EXP + 2
+        else:  # rfloat.isnan(x):
+            mant = r_ulonglong(1) << (MANT_DIG-2) # other values possible
+            exp = MAX_EXP - MIN_EXP + 2
+    elif x == 0.0:
+        mant = r_ulonglong(0)
+        exp = 0
+    else:
+        m, e = math.frexp(abs(x))  # abs(x) == m * 2**e
+        exp = e - (MIN_EXP - 1)
+        if exp > 0:
+            # Normal case. Avoid uint64 overflow by using MANT_DIG-1
+            mant = round_to_nearest(m * (r_ulonglong(1) << MANT_DIG - 1))
+        else:
+            # Subnormal case.
+            if exp + MANT_DIG - 1 >= 0:
+                mant = round_to_nearest(m * (r_ulonglong(1) << exp + MANT_DIG - 1))
+            else:
+                mant = r_ulonglong(0)
+            exp = 0
+
+        # Special case: rounding produced a MANT_DIG-bit mantissa.
+        if mant == r_ulonglong(1) << MANT_DIG - 1:
+            mant = r_ulonglong(0)
+            exp += 1
+
+        # Raise on overflow (in some circumstances, may want to return
+        # infinity instead).
+        if exp >= MAX_EXP - MIN_EXP + 2:
+             raise OverflowError("float too large to pack in this format")
+
+    # check constraints
+    if not objectmodel.we_are_translated():
+        assert 0 <= mant < 1 << MANT_DIG - 1
+        assert 0 <= exp <= MAX_EXP - MIN_EXP + 2
+        assert 0 <= sign <= 1
+    if size==12 or size == 16:
+        mant = mant << 1
+    exp = r_ulonglong(exp)
+    sign = r_ulonglong(sign)
+    return (mant, (sign << BITS - MANT_DIG - 1) | exp)
+
 
 @jit.unroll_safe
 def pack_float(result, x, size, be):
diff --git a/pypy/rlib/rstruct/test/test_ieee.py b/pypy/rlib/rstruct/test/test_ieee.py
--- a/pypy/rlib/rstruct/test/test_ieee.py
+++ b/pypy/rlib/rstruct/test/test_ieee.py
@@ -3,7 +3,7 @@
 import struct
 
 from pypy.rlib.rfloat import isnan
-from pypy.rlib.rstruct.ieee import float_pack, float_unpack, float_pack128, float_unpack80
+from pypy.rlib.rstruct.ieee import float_pack, float_unpack, float_pack80, float_unpack80
 
 
 class TestFloatPacking:
@@ -18,10 +18,10 @@
         y = float_unpack(Q, 8)
         assert repr(x) == repr(y)
 
-        Q = float_pack128(x, 16)
+        Q = float_pack80(x, 16)
         y = float_unpack80(Q, 16)
         assert repr(x) == repr(y),'%r != %r, Q=%r'%(x, y, Q)
-        Q = float_pack128(x, 12)
+        Q = float_pack80(x, 12)
         y = float_unpack80(Q, 12)
         assert repr(x) == repr(y),'%r != %r, Q=%r'%(x, y, Q)
 
@@ -65,10 +65,10 @@
         self.check_float(-0.0)
 
     def test_nans(self):
-        Q = float_pack128(float('nan'), 16)
+        Q = float_pack80(float('nan'), 16)
         y = float_unpack80(Q, 16)
         assert repr(y) == 'nan'
-        Q = float_pack128(float('nan'), 12)
+        Q = float_pack80(float('nan'), 12)
         y = float_unpack80(Q, 12)
         assert repr(y) == 'nan'
         Q = float_pack(float('nan'), 8)