[PATCH] gfortran: PR libfortran/20124 - incorrect rounding when writing floating point numbers

[François-Xavier Coudert]

This patch, *strapped and regtested on i686-linux, fixes PR 20124.

Some details: I would have liked to be able not to use another variable 
(which I called nzero_real), but I couldn't find any simple other way to do 
it. So, we store in nzero_real the number of leading zeros in the 
fractionnal part of the floating point number, whatever the width of the 
format is.

I hereby donate this patch to the FSF.


2005-02-24  François-Xavier Coudert  <coudert@clipper.ens.fr>

     PR libfortran/20124
     * write.c (output_float): Adds a nzero_real variable to store the 
number of leading zeros whatever the format width is.  Corrects the rounding 
of numbers less than 10^(-width).  Fixes typo in an error message.  Updates 
copyright years.


Index: libgfortran/io/write.c
===================================================================
RCS file: /cvsroot/gcc/gcc/libgfortran/io/write.c,v
retrieving revision 1.22
diff -u -p -r1.22 write.c
--- libgfortran/io/write.c      12 Jan 2005 21:27:31 -0000      1.22
+++ libgfortran/io/write.c      24 Feb 2005 17:33:05 -0000
@@ -1,4 +1,4 @@
-/* Copyright (C) 2002, 2003, 2004 Free Software Foundation, Inc.
+/* Copyright (C) 2002, 2003, 2004, 2005 Free Software Foundation, Inc.
     Contributed by Andy Vaught

  This file is part of the GNU Fortran 95 runtime library (libgfortran).
@@ -286,6 +286,8 @@ output_float (fnode *f, double value, in
    int nzero;
    /* Number of digits after the decimal point.  */
    int nafter;
+  /* Number of zeros after the decimal point, whatever the precision.  */
+  int nzero_real;
    int leadzero;
    int nblanks;
    int i;
@@ -295,9 +297,11 @@ output_float (fnode *f, double value, in
    w = f->u.real.w;
    d = f->u.real.d;

+  nzero_real = -1;
+
    /* We should always know the field width and precision.  */
    if (d < 0)
-    internal_error ("Uspecified precision");
+    internal_error ("Unspecified precision");

    /* Use sprintf to print the number in the format +D.DDDDe+ddd
       For an N digit exponent, this gives us (32-6)-N digits after the
@@ -359,6 +363,7 @@ output_float (fnode *f, double value, in
        if (nbefore < 0)
         {
           nzero = -nbefore;
+         nzero_real = nzero;
           if (nzero > d)
             nzero = d;
           nafter = d - nzero;
@@ -435,7 +440,17 @@ output_float (fnode *f, double value, in

    /* Round the value.  */
    if (nbefore + nafter == 0)
+  {
      ndigits = 0;
+    if (nzero_real == d && digits[0] >= '5')
+      {
+       /* We rounded to zero but shouldn't have */
+       nzero--;
+       nafter = 1;
+       digits[0] = '1';
+       ndigits = 1;
+      }
+  }
    else if (nbefore + nafter < ndigits)
      {
        ndigits = nbefore + nafter;