X-Git-Url: https://git.saurik.com/apple/libc.git/blobdiff_plain/224c70764cab4e0e39a26aaf3ad3016552f62f55..2650fa9ee9806a25904566dea091b1225d74f063:/gdtoa/_ldbl_util.c

diff --git a/gdtoa/_ldbl_util.c b/gdtoa/_ldbl_util.c
index 1c55d27..5028b20 100644
--- 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;
 }