This is the mail archive of the libc-alpha@sourceware.org mailing list for the glibc project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

A little more precise %Lf for IBM long double


I'd like to print IBM long doubles with a few extra bits of precision,
because IBM long double has variable precision and we often have
values with 107 bits during the normal course of calculations.  It's
possible to have an IBM long double with 2k bits of precision, and
while some might argue we ought to print to the full precision, I'm
not inclined to try to implement that.

So this patch gives the mpn values returned from ldbl2mpn() an extra
11 bits of precision.  To do that it's necessary to inform the caller
of __mpn_extract_long_double() exactly how many bits are returned
rather than assuming a fixed LDBL_MANT_DIG bits.  I did that
indirectly by returning the number of zero bits in the last mp_limb,
which turns out to be convenient both in __mpn_extract_long_double(),
and it's caller.  Given the extra parameter it then becomes possible
to omit shifting in __mpn_extract_long_double(), tidying the code a
little, since the caller does shifting anyway.

In the following, note that 1 - IEEE854_LONG_DOUBLE_BIAS is equal to
LDBL_MIN_EXP - 1.  I prefer the former since we already use
IEEE854_LONG_DOUBLE_BIAS in the exponent calculation for normalised
values, and it reinforces the fact that denormals are treated as if
their unbiased exponent was 1.

	* include/gmp.h (__mpn_extract_double, __mpn_extract_long_double):
	Update prototypes.
	* stdio-common/printf_fp.c: Likewise.
	(__printf_fp): Use returned zero_bits.
	* stdlib/dbl2mpn.c (__mpn_extract_double): Update funtion arguments.
	* sysdeps/i386/ldbl2mpn.c (__mpn_extract_long_double): Don't perform
	shifts for denormals here.  Instead return leading zero bit count
	and adjust return value of function.
	* sysdeps/ieee754/dbl-64/dbl2mpn.c (__mpn_extract_long_double): Likewise.
	* sysdeps/ieee754/ldbl-128/ldbl2mpn.c (__mpn_extract_long_double):
	Likewise.
	* sysdeps/ieee754/ldbl-96/ldbl2mpn.c (__mpn_extract_long_double):
	Likewise.
	* sysdeps/ieee754/ldbl-128ibm/ldbl2mpn.c (__mpn_extract_long_double):
	Likewise.  Return an extra 11 bits of precision.

diff --git a/include/gmp.h b/include/gmp.h
index b741670..bcc9360 100644
--- a/include/gmp.h
+++ b/include/gmp.h
@@ -8,12 +8,12 @@
 
 /* Now define the internal interfaces.  */
 extern mp_size_t __mpn_extract_double (mp_ptr res_ptr, mp_size_t size,
-				       int *expt, int *is_neg,
-				       double value);
+				       int *expt, int *zero_bits,
+				       int *is_neg, double value);
 
 extern mp_size_t __mpn_extract_long_double (mp_ptr res_ptr, mp_size_t size,
-					    int *expt, int *is_neg,
-					    long double value);
+					    int *expt, int *zero_bits,
+					    int *is_neg, long double value);
 
 extern float __mpn_construct_float (mp_srcptr frac_ptr, int expt, int sign);
 
diff --git a/stdio-common/printf_fp.c b/stdio-common/printf_fp.c
index 2b93e6c..c0a7a82 100644
--- a/stdio-common/printf_fp.c
+++ b/stdio-common/printf_fp.c
@@ -134,11 +134,11 @@
   (u##size > v##size || (u##size == v##size && __mpn_cmp (u, v, u##size) >= 0))
 
 extern mp_size_t __mpn_extract_double (mp_ptr res_ptr, mp_size_t size,
-				       int *expt, int *is_neg,
-				       double value);
+				       int *expt, int *zero_bits,
+				       int *is_neg, double value);
 extern mp_size_t __mpn_extract_long_double (mp_ptr res_ptr, mp_size_t size,
-					    int *expt, int *is_neg,
-					    long double value);
+					    int *expt, int *zero_bits,
+					    int *is_neg, long double value);
 extern unsigned int __guess_grouping (unsigned int intdig_max,
 				      const char *grouping);
 
