Re: real_to_decimal rounds incorrectly

[Richard Henderson <rth@redhat.com>]

I've forward-ported the logic from the previous real.c, and added
some commentary so that the next person forced to maintain this
code won't have to spend quite so long trying to figure out what's
going on.

If you have a literature reference to the algorithm being used here,
it would be nice to add a pointer to it.



r~


        * real.c (cmp_significand_0, rtd_divmod, ten_to_mptwo): New.
        (real_to_decimal): Re-implement using the logic from the
        gcc 3.2 etoasc.  Comment heavily.
        (div_significands): Simplify loop startup and comparison logic.

Index: real.c
===================================================================
RCS file: /cvs/gcc/gcc/gcc/real.c,v
retrieving revision 1.98
diff -c -p -d -b -r1.98 real.c
*** real.c	16 Oct 2002 00:40:27 -0000	1.98
--- real.c	18 Oct 2002 23:37:39 -0000
*************** static void neg_significand PARAMS ((REA
*** 102,107 ****
--- 102,108 ----
  				     const REAL_VALUE_TYPE *));
  static int cmp_significands PARAMS ((const REAL_VALUE_TYPE *,
  				     const REAL_VALUE_TYPE *));
+ static int cmp_significand_0 PARAMS ((const REAL_VALUE_TYPE *));
  static void set_significand_bit PARAMS ((REAL_VALUE_TYPE *, unsigned int));
  static void clear_significand_bit PARAMS ((REAL_VALUE_TYPE *, unsigned int));
  static bool test_significand_bit PARAMS ((REAL_VALUE_TYPE *, unsigned int));
*************** static void do_divide PARAMS ((REAL_VALU
*** 121,130 ****
  			       const REAL_VALUE_TYPE *));
  static int do_compare PARAMS ((const REAL_VALUE_TYPE *,
  			       const REAL_VALUE_TYPE *, int));
! static void do_fix_trunc PARAMS ((REAL_VALUE_TYPE *,
! 				  const REAL_VALUE_TYPE *));
  
  static const REAL_VALUE_TYPE * ten_to_ptwo PARAMS ((int));
  static const REAL_VALUE_TYPE * real_digit PARAMS ((int));
  static void times_pten PARAMS ((REAL_VALUE_TYPE *, int));
  
--- 122,134 ----
  			       const REAL_VALUE_TYPE *));
  static int do_compare PARAMS ((const REAL_VALUE_TYPE *,
  			       const REAL_VALUE_TYPE *, int));
! static void do_fix_trunc PARAMS ((REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *));
! 
! static unsigned long rtd_divmod PARAMS ((REAL_VALUE_TYPE *,
! 					 REAL_VALUE_TYPE *));
  
  static const REAL_VALUE_TYPE * ten_to_ptwo PARAMS ((int));
+ static const REAL_VALUE_TYPE * ten_to_mptwo PARAMS ((int));
  static const REAL_VALUE_TYPE * real_digit PARAMS ((int));
  static void times_pten PARAMS ((REAL_VALUE_TYPE *, int));
  
*************** cmp_significands (a, b)
*** 406,411 ****
--- 410,430 ----
    return 0;
  }
  
+ /* Return true if A is non-zero.  */
+ 
+ static inline int 
+ cmp_significand_0 (a)
+      const REAL_VALUE_TYPE *a;
+ {
+   int i;
+ 
+   for (i = SIGSZ - 1; i >= 0; --i)
+     if (a->sig[i])
+       return 1;
+ 
+   return 0;
+ }
+ 
  /* Set bit N of the significand of R.  */
  
  static inline void
