* stdlib/strtod.c: Fixes from drepper. glibc-1.09 cvs/libc_1_09

	[IMPLICIT_ONE]: New macro, one for IEEE754 formats.
	* stdlib/strtold.c [IMPLICIT_ONE]: New macro, 0.
	* stdlib/strtof.c [IMPLICIT_ONE]: New macro, 1.

	* stdio/vfprintf.c: Correctly notice Z modifier.

1 files changed, 91 insertions, 60 deletions

diff --git a/stdlib/strtod.c b/stdlib/strtod.c
index 989984e66d..fa22ae9fe8 100644
--- a/stdlib/strtod.c
+++ b/stdlib/strtod.c
@@ -27,6 +27,7 @@ Cambridge, MA 02139, USA.  */
 #define	STRTOF		strtod
 #define	MPN2FLOAT	__mpn_construct_double
 #define	FLOAT_HUGE_VAL	HUGE_VAL
+#define	IMPLICIT_ONE	1
 #endif
 /* End of configuration part.  */
 
@@ -36,18 +37,19 @@ Cambridge, MA 02139, USA.  */
 #include <localeinfo.h>
 #include <math.h>
 #include <stdlib.h>
-#include "../stdio/gmp.h"
-#include "../stdio/gmp-impl.h"
+#include "gmp.h"
+#include "gmp-impl.h"
 #include <gmp-mparam.h>
-#include "../stdio/longlong.h"
-#include "../stdio/fpioconst.h"
+#include "longlong.h"
+#include "fpioconst.h"
 
-/* #define NDEBUG 1 */
+#define NDEBUG 1
 #include <assert.h>
 
 
 /* Constants we need from float.h; select the set for the FLOAT precision.  */
 #define MANT_DIG	PASTE(FLT,_MANT_DIG)
+#define	DIG		PASTE(FLT,_DIG)
 #define	MAX_EXP		PASTE(FLT,_MAX_EXP)
 #define	MIN_EXP		PASTE(FLT,_MIN_EXP)
 #define MAX_10_EXP	PASTE(FLT,_MAX_10_EXP)
@@ -75,15 +77,16 @@ extern FLOAT MPN2FLOAT (mp_srcptr mpn, int exponent, int negative);
 
 
 /* Local data structure.  */
