[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 #11 from Dominique d'Humieres <dominiq at lps dot ens.fr> ---
I have opened pr59836 for the wrong outputs in comments 0 and 9 with a patch
which fixes them plus others I have found while debugging. This new patch does
not require the kludge in comment 9. I repost the patch fixing the issue of
this PR, i.e., the inconsistent rounding, without the patch for pr59771:

--- ../_clean/libgfortran/io/write_float.def    2014-01-04 15:51:53.000000000
+0100
+++ libgfortran/io/write_float.def    2014-01-15 23:33:51.000000000 +0100
@@ -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,8 +1044,10 @@ 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))))\
     { \
@@ -1066,10 +1069,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)\
         { \

Besides some cleaning, it computes '0.1 - 0.1 * r * rexp_d' through a volatile
'temp', as it is done later, and computes 'rexp_d * (m + r) >= 1.0' as r >=
(rexp_d - m) (Note that the new rexp_d=10**d and is the inverse of the old
one).

The reason of the second change is that 'r' is 0.0, 0.5, or 1.0, i.e., stored
without rounding, as well as the new rexp_d (at least for small values of 'd',
d<11 for real(4)). So the threshold will trigger only when 'regxp_d-m' is close
to one, i.e., when the difference will lose accuracy and not be affected by
rounding.

I have collected all (most of) the tests in pr48906, pr59771, and pr59836 in a
single file and with the patch there is no differences between the output with
-m32 or -m64. The patch retested cleanly also.