*************** div_significands (r, a, b)
*** 466,500 ****
       const REAL_VALUE_TYPE *a, *b;
  {
    REAL_VALUE_TYPE u;
!   int bit = SIGNIFICAND_BITS - 1;
!   int i;
!   long inexact;
  
    u = *a;
    memset (r->sig, 0, sizeof (r->sig));
  
    goto start;
    do
      {
!       if ((u.sig[SIGSZ-1] & SIG_MSB) == 0)
! 	{
  	  lshift_significand_1 (&u, &u);
  	start:
! 	  if (cmp_significands (&u, b) >= 0)
  	    {
  	      sub_significands (&u, &u, b);
  	      set_significand_bit (r, bit);
  	    }
  	}
-       else
- 	{
- 	  /* We lose a bit here, and thus know the next quotient bit
- 	     will be one.  */
- 	  lshift_significand_1 (&u, &u);
- 	  sub_significands (&u, &u, b);
- 	  set_significand_bit (r, bit);
- 	}
-     }
    while (--bit >= 0);
  
    for (i = 0, inexact = 0; i < SIGSZ; i++)
--- 485,509 ----
       const REAL_VALUE_TYPE *a, *b;
  {
    REAL_VALUE_TYPE u;
!   int i, bit = SIGNIFICAND_BITS - 1;
!   unsigned long msb, inexact;
  
    u = *a;
    memset (r->sig, 0, sizeof (r->sig));
  
+   msb = 0;
    goto start;
    do
      {
!       msb = u.sig[SIGSZ-1] & SIG_MSB;
        lshift_significand_1 (&u, &u);
      start:
!       if (msb || cmp_significands (&u, b) >= 0)
  	{
  	  sub_significands (&u, &u, b);
  	  set_significand_bit (r, bit);
  	}
      }
    while (--bit >= 0);
  
    for (i = 0, inexact = 0; i < SIGSZ; i++)
*************** do_fix_trunc (r, a)
*** 972,978 ****
  {
    *r = *a;
  
!   switch (a->class)
      {
      case rvc_zero:
      case rvc_inf:
--- 981,987 ----
  {
    *r = *a;
  
!   switch (r->class)
      {
      case rvc_zero:
      case rvc_inf:
*************** real_to_integer2 (plow, phigh, r)
*** 1396,1401 ****
--- 1405,1447 ----
    *phigh = high;
  }
  
+ /* A subroutine of real_to_decimal.  Compute the quotient and remainder
+    of NUM / DEN.  Return the quotient and place the remainder in NUM.
+    It is expected that NUM / DEN are close enough that the quotient is
+    small.  */
+ 
+ static unsigned long
+ rtd_divmod (num, den)
+      REAL_VALUE_TYPE *num, *den;
+ {
+   unsigned long q, msb;
+   int expn = num->exp, expd = den->exp;
+ 
+   if (expn < expd)
+     return 0;
+ 
+   q = msb = 0;
+   goto start;
+   do
+     {
+       msb = num->sig[SIGSZ-1] & SIG_MSB;
+       q <<= 1;
+       lshift_significand_1 (num, num);
+     start:
+       if (msb || cmp_significands (num, den) >= 0)
+ 	{
+ 	  sub_significands (num, num, den);
+ 	  q |= 1;
+ 	}
+     }
+   while (--expn >= expd);
+ 
+   num->exp = expd;
+   normalize (num);
+ 
+   return q;
+ }
+ 
  /* Render R as a decimal floating point constant.  Emit DIGITS significant
     digits in the result, bounded by BUF_SIZE.  If DIGITS is 0, choose the
     maximum for the representation.  If CROP_TRAILING_ZEROS, strip trailing
*************** real_to_decimal (str, r_orig, buf_size, 
*** 1410,1421 ****
       size_t buf_size, digits;
       int crop_trailing_zeros;
  {
-   REAL_VALUE_TYPE r;
    const REAL_VALUE_TYPE *one, *ten;
!   int dec_exp, d, cmp_half;
    size_t max_digits;
    char *p, *first, *last;
-   char exp_buf[16];
    bool sign;
  
    r = *r_orig;
--- 1456,1466 ----
       size_t buf_size, digits;
       int crop_trailing_zeros;
  {
    const REAL_VALUE_TYPE *one, *ten;
!   REAL_VALUE_TYPE r, pten, u, v;
!   int dec_exp, cmp_one, digit;
    size_t max_digits;
    char *p, *first, *last;
    bool sign;
  
    r = *r_orig;
*************** real_to_decimal (str, r_orig, buf_size, 
*** 1437,1465 ****
        abort ();
      }
  
    one = real_digit (1);
    ten = ten_to_ptwo (0);
  
    sign = r.sign;
    r.sign = 0;
  
!   /* Estimate the decimal exponent.  */
!   dec_exp = r.exp * M_LOG10_2;
    
!   /* Scale the number such that it is in [1, 10).  */
!   times_pten (&r, (dec_exp > 0 ? -dec_exp : -(--dec_exp)));
  
!   /* Assert that the number is in the proper range.  Round-off can
!      prevent the above from working exactly.  */
!   if (do_compare (&r, one, -1) < 0)
      {
!       do_multiply (&r, &r, ten);
!       dec_exp--;
      }
!   else if (do_compare (&r, ten, 1) >= 0)
      {
!       do_divide (&r, &r, ten);
!       dec_exp++;
      }
  
    p = str;
--- 1482,1612 ----
        abort ();
      }
  
+   /* Estimate the decimal exponent, and compute the length of the string it
+      will print as.  Be conservative and add one to account for possible
+      overflow or rounding error.  */
+   dec_exp = r.exp * M_LOG10_2;
+   for (max_digits = 1; dec_exp ; max_digits++)
+     dec_exp /= 10;
+ 
+   /* Bound the number of digits printed by the size of the output buffer.  */
+   max_digits = buf_size - 1 - 1 - 2 - max_digits - 1;
+   if (max_digits > buf_size)
+     abort ();
+   if (digits > max_digits)
+     digits = max_digits;
+ 
+   /* Bound the number of digits printed by the size of the representation.  */
+   max_digits = SIGNIFICAND_BITS * M_LOG10_2;
+   if (digits == 0 || digits > max_digits)
+     digits = max_digits;
+ 
    one = real_digit (1);
    ten = ten_to_ptwo (0);
  
    sign = r.sign;
    r.sign = 0;
  
!   dec_exp = 0;
!   pten = *one;
  
!   cmp_one = do_compare (&r, one, 0);
!   if (cmp_one > 0)
!     {
!       int m;
  
!       /* Number is greater than one.  Convert significand to an integer
! 	 and strip trailing decimal zeros.  */
! 
!       u = r;
!       u.exp = SIGNIFICAND_BITS - 1;
! 
!       /* Largest M, such that 10**2**M fits within SIGNIFICAND_BITS.  */
!       m = floor_log2 (max_digits);
! 
!       /* Iterate over the bits of the possible powers of 10 that might
! 	 be present in U and eliminate them.  That is, if we find that
! 	 10**2**M divides U evenly, keep the division and increase 
! 	 DEC_EXP by 2**M.  */
!       do
  	{
! 	  REAL_VALUE_TYPE t;
! 
! 	  do_divide (&t, &u, ten_to_ptwo (m));
! 	  do_fix_trunc (&v, &t);
! 	  if (cmp_significands (&v, &t) == 0)
! 	    {
! 	      u = t;
! 	      dec_exp += 1 << m;
  	    }
! 	}
!       while (--m >= 0);
! 
!       /* Revert the scaling to integer that we performed earlier.  */
!       u.exp += r.exp - (SIGNIFICAND_BITS - 1);
!       r = u;
! 
!       /* Find power of 10.  Do this by dividing out 10**2**M when
! 	 this is larger than the current remainder.  Fill PTEN with 
! 	 the power of 10 that we compute.  */
!       m = floor_log2 ((int)(r.exp * M_LOG10_2)) + 1;
!       do
  	{
! 	  const REAL_VALUE_TYPE *ptentwo = ten_to_ptwo (m);
! 	  if (do_compare (&u, ptentwo, 0) >= 0)
! 	    {
! 	      do_divide (&u, &u, ptentwo);
! 	      do_multiply (&pten, &pten, ptentwo);
! 	      dec_exp += 1 << m;
! 	    }
! 	}
!       while (--m >= 0);
!     }
!   else if (cmp_one < 0)
!     {
!       int m;
! 
!       /* Number is less than one.  Pad significand with leading
! 	 decimal zeros.  */
! 
!       v = r;
!       while (1)
! 	{
! 	  /* Stop if we'd shift bits off the bottom.  */
! 	  if (v.sig[0] & 7)
! 	    break;
! 
! 	  do_multiply (&u, &v, ten);
! 
! 	  /* Stop if we're now >= 1.  */
! 	  if (u.exp > 0)
! 	    break;
! 
! 	  v = u;
! 	  dec_exp -= 1;
! 	}
!       r = v;
! 
!       /* Find power of 10.  Do this by multiplying in P=10**2**M when
! 	 the current remainder is smaller than 1/P.  Fill PTEN with the
! 	 power of 10 that we compute.  */
!       m = floor_log2 ((int)(-r.exp * M_LOG10_2)) + 1;
!       do
! 	{
! 	  const REAL_VALUE_TYPE *ptentwo = ten_to_ptwo (m);
! 	  const REAL_VALUE_TYPE *ptenmtwo = ten_to_mptwo (m);
! 
! 	  if (do_compare (&v, ptenmtwo, 0) <= 0)
! 	    {
! 	      do_multiply (&v, &v, ptentwo);
! 	      do_multiply (&pten, &pten, ptentwo);
! 	      dec_exp -= 1 << m;
! 	    }
! 	}
!       while (--m >= 0);
! 
!       /* Invert the positive power of 10 that we've collected so far.  */
!       do_divide (&pten, one, &pten);
      }
  
    p = str;
*************** real_to_decimal (str, r_orig, buf_size, 
*** 1467,1518 ****
      *p++ = '-';
    first = p++;
  
!   sprintf (exp_buf, "e%+d", dec_exp);
  
!   /* Bound the number of digits printed by the size of the representation.  */
!   max_digits = SIGNIFICAND_BITS * M_LOG10_2;
!   if (digits == 0 || digits > max_digits)
!     digits = max_digits;
  
!   /* Bound the number of digits printed by the size of the output buffer.  */
!   max_digits = buf_size - strlen (exp_buf) - sign - 1;
!   if (max_digits > buf_size)
      abort ();
!   if (digits > max_digits)
!     digits = max_digits;
  
!   while (1)
      {
!       d = real_to_integer ((const REAL_VALUE_TYPE *) &r);
!       do_add (&r, &r, real_digit (d), 1);
  
!       *p++ = d + '0';
!       if (--digits == 0)
! 	break;
        do_multiply (&r, &r, ten);
      }
    last = p;
  
!   /* Round the result.  Compare R vs 0.5 by doing R*2 vs 1.0.  */
!   r.exp += 1;
!   cmp_half = do_compare (&r, one, -1);
!   if (cmp_half == 0)
      /* Round to even.  */
!     cmp_half += d & 1;
!   if (cmp_half > 0)
      {
        while (p > first)
  	{
! 	  d = *--p;
! 	  if (d == '9')
  	    *p = '0';
  	  else
  	    {
! 	      *p = d + 1;
  	      break;
  	    }
  	}
  
        if (p == first)
  	{
  	  first[1] = '1';
--- 1614,1693 ----
      *p++ = '-';
    first = p++;
  
!   /* At this point, PTEN should contain the nearest power of 10 smaller
!      than R, such that this division produces the first digit.
  
!      Using a divide-step primitive that returns the complete integral
!      remainder avoids the rounding error that would be produced if
!      we were to use do_divide here and then simply multiply by 10 for
!      each subsequent digit.  */
  
!   digit = rtd_divmod (&r, &pten);
! 
!   /* Be prepared for error in that division via underflow ... */
!   if (digit == 0 && cmp_significand_0 (&r))
!     {
!       /* Multiply by 10 and try again.  */
!       do_multiply (&r, &r, ten);
!       digit = rtd_divmod (&r, &pten);
!       dec_exp -= 1;
!       if (digit == 0)
  	abort ();
!     }
  
!   /* ... or overflow.  */
!   if (digit == 10)
      {
!       *p++ = '1';
!       if (--digits > 0)
! 	*p++ = '0';
!       dec_exp += 1;
!     }
!   else if (digit > 10)
!     abort ();
!   else
!     *p++ = digit + '0';
  
!   /* Generate subsequent digits.  */
!   while (--digits > 0)
!     {
        do_multiply (&r, &r, ten);
+       digit = rtd_divmod (&r, &pten);
+       *p++ = digit + '0';
      }
    last = p;
  
!   /* Generate one more digit with which to do rounding.  */
!   do_multiply (&r, &r, ten);
!   digit = rtd_divmod (&r, &pten);
! 
!   /* Round the result.  */
!   if (digit == 5)
!     {
!       /* Round to nearest.  If R is non-zero there are additional
! 	 non-zero digits to be extracted.  */
!       if (cmp_significand_0 (&r))
! 	digit++;
        /* Round to even.  */
!       else if ((p[-1] - '0') & 1)
! 	digit++;
!     }
!   if (digit > 5)
      {
        while (p > first)
  	{
! 	  digit = *--p;
! 	  if (digit == '9')
  	    *p = '0';
  	  else
  	    {
! 	      *p = digit + 1;
  	      break;
  	    }
  	}
  
+       /* Carry out of the first digit.  This means we had all 9's and
+ 	 now have all 0's.  "Prepend" a 1 by overwriting the first 0.  */
        if (p == first)
  	{
  	  first[1] = '1';
*************** real_to_decimal (str, r_orig, buf_size, 
*** 1520,1533 ****
  	}
      }
    
    first[0] = first[1];
    first[1] = '.';
  
    if (crop_trailing_zeros)
      while (last > first + 3 && last[-1] == '0')
        last--;
  
!   strcpy (last, exp_buf);
  }
  
  /* Render R as a hexadecimal floating point constant.  Emit DIGITS
--- 1695,1711 ----
  	}
      }
    
+   /* Insert the decimal point.  */
    first[0] = first[1];
    first[1] = '.';
  
+   /* If requested, drop trailing zeros.  Never crop past "1.0".  */
    if (crop_trailing_zeros)
      while (last > first + 3 && last[-1] == '0')
        last--;
  
!   /* Append the exponent.  */
!   sprintf (last, "e%+d", dec_exp);
  }
  
  /* Render R as a hexadecimal floating point constant.  Emit DIGITS
*************** real_from_integer (r, mode, low, high, u
*** 1857,1863 ****
      real_convert (r, mode, r);
  }
  
! /* Returns 10**2**n.  */
  
  static const REAL_VALUE_TYPE *
  ten_to_ptwo (n)
--- 2035,2041 ----
      real_convert (r, mode, r);
  }
  
! /* Returns 10**2**N.  */
  
  static const REAL_VALUE_TYPE *
  ten_to_ptwo (n)
*************** ten_to_ptwo (n)
*** 1890,1895 ****
--- 2068,2090 ----
    return &tens[n];
  }
  
+ /* Returns 10**(-2**N).  */
+ 
+ static const REAL_VALUE_TYPE *
+ ten_to_mptwo (n)
+      int n;
+ {
+   static REAL_VALUE_TYPE tens[EXP_BITS];
+ 
+   if (n < 0 || n >= EXP_BITS)
+     abort ();
+ 
+   if (tens[n].class == rvc_zero)
+     do_divide (&tens[n], real_digit (1), ten_to_ptwo (n));
+ 
+   return &tens[n];
+ }
+ 
  /* Returns N.  */
  
  static const REAL_VALUE_TYPE *