Libc-1158.50.2.tar.gz

[apple/libc.git] / gdtoa / _ldbl_util.c

diff --git a/gdtoa/_ldbl_util.c b/gdtoa/_ldbl_util.c
index 1c55d27cb4c6582b24fc981437aff10c9f6a11a0..5028b20356e781422b9f472db711be5aae98261c 100644 (file)
--- a/gdtoa/_ldbl_util.c
+++ b/gdtoa/_ldbl_util.c
@@ -1,5 +1,5 @@
 /*
- * Copyright (c) 2004 Apple Computer, Inc. All rights reserved.
+ * Copyright (c) 2004, 2008 Apple Inc. All rights reserved.
  *
  * @APPLE_LICENSE_HEADER_START@
  * 
@@ -26,7 +26,6 @@
 #include <strings.h>
 #include <float.h>
 #include <math.h>
-#include <alloca.h>
 
 #include "fpmath.h"
 
@@ -36,23 +35,40 @@
 #define        LL_BITS         (8 * sizeof(int64_t))
 #define        LL_HIGHBIT      (1LL << 63)
 
-__private_extern__ void
+__private_extern__ int
 _ldbl2array32dd(union IEEEl2bits u, uint32_t *a)
 {
        int bit, shift, highbit, dexp;
        uint64_t a64[2];
        int64_t t64;
+       int extrabit = 0;
+
+       if(u.d[0] == 0.0) {
+           a[0] = a[1] = a[2] = a[3] = 0;
+           return 0;
+       }
 
        bzero(a64, sizeof(a64));
 
        switch (__fpclassifyd(u.d[0])) {
        case FP_NORMAL:
-               a64[1] = (1LL << (LDBL_MANT_DIG - BITS64 - 1));
-               /* drop through */
+               /*
+                * special case: if the head double only has the high (hidden)
+                * bit set, and the tail double is non-zero and is opposite
+                * in sign, then we increment extrabit to keep 106 bit
+                * precision in the results.
+                */
+               if(u.bits.manh == 0 && u.d[1] != 0 && u.bits.sign != u.bits.sign2)
+                   extrabit++;
+               a64[1] = (1LL << (LDBL_MANT_DIG - BITS64 - 1 + extrabit));
+               a64[1] |= ((uint64_t)u.bits.manh >> (BITS64 - LDBL_MANL_SIZE - extrabit));
+               a64[0] = ((uint64_t)u.bits.manh << (LDBL_MANL_SIZE + extrabit));
+               break;
        case FP_SUBNORMAL:
                a64[1] |= ((uint64_t)u.bits.manh >> (BITS64 - LDBL_MANL_SIZE));
                a64[0] = ((uint64_t)u.bits.manh << LDBL_MANL_SIZE);
-               break;
+               /* the tail double will be zero, so we are done */
+               goto done;
        default:
                goto done;
        }
@@ -62,20 +78,30 @@ _ldbl2array32dd(union IEEEl2bits u, uint32_t *a)
         * if the tail double is so small to not fit in LDBL_MANT_DIG bits,
         * then just skip it.
         */
-       if (dexp >= LDBL_MANT_DIG)
+       if (dexp >= LDBL_MANT_DIG + extrabit) {
+               reshift:
+               if (extrabit) {
+                       bit = a64[1] & 1;
+                       a64[1] >>= 1;
+                       a64[0] >>= 1;
+                       a64[0] |= ((uint64_t)bit) << (BITS64 - 1);
+                       extrabit = 0;
+               }
                goto done;
+       }
 
        switch (__fpclassifyd(u.d[1])) {
        case FP_NORMAL:
-               bit = LDBL_MANT_DIG - dexp - 1;
+               bit = LDBL_MANT_DIG - dexp - 1 + extrabit;
                t64 = (1LL << bit);
                break;
        case FP_SUBNORMAL:
-               bit = LDBL_MANT_DIG - (int)u.bits.exp;
+               bit = LDBL_MANT_DIG - (int)u.bits.exp + extrabit;
                t64 = 0;
                break;
        default:
-               goto done;
+               /* should never get here */
+               goto reshift;
        }
        shift = LDBL_MANL_SIZE - bit - 1;
        if (shift >= 0)
@@ -98,4 +124,5 @@ _ldbl2array32dd(union IEEEl2bits u, uint32_t *a)
        a[1] = (uint32_t)(a64[0] >> 32);
        a[2] = (uint32_t)a64[1];
        a[3] = (uint32_t)(a64[1] >> 32);
+       return extrabit;
 }