real_to_decimal rounds incorrectly

real_to_decimal rounds incorrectly

[Stephen L Moshier]

/* Link this test case with real.o and libiberty.a:
  ./xgcc -B./ tst.c real.o ../libiberty/libiberty.a

  (The test case assumes 'long' is 32 bits, for real_from_target_fmt.)   */

#include <stdio.h>
extern const struct real_format ieee_double_format;
extern void real_from_target_fmt (unsigned long *r, const long *buf,
					     const struct real_format
					     *fmt);
extern void real_to_decimal (char *str, const unsigned long *r_orig,
			         int digits);
void fancy_abort () {}

double td = 1.975949500000000000000000E6;

main ()
{
  char s[80];
  unsigned long r[6];
  union
  {
    unsigned long l[4];
    double d;
  } u;

  u.d = td;
  printf ("Compiled double value %.16e\n", u.d);
  real_from_target_fmt (r, &u.l[0], &ieee_double_format);
  real_to_decimal (s, r, 7);
  printf ("real_to_decimal value %s, should be 1.975950e6\n", s);
}

Re: real_to_decimal rounds incorrectly

[Richard Henderson]

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 *

Re: real_to_decimal rounds incorrectly

[Stephen L Moshier]

> 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.

IEEE-specific modifications were taken from the IEEE spec.
Other than that, I draw a blank on references.  It's just arithmetic.
Maybe Brad, the Purdue Profossor, can suggest some reading material.

Purdue is a world hotbed of numerical analysis work, by the way.  It
would be wonderful if you could talk them into mounting a GNU support
project.  But I'm afraid that mathematicians usually consider this
stuff too trivial and uninteresting to spend much time on.

Re: real_to_decimal rounds incorrectly

[Richard Henderson]

On Sat, Oct 19, 2002 at 01:54:47PM -0400, Stephen L Moshier wrote:
> IEEE-specific modifications were taken from the IEEE spec.
> Other than that, I draw a blank on references.  It's just arithmetic.

By inference, it appears to be (a variation of) Steele & White's 1990
algorithm.

I havn't been able to get hold of the original paper, not being an
ACM member, but I've found a few papers that expand on Steele & White's,
and which give a description good enough to recognize the salient points.


r~

Re: real_to_decimal rounds incorrectly

[Stephen L Moshier]

On Sat, 19 Oct 2002, Richard Henderson wrote:

> On Sat, Oct 19, 2002 at 01:54:47PM -0400, Stephen L Moshier wrote:
> > IEEE-specific modifications were taken from the IEEE spec.
> > Other than that, I draw a blank on references.  It's just arithmetic.
>
> By inference, it appears to be (a variation of) Steele & White's 1990
> algorithm.
>
> I havn't been able to get hold of the original paper, not being an
> ACM member, but I've found a few papers that expand on Steele & White's,
> and which give a description good enough to recognize the salient points.


I really can't say.  The program was written in 1978.  The remarks in
the standard and the thesis of Coonen were helpful for the
modifcations that were made when the IEEE standard came out.
Those are the suggestions like padding zeros to minimize the power of 10.
I've just been looking through the books I have on numerical
methods, but there seems to be nothing at all in them about
radix conversion techniques.  Maybe Knuth's book has something.

Re: real_to_decimal rounds incorrectly

[Stephen L Moshier]

Yesterday I succeeded in getting an old test for real.c running on a
sparcstation, in order to check the 128-bit quad precision decimal
conversions.  They seemed accurate enough to satisfy IEEE, which
as far as I know is the policy we are supposed to meet.
The arithmetic operations, in standard precisions, were already
shown to be probably OK.  So I do not have any new bugs to
report.  I do have more old tests, maybe will get to them next weekend.