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

[FranÃois-Xavier Coudert <Francois-Xavier dot Coudert at lcp dot u-psud dot fr>]

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;