-static const mp_limb _tens_in_limb[MAX_DIG_PER_LIMB] =
-{    0,                  10,                100,
-     1000,               10000,             100000,
-     1000000,            10000000,          100000000
+static const mp_limb _tens_in_limb[MAX_DIG_PER_LIMB + 1] =
+{    0,                   10,                   100,
+     1000,                10000,                100000,
+     1000000,             10000000,             100000000,
+     1000000000
 #if BITS_PER_MP_LIMB > 32
-   , 1000000000,         10000000000,       100000000000,
-     1000000000000,      10000000000000,    100000000000000,
-     1000000000000000,   10000000000000000, 100000000000000000,
-     1000000000000000000
+	       ,	  10000000000,          100000000000,
+     1000000000000,       10000000000000,       100000000000000,
+     1000000000000000,    10000000000000000,    100000000000000000,
+     1000000000000000000, 10000000000000000000
 #endif
 #if BITS_PER_MP_LIMB > 64
   #error "Need to expand tens_in_limb table to" MAX_DIG_PER_LIMB
@@ -102,7 +105,7 @@ static const mp_limb _tens_in_limb[MAX_DIG_PER_LIMB] =
 	do { if (endptr != 0) *endptr = (char *) end; return val; } while (0)
 
 /* Maximum size necessary for mpn integers to hold floating point numbers.  */ 
-#define	MPNSIZE		(howmany (MAX_EXP + MANT_DIG, BITS_PER_MP_LIMB) + 1)
+#define	MPNSIZE		(howmany (MAX_EXP + 2 * MANT_DIG, BITS_PER_MP_LIMB) + 2)
 /* Declare an mpn integer variable that big.  */
 #define	MPN_VAR(name)	mp_limb name[MPNSIZE]; mp_size_t name##size
 /* Copy an mpn integer value.  */
@@ -116,9 +119,9 @@ static inline FLOAT
 round_and_return (mp_limb *retval, int exponent, int negative,
 		  mp_limb round_limb, mp_size_t round_bit, int more_bits)
 {
-  if (exponent < MIN_EXP)
+  if (exponent < MIN_EXP - 2 + IMPLICIT_ONE)
     {
-      mp_size_t shift = MIN_EXP - 1 - exponent;
+      mp_size_t shift = MIN_EXP - 2 + IMPLICIT_ONE - exponent;
 
       if (shift >= MANT_DIG)
 	{
@@ -131,44 +134,43 @@ round_and_return (mp_limb *retval, int exponent, int negative,
 	{
 	  round_limb = retval[(shift - 1) / BITS_PER_MP_LIMB];
 	  round_bit = (shift - 1) % BITS_PER_MP_LIMB;
-#if RETURN_LIMB_SIZE <= 2
-	  assert (RETURN_LIMB_SIZE == 2);
-	  more_bits |= retval[0] != 0;
-	  retval[0] = retval[1];
-	  retval[1] = 0;
-#else
-	  int disp = shift / BITS_PER_MP_LIMB;
-	  int i = 0;
-	  while (retval[i] == 0 && i < disp)
-	    ++i;
-	  more_bits |= i < disp;
-	  for (i = disp; i < RETURN_LIMB_SIZE; ++i)
-	    retval[i - disp] = retval[i];
-	  MPN_ZERO (&retval[RETURN_LIMB_SIZE - disp], disp);
-#endif
-	  shift %= BITS_PER_MP_LIMB;
+
+	  (void) __mpn_rshift (retval, &retval[shift / BITS_PER_MP_LIMB],
+                               RETURN_LIMB_SIZE - (shift / BITS_PER_MP_LIMB),
+                               shift % BITS_PER_MP_LIMB);
+          MPN_ZERO (&retval[RETURN_LIMB_SIZE - (shift / BITS_PER_MP_LIMB)],
+                    shift / BITS_PER_MP_LIMB);
 	}
-      else
+      else if (shift > 0)
 	{
           round_limb = retval[0];
           round_bit = shift - 1;
+	  (void) __mpn_rshift (retval, retval, RETURN_LIMB_SIZE, shift);
 	}
-      (void) __mpn_rshift (retval, retval, RETURN_LIMB_SIZE, shift);
       exponent = MIN_EXP - 2;
     }
 
-  if ((round_limb & (1 << round_bit)) != 0 &&
-      (more_bits || (retval[0] & 1) != 0 ||
-       (round_limb & ((1 << round_bit) - 1)) != 0))
+  if ((round_limb & (1 << round_bit)) != 0
+      && (more_bits || (retval[0] & 1) != 0
+          || (round_limb & ((1 << round_bit) - 1)) != 0))
     {
       mp_limb cy = __mpn_add_1 (retval, retval, RETURN_LIMB_SIZE, 1);
-      if (cy || (retval[RETURN_LIMB_SIZE - 1]
-		 & (1 << (MANT_DIG % BITS_PER_MP_LIMB))) != 0)
+
+      if (((MANT_DIG % BITS_PER_MP_LIMB) == 0 && cy) ||
+          ((MANT_DIG % BITS_PER_MP_LIMB) != 0 &&
+           (retval[RETURN_LIMB_SIZE - 1]
+            & (1 << (MANT_DIG % BITS_PER_MP_LIMB))) != 0))
 	{
 	  ++exponent;
 	  (void) __mpn_rshift (retval, retval, RETURN_LIMB_SIZE, 1);
-	  retval[RETURN_LIMB_SIZE - 1] |= 1 << (MANT_DIG % BITS_PER_MP_LIMB);
+	  retval[RETURN_LIMB_SIZE - 1] |= 1 << ((MANT_DIG - 1)
+						% BITS_PER_MP_LIMB);
 	}
+      else if (IMPLICIT_ONE && exponent == MIN_EXP - 2
+	       && (retval[RETURN_LIMB_SIZE - 1]
+		   & (1 << ((MANT_DIG - 1) % BITS_PER_MP_LIMB))) != 0)
+	  /* The number was denormalized but now normalized.  */
+	exponent = MIN_EXP - 1;
     }
 
   if (exponent > MAX_EXP)
@@ -305,7 +307,7 @@ STRTOF (nptr, endptr)
   /* Points at the character following the integer and fractional digits.  */
   const char *expp;
   /* Total number of digit and number of digits in integer part.  */
-  int dig_no, int_no;
+  int dig_no, int_no, lead_zero;
   /* Contains the last character read.  */
   char c;
 
@@ -440,15 +442,18 @@ STRTOF (nptr, endptr)
   /* We have the number digits in the integer part.  Whether these are all or
      any is really a fractional digit will be decided later.  */
   int_no = dig_no;
+  lead_zero = int_no == 0 ? -1 : 0;
 
   /* Read the fractional digits.  */
   if (c == decimal)
     {
       if (isdigit (cp[1]))
 	{
-	  ++cp;
+	  c = *++cp;
 	  do
 	    {
+	      if (c != '0' && lead_zero == -1)
+		lead_zero = dig_no - int_no;
 	      ++dig_no;
 	      c = *++cp;
 	    }
@@ -547,7 +552,7 @@ STRTOF (nptr, endptr)
       return negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL;
     }
 
-  if (int_no - dig_no + exponent < MIN_10_EXP - MANT_DIG)
+  if (exponent - MAX(0, lead_zero) < MIN_10_EXP - (DIG + 1))
     {
       errno = ERANGE;
       return 0.0;
@@ -607,20 +612,26 @@ STRTOF (nptr, endptr)
 	 to determine the result.  */
       if (bits > MANT_DIG)
 	{
+	  int i;
 	  const mp_size_t least_idx = (bits - MANT_DIG) / BITS_PER_MP_LIMB;
 	  const mp_size_t least_bit = (bits - MANT_DIG) % BITS_PER_MP_LIMB;
 	  const mp_size_t round_idx = least_bit == 0 ? least_idx - 1
 						     : least_idx;
 	  const mp_size_t round_bit = least_bit == 0 ? BITS_PER_MP_LIMB - 1
-						     : least_idx - 1;
-	  int i;
+						     : least_bit - 1;
 
 	  if (least_bit == 0)
 	    memcpy (retval, &num[least_idx],
 		    RETURN_LIMB_SIZE * sizeof (mp_limb));
 	  else
-	    (void) __mpn_rshift (retval, &num[least_idx],
-				 numsize - least_idx + 1, least_bit);
+            {
+              for (i = least_idx; i < numsize - 1; ++i)
+                retval[i - least_idx] = (num[i] >> least_bit)
+                                        | (num[i + 1]
+                                           << (BITS_PER_MP_LIMB - least_bit));
+              if (i - least_idx < RETURN_LIMB_SIZE)
+                retval[RETURN_LIMB_SIZE - 1] = num[i] >> least_bit;
+            }
 
 	  /* Check whether any limb beside the ones in RETVAL are non-zero.  */
 	  for (i = 0; num[i] == 0; ++i)
@@ -686,14 +697,32 @@ STRTOF (nptr, endptr)
 
     int expbit;
     int cnt;
+    int neg_exp;
+    int more_bits;
     mp_limb cy;
     mp_limb *psrc = den;
     mp_limb *pdest = num;
-    int neg_exp = dig_no - int_no - exponent;
     const struct mp_power *ttab = &_fpioconst_pow10[0];
 
     assert (dig_no > int_no && exponent <= 0);
 
+
+    /* For the fractional part we need not process too much digits.  One
+       decimal digits gives us log_2(10) ~ 3.32 bits.  If we now compute
+                        ceil(BITS / 3) =: N
+       digits we should have enough bits for the result.  The remaining
+       decimal digits give us the information that more bits are following.
+       This can be used while rounding.  (One added as a safety margin.)  */
+    if (dig_no - int_no > (MANT_DIG - bits + 2) / 3 + 1)
+      {
+        dig_no = int_no + (MANT_DIG - bits + 2) / 3 + 1;
+        more_bits = 1;
+      }
+    else
+      more_bits = 0;
+
+    neg_exp = dig_no - int_no - exponent;
+
     /* Construct the denominator.  */
     densize = 0;
     expbit = 1;
@@ -782,36 +811,38 @@ STRTOF (nptr, endptr)
 		  exponent -= cnt;					      \
 		  if (BITS_PER_MP_LIMB - cnt > MANT_DIG)		      \
 		    {							      \
-		      used = cnt + MANT_DIG;				      \
+		      used = MANT_DIG + cnt;				      \
 		      retval[0] = quot >> (BITS_PER_MP_LIMB - used);	      \
-		      bits -= BITS_PER_MP_LIMB - used;			      \
+		      bits = MANT_DIG + 1;				      \
 		    }							      \
 		  else							      \
 		    {							      \
 		      /* Note that we only clear the second element.  */      \
-		      retval[1] = 0;					      \
+		      /* The conditional is determined at compile time.  */   \
+		      if (RETURN_LIMB_SIZE > 1)				      \
+			retval[1] = 0;					      \
 		      retval[0] = quot;					      \
-		      bits -= cnt;					      \
+		      bits = -cnt;					      \
 		    }							      \
 		}							      \
 	      else if (bits + BITS_PER_MP_LIMB <= MANT_DIG)		      \
-		__mpn_lshift_1 (retval, RETURN_LIMB_SIZE, BITS_PER_MP_LIMB,    \
+		__mpn_lshift_1 (retval, RETURN_LIMB_SIZE, BITS_PER_MP_LIMB,   \
 				quot);					      \
 	      else							      \
 		{							      \
 		  used = MANT_DIG - bits;				      \
 		  if (used > 0)						      \
-		    __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, quot);     \
+		    __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, quot);    \
 		}							      \
 	      bits += BITS_PER_MP_LIMB
 
-              got_limb;
+	      got_limb;
 	    }
 	  while (bits <= MANT_DIG);
 
 	  return round_and_return (retval, exponent - 1, negative,
 				   quot, BITS_PER_MP_LIMB - 1 - used,
-				   n != 0);
+				   more_bits || n != 0);
 	}
       case 2:
 	{
@@ -893,7 +924,7 @@ STRTOF (nptr, endptr)
 	    
 	  return round_and_return (retval, exponent - 1, negative,
 				   quot, BITS_PER_MP_LIMB - 1 - used,
-				   n1 != 0 || n0 != 0);
+				   more_bits || n1 != 0 || n0 != 0);
 	}
       default:
 	{
@@ -907,7 +938,7 @@ STRTOF (nptr, endptr)
 
 	  /* The division does not work if the upper limb of the two-limb
 	     numerator is greater than the denominator.  */
-	  if (num[numsize - 1] > dX)
+	  if (__mpn_cmp (num, &den[densize - numsize], numsize) > 0)
 	    num[numsize++] = 0;
 
 	  if (numsize < densize)
@@ -1017,11 +1048,11 @@ STRTOF (nptr, endptr)
 	      got_limb;
 	    }
 
-	  for (i = densize - 1; num[i] != 0 && i >= 0; --i)
+	  for (i = densize; num[i] == 0 && i >= 0; --i)
 	    ;
 	  return round_and_return (retval, exponent - 1, negative,
 				   quot, BITS_PER_MP_LIMB - 1 - used,
-				   i >= 0);
+				   more_bits || i >= 0);
 	}
       }
   }