[commit: integer-gmp] master: Eliminate intermediate overflow for encodeFloat, fixes #5524 (3b9a78d)

[Daniel Fischer]

Repository : ssh://darcs.haskell.org//srv/darcs/packages/integer-gmp

On branch  : master

http://hackage.haskell.org/trac/ghc/changeset/3b9a78d45eab41ecd616e805a3adc25d5d3c8332

>---------------------------------------------------------------

commit 3b9a78d45eab41ecd616e805a3adc25d5d3c8332
Author: Daniel Fischer <daniel.is.fischer at googlemail.com>
Date:   Sat Oct 1 13:59:55 2011 +0200

    Eliminate intermediate overflow for encodeFloat, fixes #5524

>---------------------------------------------------------------

 cbits/float.c |   94 ++++++++++++++++++++++++++++++++++++++++++++++++++++++--
 1 files changed, 90 insertions(+), 4 deletions(-)

diff --git a/cbits/float.c b/cbits/float.c
index e2ef91a..ec2ec13 100644
--- a/cbits/float.c
+++ b/cbits/float.c
@@ -33,8 +33,10 @@
 
 #if SIZEOF_LIMB_T == 4
 #define GMP_BASE 4294967296.0
+#define LIMBBITS_LOG_2 5
 #elif SIZEOF_LIMB_T == 8
 #define GMP_BASE 18446744073709551616.0
+#define LIMBBITS_LOG_2 6
 #else
 #error Cannot cope with SIZEOF_LIMB_T -- please add definition of GMP_BASE
 #endif
@@ -71,8 +73,49 @@ integer_cbits_encodeDouble (I_ size, StgByteArray ba, I_ e) /* result = s * 2^e
     I_ i;
 
     /* Convert MP_INT to a double; knows a lot about internal rep! */
-    for(r = 0.0, i = __abs(size)-1; i >= 0; i--)
-        r = (r * GMP_BASE) + arr[i];
+    i = __abs(size)-1;
+    if ((i < 15) || (e >= 0)) /* overflows only if the final result does */
+    {
+        /* This would cause overflow if a large MP_INT is passed, even if the
+         * exponent would scale it back into range, so we do it only when it's safe. */
+        for(r = 0.0; i >= 0; i--)
+            r = (r * GMP_BASE) + arr[i];
+
+    } else { /* possibly more than 1024 bits in the MP_INT, but gets scaled down */
+
+        /* Find the first nonzero limb; normally it would be the first */
+        r = 0.0;
+        while((i >= 0) && (r == 0.0))
+        {
+            r = arr[i--];
+        }
+        if (i >= 0)
+            r = (r * GMP_BASE) + arr[i];
+#if SIZEOF_LIMB_T < 8
+        if (i > 0)
+            r = (r * GMP_BASE) + arr[--i];
+#endif
+        /* Now we have at least the 65 leading bits of the MP_INT or all of it.
+         * Any further bits would be rounded down, so from now on everything is
+         * multiplication by powers of 2.
+         * If i is positive, arr contains i limbs we haven't looked at yet, so
+         * adjust the exponent by i*8*SIZEOF_LIMB_T. Unfortunately, we must
+         * beware of overflow, so we can't simply add this to e. */
+        if (i > 0)
+        {
+            /* first add the number of whole limbs that would be cancelled */
+            i = i + e / (8 * SIZEOF_LIMB_T);
+            /* check for overflow */
+            if ((i > 0) && ((i >> (8*sizeof(I_) - 1 - LIMBBITS_LOG_2)) > 0))
+            {
+                /* overflow, give e a large dummy value */
+                e = 2147483647;
+            } else {
+                /* no overflow, get the exact value */
+                e = i * (8 * SIZEOF_LIMB_T) + (e % (8 * SIZEOF_LIMB_T));
+            }
+        }
+    }
 
     /* Now raise to the exponent */
     if ( r != 0.0 ) /* Lennart suggests this avoids a bug in MIPS's ldexp */
@@ -93,8 +136,51 @@ integer_cbits_encodeFloat (I_ size, StgByteArray ba, I_ e) /* result = s * 2^e *
     I_ i;
 
     /* Convert MP_INT to a float; knows a lot about internal rep! */
-    for(r = 0.0, i = __abs(size)-1; i >= 0; i--)
-        r = (r * GMP_BASE) + arr[i];
+    i = __abs(size)-1;
+    /* just in case StgFloat is a double, check sizes */
+#if SIZEOF_FLOAT == 4
+    if ((i < 2) || (e >= 0))
+#else
+    if ((i < 15) || (e >= 0))
+#endif
+    {
+        for(r = 0.0; i >= 0; i--)
+            r = (r * GMP_BASE) + arr[i];
+    } else {
+
+        /* Find the first nonzero limb; normally it would be the first */
+        r = 0.0;
+        while((i >= 0) && (r == 0.0))
+        {
+            r = arr[i--];
+        }
+        if (i >= 0)
+            r = (r * GMP_BASE) + arr[i];
+#if (SIZEOF_LIMB_T < 8) && (SIZEOF_FLOAT > 4)
+        if (i > 0)
+            r = (r * GMP_BASE) + arr[--i];
+#endif
+        /* Now we have enough leading bits of the MP_INT.
+         * Any further bits would be rounded down, so from now on everything is
+         * multiplication by powers of 2.
+         * If i is positive, arr contains i limbs we haven't looked at yet, so
+         * adjust the exponent by i*8*SIZEOF_LIMB_T. Unfortunately, we must
+         * beware of overflow, so we can't simply add this to e. */
+        if (i > 0)
+        {
+            /* first add the number of whole limbs that would be cancelled */
+            i = i + e / (8 * SIZEOF_LIMB_T);
+            /* check for overflow */
+            if ((i > 0) && ((i >> (8*sizeof(I_) - 1 - LIMBBITS_LOG_2)) > 0))
+            {
+                /* overflow, give e a large dummy value */
+                e = 2147483647;
+            } else {
+                /* no overflow, get the exact value */
+                e = i * (8 * SIZEOF_LIMB_T) + (e % (8 * SIZEOF_LIMB_T));
+            }
+        }
+    }
 
     /* Now raise to the exponent */
     if ( r != 0.0 ) /* Lennart suggests this avoids a bug in MIPS's ldexp */