@@ -360,12 +360,13 @@ ___printf_fp (FILE *fp,
 	}
       else
 	{
+	  int zero_bits;
 	  fracsize = __mpn_extract_long_double (fp_input,
 						(sizeof (fp_input) /
 						 sizeof (fp_input[0])),
-						&exponent, &is_neg,
+						&exponent, &zero_bits, &is_neg,
 						fpnum.ldbl);
-	  to_shift = 1 + fracsize * BITS_PER_MP_LIMB - LDBL_MANT_DIG;
+	  to_shift = 1 + zero_bits;
 	}
     }
   else
@@ -406,11 +407,13 @@ ___printf_fp (FILE *fp,
 	}
       else
 	{
+	  int zero_bits;
 	  fracsize = __mpn_extract_double (fp_input,
 					   (sizeof (fp_input)
 					    / sizeof (fp_input[0])),
-					   &exponent, &is_neg, fpnum.dbl);
-	  to_shift = 1 + fracsize * BITS_PER_MP_LIMB - DBL_MANT_DIG;
+					   &exponent, &zero_bits, &is_neg,
+					   fpnum.dbl);
+	  to_shift = 1 + zero_bits;
 	}
     }
 
diff --git a/stdlib/dbl2mpn.c b/stdlib/dbl2mpn.c
index 429e20a..dd835f3 100644
--- a/stdlib/dbl2mpn.c
+++ b/stdlib/dbl2mpn.c
@@ -24,7 +24,7 @@
 
 mp_size_t
 __mpn_extract_double (mp_ptr res_ptr, mp_size_t size,
-		      int *expt, int *is_neg,
+		      int *expt, int *zero_bits, int *is_neg,
 		      double value)
 {
 #error "__mpn_extract_double is not implemented for this floating point format"
diff --git a/sysdeps/i386/ldbl2mpn.c b/sysdeps/i386/ldbl2mpn.c
index c7b322b..968aea5 100644
--- a/sysdeps/i386/ldbl2mpn.c
+++ b/sysdeps/i386/ldbl2mpn.c
@@ -29,7 +29,7 @@
 
 mp_size_t
 __mpn_extract_long_double (mp_ptr res_ptr, mp_size_t size,
-			   int *expt, int *is_neg,
+			   int *expt, int *zero_bits, int *is_neg,
 			   long double value)
 {
   union ieee854_long_double u;
@@ -51,18 +51,26 @@ __mpn_extract_long_double (mp_ptr res_ptr, mp_size_t size,
   #error "mp_limb size " BITS_PER_MP_LIMB "not accounted for"
 #endif
 
+/* The format does not fill the last limb.  There are some zeros.  */
+#define NUM_LEADING_ZEROS (N * BITS_PER_MP_LIMB - LDBL_MANT_DIG)
+  *zero_bits = NUM_LEADING_ZEROS;
+
   if (u.ieee.exponent == 0)
     {
       /* A biased exponent of zero is a special case.
 	 Either it is a zero or it is a denormal number.  */
       if (res_ptr[0] == 0 && res_ptr[N - 1] == 0) /* Assumes N<=2.  */
-	/* It's zero.  */
-	*expt = 0;
+	{
+	  /* It's zero.  */
+	  *expt = 0;
+	  return 1;
+	}
       else
 	{
 	  /* It is a denormal number, meaning it has no implicit leading
 	     one bit, and its exponent is in fact the format minimum.  */
 	  int cnt;
+	  int exp = 1 - IEEE854_LONG_DOUBLE_BIAS;
 
 	  /* One problem with Intel's 80-bit format is that the explicit
 	     leading one in the normalized representation has to be zero
@@ -74,24 +82,17 @@ __mpn_extract_long_double (mp_ptr res_ptr, mp_size_t size,
 	  if (res_ptr[N - 1] != 0)
 	    {
 	      count_leading_zeros (cnt, res_ptr[N - 1]);
-	      if (cnt != 0)
-		{
-#if N == 2
-		  res_ptr[N - 1] = res_ptr[N - 1] << cnt
-				   | (res_ptr[0] >> (BITS_PER_MP_LIMB - cnt));
-		  res_ptr[0] <<= cnt;
-#else
-		  res_ptr[N - 1] <<= cnt;
-#endif
-		}
-	      *expt = LDBL_MIN_EXP - 1 - cnt;
+	      *zero_bits = cnt;
+	      *expt = exp + NUM_LEADING_ZEROS - cnt;
+	      return N;
 	    }
 	  else if (res_ptr[0] != 0)
 	    {
 	      count_leading_zeros (cnt, res_ptr[0]);
-	      res_ptr[N - 1] = res_ptr[0] << cnt;
-	      res_ptr[0] = 0;
-	      *expt = LDBL_MIN_EXP - 1 - BITS_PER_MP_LIMB - cnt;
+	      exp -= BITS_PER_MP_LIMB;
+	      *zero_bits = cnt;
+	      *expt = exp + NUM_LEADING_ZEROS - cnt;
+	      return 1;
 	    }
 	  else
 	    {
@@ -104,7 +105,7 @@ __mpn_extract_long_double (mp_ptr res_ptr, mp_size_t size,
 #else
 	      res_ptr[0] = 0x8000000000000000ul;
 #endif
-	      *expt = LDBL_MIN_EXP - 1;
+	      *expt = 1 - IEEE854_LONG_DOUBLE_BIAS;
 	    }
 	}
     }
diff --git a/sysdeps/ieee754/dbl-64/dbl2mpn.c b/sysdeps/ieee754/dbl-64/dbl2mpn.c
index f2294de..cae062a 100644
--- a/sysdeps/ieee754/dbl-64/dbl2mpn.c
+++ b/sysdeps/ieee754/dbl-64/dbl2mpn.c
@@ -28,7 +28,7 @@
 
 mp_size_t
 __mpn_extract_double (mp_ptr res_ptr, mp_size_t size,
-		      int *expt, int *is_neg,
+		      int *expt, int *zero_bits, int *is_neg,
 		      double value)
 {
   union ieee754_double u;
@@ -49,9 +49,10 @@ __mpn_extract_double (mp_ptr res_ptr, mp_size_t size,
 #else
   #error "mp_limb size " BITS_PER_MP_LIMB "not accounted for"
 #endif
+
 /* The format does not fill the last limb.  There are some zeros.  */
-#define NUM_LEADING_ZEROS (BITS_PER_MP_LIMB \
-			   - (DBL_MANT_DIG - ((N - 1) * BITS_PER_MP_LIMB)))
+#define NUM_LEADING_ZEROS (N * BITS_PER_MP_LIMB - DBL_MANT_DIG)
+  *zero_bits = NUM_LEADING_ZEROS;
 
   if (u.ieee.exponent == 0)
     {
@@ -65,37 +66,20 @@ __mpn_extract_double (mp_ptr res_ptr, mp_size_t size,
           /* It is a denormal number, meaning it has no implicit leading
 	     one bit, and its exponent is in fact the format minimum.  */
 	  int cnt;
+	  int exp = 1 - IEEE754_DOUBLE_BIAS;
+	  int n = N;
 
 	  if (res_ptr[N - 1] != 0)
-	    {
-	      count_leading_zeros (cnt, res_ptr[N - 1]);
-	      cnt -= NUM_LEADING_ZEROS;
-#if N == 2
-	      res_ptr[N - 1] = res_ptr[1] << cnt
-			       | (N - 1)
-			         * (res_ptr[0] >> (BITS_PER_MP_LIMB - cnt));
-	      res_ptr[0] <<= cnt;
-#else
-	      res_ptr[N - 1] <<= cnt;
-#endif
-	      *expt = DBL_MIN_EXP - 1 - cnt;
-	    }
+	    count_leading_zeros (cnt, res_ptr[N - 1]);
 	  else
 	    {
 	      count_leading_zeros (cnt, res_ptr[0]);
-	      if (cnt >= NUM_LEADING_ZEROS)
-		{
-		  res_ptr[N - 1] = res_ptr[0] << (cnt - NUM_LEADING_ZEROS);
-		  res_ptr[0] = 0;
-		}
-	      else
-		{
-		  res_ptr[N - 1] = res_ptr[0] >> (NUM_LEADING_ZEROS - cnt);
-		  res_ptr[0] <<= BITS_PER_MP_LIMB - (NUM_LEADING_ZEROS - cnt);
-		}
-	      *expt = DBL_MIN_EXP - 1
-		      - (BITS_PER_MP_LIMB - NUM_LEADING_ZEROS) - cnt;
+	      exp -= BITS_PER_MP_LIMB;
+	      n = 1;
 	    }
+	  *zero_bits = cnt;
+	  *expt = exp + NUM_LEADING_ZEROS - cnt;
+	  return n;
 	}
     }
   else
diff --git a/sysdeps/ieee754/ldbl-128/ldbl2mpn.c b/sysdeps/ieee754/ldbl-128/ldbl2mpn.c
index 64f98a5..8e2145b 100644
--- a/sysdeps/ieee754/ldbl-128/ldbl2mpn.c
+++ b/sysdeps/ieee754/ldbl-128/ldbl2mpn.c
@@ -30,7 +30,7 @@
 
 mp_size_t
 __mpn_extract_long_double (mp_ptr res_ptr, mp_size_t size,
-			   int *expt, int *is_neg,
+			   int *expt, int *zero_bits, int *is_neg,
 			   long double value)
 {
   union ieee854_long_double u;
@@ -54,9 +54,10 @@ __mpn_extract_long_double (mp_ptr res_ptr, mp_size_t size,
 #else
   #error "mp_limb size " BITS_PER_MP_LIMB "not accounted for"
 #endif
+
 /* The format does not fill the last limb.  There are some zeros.  */
-#define NUM_LEADING_ZEROS (BITS_PER_MP_LIMB \
-			   - (LDBL_MANT_DIG - ((N - 1) * BITS_PER_MP_LIMB)))
+#define NUM_LEADING_ZEROS (N * BITS_PER_MP_LIMB - LDBL_MANT_DIG)
+  *zero_bits = NUM_LEADING_ZEROS;
 
   if (u.ieee.exponent == 0)
     {
@@ -64,70 +65,42 @@ __mpn_extract_long_double (mp_ptr res_ptr, mp_size_t size,
 	 Either it is a zero or it is a denormal number.  */
       if (res_ptr[0] == 0 && res_ptr[1] == 0
           && res_ptr[N - 2] == 0 && res_ptr[N - 1] == 0) /* Assumes N<=4.  */
-	/* It's zero.  */
-	*expt = 0;
+	{
+	  /* It's zero.  */
+	  *expt = 0;
+	  return 1;
+	}
       else
 	{
           /* It is a denormal number, meaning it has no implicit leading
 	     one bit, and its exponent is in fact the format minimum.  */
 	  int cnt;
+	  int exp = 1 - IEEE854_LONG_DOUBLE_BIAS;
+	  int n = N;
 
 #if N == 2
 	  if (res_ptr[N - 1] != 0)
-	    {
-	      count_leading_zeros (cnt, res_ptr[N - 1]);
-	      cnt -= NUM_LEADING_ZEROS;
-	      res_ptr[N - 1] = res_ptr[N - 1] << cnt
-			       | (res_ptr[0] >> (BITS_PER_MP_LIMB - cnt));
-	      res_ptr[0] <<= cnt;
-	      *expt = LDBL_MIN_EXP - 1 - cnt;
-	    }
+	    count_leading_zeros (cnt, res_ptr[N - 1]);
 	  else
 	    {
 	      count_leading_zeros (cnt, res_ptr[0]);
-	      if (cnt >= NUM_LEADING_ZEROS)
-		{
-		  res_ptr[N - 1] = res_ptr[0] << (cnt - NUM_LEADING_ZEROS);
-		  res_ptr[0] = 0;
-		}
-	      else
-		{
-		  res_ptr[N - 1] = res_ptr[0] >> (NUM_LEADING_ZEROS - cnt);
-		  res_ptr[0] <<= BITS_PER_MP_LIMB - (NUM_LEADING_ZEROS - cnt);
-		}
-	      *expt = LDBL_MIN_EXP - 1
-		- (BITS_PER_MP_LIMB - NUM_LEADING_ZEROS) - cnt;
+	      exp -= BITS_PER_MP_LIMB;
+	      n = 1;
 	    }
 #else
-	  int j, k, l;
+	  int j;
 
 	  for (j = N - 1; j > 0; j--)
 	    if (res_ptr[j] != 0)
 	      break;
 
 	  count_leading_zeros (cnt, res_ptr[j]);
-	  cnt -= NUM_LEADING_ZEROS;
-	  l = N - 1 - j;
-	  if (cnt < 0)
-	    {
-	      cnt += BITS_PER_MP_LIMB;
-	      l--;
-	    }
-	  if (!cnt)
-	    for (k = N - 1; k >= l; k--)
-	      res_ptr[k] = res_ptr[k-l];
-	  else
-	    {
-	      for (k = N - 1; k > l; k--)
-		res_ptr[k] = res_ptr[k-l] << cnt
-			     | res_ptr[k-l-1] >> (BITS_PER_MP_LIMB - cnt);
-	      res_ptr[k--] = res_ptr[0] << cnt;
-	    }
-
-	  for (; k >= 0; k--)
-	    res_ptr[k] = 0;
-	  *expt = LDBL_MIN_EXP - 1 - l * BITS_PER_MP_LIMB - cnt;
+	  exp -= (N - 1 - j) * BITS_PER_MP_LIMB;
+	  n = j + 1;
 #endif
+	  *zero_bits = cnt;
+	  *expt = exp + NUM_LEADING_ZEROS - cnt;
+	  return n;
 	}
     }
   else
diff --git a/sysdeps/ieee754/ldbl-128ibm/ldbl2mpn.c b/sysdeps/ieee754/ldbl-128ibm/ldbl2mpn.c
index e46fde7..14f9ed3 100644
--- a/sysdeps/ieee754/ldbl-128ibm/ldbl2mpn.c
+++ b/sysdeps/ieee754/ldbl-128ibm/ldbl2mpn.c
@@ -30,157 +30,121 @@
 
 mp_size_t
 __mpn_extract_long_double (mp_ptr res_ptr, mp_size_t size,
-			   int *expt, int *is_neg,
+			   int *expt, int *zero_bits, int *is_neg,
 			   long double value)
 {
   union ibm_extended_long_double u;
-  unsigned long long hi, lo;
+  uint64_t hi, lo;
   int ediff;
 
   u.ld = value;
 
   *is_neg = u.d[0].ieee.negative;
   *expt = (int) u.d[0].ieee.exponent - IEEE754_DOUBLE_BIAS;
+#define NUM_LEADING_ZEROS (128 - (LDBL_MANT_DIG + 11))
+  *zero_bits = NUM_LEADING_ZEROS;
 
-  lo = ((long long) u.d[1].ieee.mantissa0 << 32) | u.d[1].ieee.mantissa1;
-  hi = ((long long) u.d[0].ieee.mantissa0 << 32) | u.d[0].ieee.mantissa1;
+  lo = ((uint64_t) u.d[1].ieee.mantissa0 << 32) | u.d[1].ieee.mantissa1;
+  hi = ((uint64_t) u.d[0].ieee.mantissa0 << 32) | u.d[0].ieee.mantissa1;
 
-  /* If the lower double is not a denormal or zero then set the hidden
-     53rd bit.  */
-  if (u.d[1].ieee.exponent != 0)
-    lo |= 1ULL << 52;
-  else
-    lo = lo << 1;
-
-  /* The lower double is normalized separately from the upper.  We may
-     need to adjust the lower manitissa to reflect this.  */
-  ediff = u.d[0].ieee.exponent - u.d[1].ieee.exponent - 53;
-  if (ediff > 0)
+  if (u.d[0].ieee.exponent != 0)
     {
-      if (ediff < 64)
-	lo = lo >> ediff;
+      /* If the high double is not a denormal or zero then set the hidden
+	 53rd bit.  */
+      hi |= (uint64_t) 1 << 52;
+
+      /* If the lower double is not a denormal or zero then set the hidden
+	 53rd bit.  */
+      if (u.d[1].ieee.exponent != 0)
+	lo |= (uint64_t) 1 << 52;
       else
-	lo = 0;
-    }
-  else if (ediff < 0)
-    lo = lo << -ediff;
-
-  /* The high double may be rounded and the low double reflects the
-     difference between the long double and the rounded high double
-     value.  This is indicated by a differnce between the signs of the
-     high and low doubles.  */
-  if (u.d[0].ieee.negative != u.d[1].ieee.negative
-      && lo != 0)
-    {
-      lo = (1ULL << 53) - lo;
-      if (hi == 0)
+	lo = lo << 1;
+
+      /* We currently only have 53 bits in lo.  Gain a few more bits
+	 of precision.  */
+      lo = lo << 11;
+
+      /* The lower double is normalized separately from the upper.  We may
+	 need to adjust the lower manitissa to reflect this.  */
+      ediff = u.d[0].ieee.exponent - u.d[1].ieee.exponent - 53;
+      if (ediff > 0)
 	{
-	  /* we have a borrow from the hidden bit, so shift left 1.  */
-	  hi = 0x0ffffffffffffeLL | (lo >> 51);
-	  lo = 0x1fffffffffffffLL & (lo << 1);
-	  (*expt)--;
+	  if (ediff < 64)
+	    lo = lo >> ediff;
+	  else
+	    lo = 0;
 	}
-      else
-	hi--;
-    }
+      else if (ediff < 0)
+	lo = lo << -ediff;
+
+      /* The high double may be rounded and the low double reflects the
+	 difference between the long double and the rounded high double
+	 value.  This is indicated by a differnce between the signs of the
+	 high and low doubles.  */
+      if (u.d[0].ieee.negative != u.d[1].ieee.negative
+	  && lo != 0)
+	{
+	  hi--;
+	  lo = -lo;
+	  if (hi < (uint64_t) 1 << 52)
+	    {
+	      /* We have a borrow from the hidden bit, so shift left 1.  */
+	      hi = (hi << 1) | (lo >> 63);
+	      lo = lo << 1;
+	      (*expt)--;
+	    }
+	}
+
 #if BITS_PER_MP_LIMB == 32
-  /* Combine the mantissas to be contiguous.  */
-  res_ptr[0] = lo;
-  res_ptr[1] = (hi << (53 - 32)) | (lo >> 32);
-  res_ptr[2] = hi >> 11;
-  res_ptr[3] = hi >> (32 + 11);
-  #define N 4
+      res_ptr[0] = lo;
+      res_ptr[1] = lo >> 32;
+      res_ptr[2] = hi;
+      res_ptr[3] = hi >> 32;
+      return 4;
 #elif BITS_PER_MP_LIMB == 64
-  /* Combine the two mantissas to be contiguous.  */
-  res_ptr[0] = (hi << 53) | lo;
-  res_ptr[1] = hi >> 11;
-  #define N 2
+      res_ptr[0] = lo;
+      res_ptr[1] = hi;
+      return 2;
 #else
-  #error "mp_limb size " BITS_PER_MP_LIMB "not accounted for"
+# error "mp_limb size " BITS_PER_MP_LIMB "not accounted for"
 #endif
-/* The format does not fill the last limb.  There are some zeros.  */
-#define NUM_LEADING_ZEROS (BITS_PER_MP_LIMB \
-			   - (LDBL_MANT_DIG - ((N - 1) * BITS_PER_MP_LIMB)))
+    }
 
-  if (u.d[0].ieee.exponent == 0)
+  /* The high double is a denormal or zero.  The low double must
+     be zero.  A denormal is interpreted as having a biased
+     exponent of 1.  */
+  res_ptr[0] = hi;
+#if BITS_PER_MP_LIMB == 32
+  res_ptr[1] = hi >> 32;
+#endif
+  if (hi == 0)
     {
-      /* A biased exponent of zero is a special case.
-	 Either it is a zero or it is a denormal number.  */
-      if (res_ptr[0] == 0 && res_ptr[1] == 0
-	  && res_ptr[N - 2] == 0 && res_ptr[N - 1] == 0) /* Assumes N<=4.  */
-	/* It's zero.  */
-	*expt = 0;
+      /* It's zero.  */
+      *expt = 0;
+      return 1;
+    }
+  else
+    {
+      int cnt;
+      int exp = 1 - IEEE754_DOUBLE_BIAS;
+      int n = 1;
+
+#if BITS_PER_MP_LIMB == 32
+      if (res_ptr[1] != 0)
+	{
+	  count_leading_zeros (cnt, res_ptr[1]);
+	  n = 2;
+	}
       else
 	{
-	  /* It is a denormal number, meaning it has no implicit leading
-	     one bit, and its exponent is in fact the format minimum.  We
-	     use DBL_MIN_EXP instead of LDBL_MIN_EXP below because the
-	     latter describes the properties of both parts together, but
-	     the exponent is computed from the high part only.  */
-	  int cnt;
-
-#if N == 2
-	  if (res_ptr[N - 1] != 0)
-	    {
-	      count_leading_zeros (cnt, res_ptr[N - 1]);
-	      cnt -= NUM_LEADING_ZEROS;
-	      res_ptr[N - 1] = res_ptr[N - 1] << cnt
-			       | (res_ptr[0] >> (BITS_PER_MP_LIMB - cnt));
-	      res_ptr[0] <<= cnt;
-	      *expt = DBL_MIN_EXP - 1 - cnt;
-	    }
-	  else
-	    {
-	      count_leading_zeros (cnt, res_ptr[0]);
-	      if (cnt >= NUM_LEADING_ZEROS)
-		{
-		  res_ptr[N - 1] = res_ptr[0] << (cnt - NUM_LEADING_ZEROS);
-		  res_ptr[0] = 0;
-		}
-	      else
-		{
-		  res_ptr[N - 1] = res_ptr[0] >> (NUM_LEADING_ZEROS - cnt);
-		  res_ptr[0] <<= BITS_PER_MP_LIMB - (NUM_LEADING_ZEROS - cnt);
-		}
-	      *expt = DBL_MIN_EXP - 1
-		- (BITS_PER_MP_LIMB - NUM_LEADING_ZEROS) - cnt;
-	    }
+	  count_leading_zeros (cnt, res_ptr[0]);
+	  exp -= BITS_PER_MP_LIMB;
+	}
 #else
-	  int j, k, l;
-
-	  for (j = N - 1; j > 0; j--)
-	    if (res_ptr[j] != 0)
-	      break;
-
-	  count_leading_zeros (cnt, res_ptr[j]);
-	  cnt -= NUM_LEADING_ZEROS;
-	  l = N - 1 - j;
-	  if (cnt < 0)
-	    {
-	      cnt += BITS_PER_MP_LIMB;
-	      l--;
-	    }
-	  if (!cnt)
-	    for (k = N - 1; k >= l; k--)
-	      res_ptr[k] = res_ptr[k-l];
-	  else
-	    {
-	      for (k = N - 1; k > l; k--)
-		res_ptr[k] = res_ptr[k-l] << cnt
-			     | res_ptr[k-l-1] >> (BITS_PER_MP_LIMB - cnt);
-	      res_ptr[k--] = res_ptr[0] << cnt;
-	    }
-
-	  for (; k >= 0; k--)
-	    res_ptr[k] = 0;
-	  *expt = DBL_MIN_EXP - 1 - l * BITS_PER_MP_LIMB - cnt;
+      count_leading_zeros (cnt, hi);
 #endif
-	}
+      *zero_bits = cnt;
+      *expt = exp + NUM_LEADING_ZEROS - cnt;
+      return n;
     }
-  else
-    /* Add the implicit leading one bit for a normalized number.  */
-    res_ptr[N - 1] |= (mp_limb_t) 1 << (LDBL_MANT_DIG - 1
-					- ((N - 1) * BITS_PER_MP_LIMB));
-
-  return N;
 }
diff --git a/sysdeps/ieee754/ldbl-96/ldbl2mpn.c b/sysdeps/ieee754/ldbl-96/ldbl2mpn.c
index 3f85e0a..e83171b 100644
--- a/sysdeps/ieee754/ldbl-96/ldbl2mpn.c
+++ b/sysdeps/ieee754/ldbl-96/ldbl2mpn.c
@@ -30,7 +30,7 @@
 
 mp_size_t
 __mpn_extract_long_double (mp_ptr res_ptr, mp_size_t size,
-			   int *expt, int *is_neg,
+			   int *expt, int *zero_bits, int *is_neg,
 			   long double value)
 {
   union ieee854_long_double u;
@@ -52,41 +52,38 @@ __mpn_extract_long_double (mp_ptr res_ptr, mp_size_t size,
   #error "mp_limb size " BITS_PER_MP_LIMB "not accounted for"
 #endif
 
+/* The format does not fill the last limb.  There are some zeros.  */
+#define NUM_LEADING_ZEROS (N * BITS_PER_MP_LIMB - LDBL_MANT_DIG)
+  *zero_bits = NUM_LEADING_ZEROS;
+
   if (u.ieee.exponent == 0)
     {
       /* A biased exponent of zero is a special case.
 	 Either it is a zero or it is a denormal number.  */
       if (res_ptr[0] == 0 && res_ptr[N - 1] == 0) /* Assumes N<=2.  */
-	/* It's zero.  */
-	*expt = 0;
+	{
+	  /* It's zero.  */
+	  *expt = 0;
+	  return 1;
+	}
       else
 	{
           /* It is a denormal number, meaning it has no implicit leading
 	     one bit, and its exponent is in fact the format minimum.  */
 	  int cnt;
+	  int exp = 1 - IEEE854_LONG_DOUBLE_BIAS;
+	  int n = N;
 
 	  if (res_ptr[N - 1] != 0)
-	    {
-	      count_leading_zeros (cnt, res_ptr[N - 1]);
-	      if (cnt != 0)
-		{
-#if N == 2
-	          res_ptr[N - 1] = res_ptr[N - 1] << cnt
-			           | (res_ptr[0] >> (BITS_PER_MP_LIMB - cnt));
-	          res_ptr[0] <<= cnt;
-#else
-	          res_ptr[N - 1] <<= cnt;
-#endif
-		}
-	      *expt = LDBL_MIN_EXP - 1 - cnt;
-	    }
+	    count_leading_zeros (cnt, res_ptr[N - 1]);
 	  else
 	    {
 	      count_leading_zeros (cnt, res_ptr[0]);
-	      res_ptr[N - 1] = res_ptr[0] << cnt;
-	      res_ptr[0] = 0;
-	      *expt = LDBL_MIN_EXP - 1 - BITS_PER_MP_LIMB - cnt;
+	      exp -= BITS_PER_MP_LIMB;
+	      n = 1;
 	    }
+	  *zero_bits = cnt;
+	  *expt = exp + NUM_LEADING_ZEROS - cnt;
 	}
     }
 

-- 
Alan Modra
Australia Development Lab, IBM


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]