[Bug libfortran/59774] [Regression] Inconsistent rounding between -m32 and -m64

["dominiq at lps dot ens.fr" <gcc-bugzilla@gcc.gnu.org>]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=59774

--- Comment #9 from Dominique d'Humieres <dominiq at lps dot ens.fr> ---
I have understood the problem in comment 8. It is illustrated by the following
code

 print "(ru,g45.3)", 891.1
 print "(rd,g45.3)", -891.1
end

which gives the output

                                       9.    
                                      -9.    

with current releases and trunk. The problem comes from the lines

  newf.u.real.d = m == 0.0 ? d - 1 : -(mid - d - 1) ;\

and

  if (w > 0 && d == 0 && p == 0)
    nbefore = 1;

in libgfortran/io/write_float.def when mid==d+1.

I have also noticed the sentence "the asm volatile is required for 32-bit x86
platforms" which seems to answer my question in comment 6. These remarks led me
to the following patch

--- ../_clean/libgfortran/io/write_float.def    2014-01-04 15:51:53.000000000
+0100
+++ libgfortran/io/write_float.def    2014-01-14 22:55:24.000000000 +0100
@@ -112,7 +112,7 @@ determine_precision (st_parameter_dt * d

 static bool
 output_float (st_parameter_dt *dtp, const fnode *f, char *buffer, size_t size,
-          int nprinted, int precision, int sign_bit, bool zero_flag)
+          int nprinted, int precision, int sign_bit, bool zero_flag, int d_o)
 {
   char *out;
   char *digits;
@@ -373,7 +373,7 @@ output_float (st_parameter_dt *dtp, cons
   updown:

   rchar = '0';
-  if (w > 0 && d == 0 && p == 0)
+  if (w > 0 && d_o == 0 && p == 0)
     nbefore = 1;
   /* Scan for trailing zeros to see if we really need to round it.  */
   for(i = nbefore + nafter; i < ndigits; i++)
@@ -1018,13 +1018,14 @@ output_float_FMT_G_ ## x (st_parameter_d
   int d = f->u.real.d;\
   int w = f->u.real.w;\
   fnode newf;\
-  GFC_REAL_ ## x rexp_d, r = 0.5;\
+  GFC_REAL_ ## x rexp_d, r = 0.5, r_sc;\
   int low, high, mid;\
   int ubound, lbound;\
   char *p, pad = ' ';\
   int save_scale_factor, nb = 0;\
   bool result;\
   int nprinted, precision;\
+  volatile GFC_REAL_ ## x temp;\
 \
   save_scale_factor = dtp->u.p.scale_factor;\
 \
@@ -1043,10 +1044,13 @@ output_float_FMT_G_ ## x (st_parameter_d
     break;\
     }\
 \
-  rexp_d = calculate_exp_ ## x (-d);\
-  if ((m > 0.0 && ((m < 0.1 - 0.1 * r * rexp_d) || (rexp_d * (m + r) >=
1.0)))\
+  rexp_d = calculate_exp_ ## x (d);\
+  r_sc = (1 - r / rexp_d);\
+  temp = 0.1 * r_sc;\
+  if ((m > 0.0 && ((m < temp) || (r >= (rexp_d - m))))\
       || ((m == 0.0) && !(compile_options.allow_std\
-              & (GFC_STD_F2003 | GFC_STD_F2008))))\
+              & (GFC_STD_F2003 | GFC_STD_F2008)))\
+      ||  d == 0)\
     { \
       newf.format = FMT_E;\
       newf.u.real.w = w;\
@@ -1066,10 +1070,9 @@ output_float_FMT_G_ ## x (st_parameter_d
 \
   while (low <= high)\
     { \
-      volatile GFC_REAL_ ## x temp;\
       mid = (low + high) / 2;\
 \
-      temp = (calculate_exp_ ## x (mid - 1) * (1 - r * rexp_d));\
+      temp = (calculate_exp_ ## x (mid - 1) * r_sc);\
 \
       if (m < temp)\
         { \
@@ -1106,7 +1109,7 @@ output_float_FMT_G_ ## x (st_parameter_d
 \
  finish:\
     result = output_float (dtp, &newf, buffer, size, nprinted, precision,\
-               sign_bit, zero_flag);\
+               sign_bit, zero_flag, d);\
   dtp->u.p.scale_factor = save_scale_factor;\
 \
 \
@@ -1240,7 +1243,7 @@ determine_en_precision (st_parameter_dt 
         else\
           nprinted = DTOA(y,precision,tmp);                    \
         output_float (dtp, f, buffer, size, nprinted, precision,\
-              sign_bit, zero_flag);\
+              sign_bit, zero_flag, f->u.real.d);\
       }\
 }\


I agree that the additional dummy argument d_o is a kludge, but I did not find
a better way to distinguish between d==0 in the format and mid==d+1. Comments
and improvements welcomed.

Regtested without regression r206590 plus the patch.