public inbox for gcc-patches@gcc.gnu.org
 help / color / mirror / Atom feed
* [patch, libgfortran] Wrong rounding results with -m32
@ 2011-06-03 12:51 jerry DeLisle
  2011-06-10 17:39 ` [patch, libgfortran] PR48906 " jerry DeLisle
  0 siblings, 1 reply; 8+ messages in thread
From: jerry DeLisle @ 2011-06-03 12:51 UTC (permalink / raw)
  To: gfortran; +Cc: gcc patches

[-- Attachment #1: Type: text/plain, Size: 1514 bytes --]

Hi,

The attached patch, which includes test cases, fixes this bug by eliminating the 
code which used floating point instructions to determine the 'r' value as 
outlined in the Fortran standard under G formatting.

Essentially, the code now examines the d and e values to determine the number of 
digits before and after the decimal point and whether or not to display the 'E' 
exponent symbol.  Adjustments are made for various corner cases, including when 
rounding has resulted in a carry. (see PR for details of the trials)

This patch is intrusive.  It results in a minor performance improvement. It 
eliminates a bit of code.

Regression tested on x86-64.

OK for trunk? then later back port to 4.6 after some proving time?

Jerry

2011-06-03  Jerry DeLisle <jvdelisle@gcc.gnu.org>

     PR libgfortran/48906
     * write.c (write_d, write_e, write_f, write_en, write_es, write_real,
     write_real_g0): Eliminate compensating flag for precision adjustment.
     * io/write_float.def (output_float): Add new comments. Add two new
     variables lblank and tblank for keeping track of required leading and
     trailing blanks. Add new code for the FMT_G case to determine digits
     before and after the decimal point. Update error messages. Reorder
     some of the existing code, separating portions handling formatting and
     portions doing the actual output.
     (CALCULATE_EXP): Delete this macro.
     (OUTPUT_FLOAT): Delete this macro.
     (write_float): Eliminate precision compensating flag.

[-- Attachment #2: pr48906-final.diff --]
[-- Type: text/x-patch, Size: 20294 bytes --]

Index: gcc/testsuite/gfortran.dg/pr20755.f
===================================================================
--- gcc/testsuite/gfortran.dg/pr20755.f	(revision 174320)
+++ gcc/testsuite/gfortran.dg/pr20755.f	(working copy)
@@ -5,8 +5,8 @@
       character*30 s
       
       write (s,2000) 0.0, 0.02
-      if (s .ne. "    0.00       2.000E-02") call abort
+      if (s .ne. "    0.00       0.200E-01") call abort
       write (s,2000) 0.01, 0.02
-      if (s .ne. "   1.000E-02   2.000E-02") call abort
+      if (s .ne. "   0.100E-01   0.200E-01") call abort
  2000 format (1PG12.3,G12.3)
       end
Index: gcc/testsuite/gfortran.dg/char4_iunit_1.f03
===================================================================
--- gcc/testsuite/gfortran.dg/char4_iunit_1.f03	(revision 174320)
+++ gcc/testsuite/gfortran.dg/char4_iunit_1.f03	(working copy)
@@ -24,7 +24,7 @@ program char4_iunit_1
   write(string, *) .true., .false. , .true.
   if (string .ne. 4_" T F T                                    ") call abort
   write(string, *) 1.2345e-06, 4.2846e+10_8
-  if (string .ne. 4_"   1.23450002E-06   42846000000.000000      ") call abort
+  if (string .ne. 4_"  1.234500019E-06   42846000000.000000      ") call abort
   write(string, *) nan, inf
   if (string .ne. 4_"              NaN         Infinity    ") call abort
   write(string, '(10x,f3.1,3x,f9.1)') nan, inf
Index: gcc/testsuite/gfortran.dg/fmt_g0_6.f08
===================================================================
--- gcc/testsuite/gfortran.dg/fmt_g0_6.f08	(revision 174320)
+++ gcc/testsuite/gfortran.dg/fmt_g0_6.f08	(working copy)
@@ -57,7 +57,7 @@ contains
             do dec = d, 0, -1
                 lower = 10.0_RT ** (d - 1 - dec) - r * 10.0_RT ** (- dec - 1)
                 upper = 10.0_RT ** (d - dec) - r * 10.0_RT ** (- dec)
-                if (lower <= mag .and. mag < upper) then
+                if (lower < mag .and. mag <= upper) then
                     write(fmt_f, "('R', a, ',F', i0, '.', i0, ',', i0, 'X')") roundmode, w - n, dec, n
                     exit
                 end if
Index: libgfortran/io/write.c
===================================================================
--- libgfortran/io/write.c	(revision 174320)
+++ libgfortran/io/write.c	(working copy)
@@ -1155,35 +1155,35 @@ write_z (st_parameter_dt *dtp, const fnode *f, con
 void
 write_d (st_parameter_dt *dtp, const fnode *f, const char *p, int len)
 {
-  write_float (dtp, f, p, len, 0);
+  write_float (dtp, f, p, len);
 }
 
 
 void
 write_e (st_parameter_dt *dtp, const fnode *f, const char *p, int len)
 {
-  write_float (dtp, f, p, len, 0);
+  write_float (dtp, f, p, len);
 }
 
 
 void
 write_f (st_parameter_dt *dtp, const fnode *f, const char *p, int len)
 {
-  write_float (dtp, f, p, len, 0);
+  write_float (dtp, f, p, len);
 }
 
 
 void
 write_en (st_parameter_dt *dtp, const fnode *f, const char *p, int len)
 {
-  write_float (dtp, f, p, len, 0);
+  write_float (dtp, f, p, len);
 }
 
 
 void
 write_es (st_parameter_dt *dtp, const fnode *f, const char *p, int len)
 {
-  write_float (dtp, f, p, len, 0);
+  write_float (dtp, f, p, len);
 }
 
 
@@ -1476,7 +1476,7 @@ write_real (st_parameter_dt *dtp, const char *sour
   int org_scale = dtp->u.p.scale_factor;
   dtp->u.p.scale_factor = 1;
   set_fnode_default (dtp, &f, length);
-  write_float (dtp, &f, source , length, 1);
+  write_float (dtp, &f, source , length);
   dtp->u.p.scale_factor = org_scale;
 }
 
@@ -1487,19 +1487,11 @@ void
 write_real_g0 (st_parameter_dt *dtp, const char *source, int length, int d)
 {
   fnode f;
-  int comp_d; 
   set_fnode_default (dtp, &f, length);
   if (d > 0)
     f.u.real.d = d;
-
-  /* Compensate for extra digits when using scale factor, d is not
-     specified, and the magnitude is such that E editing is used.  */
-  if (dtp->u.p.scale_factor > 0 && d == 0)
-    comp_d = 1;
-  else
-    comp_d = 0;
   dtp->u.p.g0_no_blanks = 1;
-  write_float (dtp, &f, source , length, comp_d);
+  write_float (dtp, &f, source , length);
   dtp->u.p.g0_no_blanks = 0;
 }
 
Index: libgfortran/io/write_float.def
===================================================================
--- libgfortran/io/write_float.def	(revision 174320)
+++ libgfortran/io/write_float.def	(working copy)
@@ -59,8 +59,22 @@ calculate_sign (st_parameter_dt *dtp, int negative
 }
 
 
-/* Output a real number according to its format which is FMT_G free.  */
+/* Output a real number according to its format.
 
+   What this does:
+
+	- Calculates the number of digits to emit before and after the
+	  decimal point.
+	- Performs the required rounding.
+	- Calculates the number of leading blanks to emit.
+	- Calculates the number of trailing blanks to emit.
+	- Calculates the exponent format if any.
+	- Determines if the sign should be displayed.
+	- Allocates the write buffer according to the width.
+	- Emits the formatted number according to the above.
+
+   buffer contains the digit string when we enter this function.  */
+
 static try
 output_float (st_parameter_dt *dtp, const fnode *f, char *buffer, size_t size, 
 	      int sign_bit, bool zero_flag, int ndigits, int edigits)
@@ -78,8 +92,10 @@ output_float (st_parameter_dt *dtp, const fnode *f
   int nafter;
   /* Number of zeros after the decimal point, whatever the precision.  */
   int nzero_real;
+  /* Flag for optional leading zero, if there is room.  */
   int leadzero;
-  int nblanks;
+  /* Number of leading and trailing blanks.  */
+  int lblanks, tblanks; 
   sign_t sign;
 
   ft = f->format;
@@ -87,7 +103,6 @@ output_float (st_parameter_dt *dtp, const fnode *f
   d = f->u.real.d;
   p = dtp->u.p.scale_factor;
 
-  rchar = '5';
   nzero_real = -1;
 
   /* We should always know the field width and precision.  */
@@ -142,19 +157,52 @@ output_float (st_parameter_dt *dtp, const fnode *f
       expchar = 0;
       break;
 
+    case FMT_G:
+      nbefore = 0;
+      nzero = 0;
+      nafter = d;
+      expchar = 0;
+
+      if (e == d)
+	{
+	  if (e > 0)
+	    {
+	      nbefore = e;
+	      nafter = d > nbefore ? d - nbefore : 0;
+	    }
+	  else
+	    expchar = 'E';
+	  break;
+	}
+      else if (e < d && e >= 0)
+	{
+	  nbefore = e;
+	  nafter = d > nbefore ? d - nbefore : 0;
+	  break;
+	}
+      else if (e >= -d && e < 0 && f->u.real.e == -1)
+	{
+	  nafter = d;
+	  if (e <= -1)
+	    expchar = 'E';
+	  break;
+	}
+
+    /* Fall through.  */
+
     case FMT_E:
     case FMT_D:
-      i = dtp->u.p.scale_factor;
       if (d <= 0 && p == 0)
 	{
 	  generate_error (&dtp->common, LIBERROR_FORMAT, "Precision not "
-			  "greater than zero in format specifier 'E' or 'D'");
+			  "greater than zero in format specifier 'D', 'E', "
+			  "or 'G'");
 	  return FAILURE;
 	}
       if (p <= -d || p >= d + 2)
 	{
 	  generate_error (&dtp->common, LIBERROR_FORMAT, "Scale factor "
-			  "out of range in format specifier 'E' or 'D'");
+			  "out of range in format specifier 'E', 'D', or 'G'");
 	  return FAILURE;
 	}
 
@@ -179,10 +227,11 @@ output_float (st_parameter_dt *dtp, const fnode *f
 	  nafter = d;
 	}
 
-      if (ft == FMT_E)
+      if (ft == FMT_E || ft == FMT_G)
 	expchar = 'E';
       else
 	expchar = 'D';
+      
       break;
 
     case FMT_EN:
@@ -224,6 +273,7 @@ output_float (st_parameter_dt *dtp, const fnode *f
 
   /* Round the value.  The value being rounded is an unsigned magnitude.
      The ROUND_COMPATIBLE is rounding away from zero when there is a tie.  */
+  rchar = '5';
   switch (dtp->u.p.current_unit->round_status)
     {
       case ROUND_ZERO: /* Do nothing and truncation occurs.  */
@@ -275,6 +325,7 @@ output_float (st_parameter_dt *dtp, const fnode *f
   rchar = '0';
   if (w > 0 && d == 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++)
     {
@@ -321,6 +372,7 @@ output_float (st_parameter_dt *dtp, const fnode *f
 	         zero.  */
 	      digits--;
 	      digits[0] = '1';
+
 	      if (ft == FMT_F)
 		{
 		  if (nzero > 0)
@@ -331,6 +383,27 @@ output_float (st_parameter_dt *dtp, const fnode *f
 		  else
 		    nbefore++;
 		}
+	      else if (ft == FMT_G)
+		{
+		  if (e < d && nbefore == 0 && e >= 0)
+		    {
+		      nbefore++;
+		      nafter--;
+		    }
+		  else if (e == d)
+		    {
+		      nbefore -= d;
+		      nafter += d;
+		      e++;
+		      expchar = 'E';
+		    }
+		  else
+		    {
+		      e++;
+		      if (e == 0)
+		        expchar = 0;
+		    }
+		}
 	      else if (ft == FMT_EN)
 		{
 		  nbefore++;
@@ -346,11 +419,69 @@ output_float (st_parameter_dt *dtp, const fnode *f
 	}
     }
 
+  /* Scan the digits string and count the number of zeros.  If we make it
+     all the way through the loop, we know the value is zero after the
+     rounding completed above. To format properly, we need to know if the
+     rounded result is zero and if so, we set the zero_flag which may have
+     been already set for actual zero.  */
+  for (i = 0; i < ndigits; i++)
+    {
+      if (digits[i] != '0')
+	break;
+    }
+  if (i == ndigits)
+    {
+      zero_flag = true;
+      /* The output is zero, so set the sign according to the sign bit unless
+	 -fno-sign-zero was specified.  */
+      if (compile_options.sign_zero == 1)
+        sign = calculate_sign (dtp, sign_bit);
+      else
+	sign = calculate_sign (dtp, 0);
+    }
+
   skip:
 
+  /* Pick a field size if none was specified, taking into account small
+     values that may have been rounded to zero.  */
+  if (w <= 0)
+    {
+      if (zero_flag)
+	w = d + (sign != S_NONE ? 2 : 1) + (d == 0 ? 1 : 0);
+      else
+	{
+	  w = nbefore + nzero + nafter + (sign != S_NONE ? 2 : 1);
+	  w = w == 1 ? 2 : w;
+	}
+    }
+  
+  /* Figure out the leading and trailing blanks.  */
+  lblanks = tblanks = 0;
+
+  /* G formatting is a special case of several things.  */
+  if (ft == FMT_G)
+    {
+      if (zero_flag && nbefore == 0 && nafter > 0)
+	{
+	  nbefore++;
+	  nafter--;
+	}
+      edigits = f->u.real.e == -1 ? 4 : f->u.real.e + 2;
+      tblanks = f->u.real.e == -1 ? 4 : f->u.real.e + 2;
+      lblanks = w - (nbefore + nzero + nafter + 1 + tblanks +
+		     (sign != S_NONE ? 1 : 0));
+    }
+  else
+    {
+      lblanks = w - (nbefore + nzero + nafter + 1 + e);
+      if (sign != S_NONE)
+	lblanks--;
+    }
+
   /* Calculate the format of the exponent field.  */
   if (expchar)
     {
+      tblanks = 0;
       edigits = 1;
       for (i = abs (e); i >= 10; i /= 10)
 	edigits++;
@@ -379,76 +510,57 @@ output_float (st_parameter_dt *dtp, const fnode *f
   else
     edigits = 0;
 
-  /* Scan the digits string and count the number of zeros.  If we make it
-     all the way through the loop, we know the value is zero after the
-     rounding completed above.  */
-  for (i = 0; i < ndigits; i++)
+  if (ft == FMT_G)
     {
-      if (digits[i] != '0')
-	break;
-    }
-
-  /* To format properly, we need to know if the rounded result is zero and if
-     so, we set the zero_flag which may have been already set for
-     actual zero.  */
-  if (i == ndigits)
-    {
-      zero_flag = true;
-      /* The output is zero, so set the sign according to the sign bit unless
-	 -fno-sign-zero was specified.  */
-      if (compile_options.sign_zero == 1)
-        sign = calculate_sign (dtp, sign_bit);
-      else
-	sign = calculate_sign (dtp, 0);
-    }
-
-  /* Pick a field size if none was specified, taking into account small
-     values that may have been rounded to zero.  */
-  if (w <= 0)
-    {
-      if (zero_flag)
-	w = d + (sign != S_NONE ? 2 : 1) + (d == 0 ? 1 : 0);
-      else
+      if (dtp->u.p.g0_no_blanks)
 	{
-	  w = nbefore + nzero + nafter + (sign != S_NONE ? 2 : 1);
-	  w = w == 1 ? 2 : w;
+	  w = w - tblanks - lblanks;
+	  lblanks = tblanks = 0;
 	}
+      else if (expchar == 0)
+	{
+	  lblanks = w - (nbefore + nzero + nafter + 1 + tblanks +
+			 (sign != S_NONE ? 1 : 0));
+	  if (dtp->u.p.no_leading_blank)
+	    {
+	      tblanks += lblanks;
+	      lblanks = 0; 
+	    }
+	}
     }
-  
-  /* Work out how much padding is needed.  */
-  nblanks = w - (nbefore + nzero + nafter + edigits + 1);
-  if (sign != S_NONE)
-    nblanks--;
+  else
+    lblanks = w - (nbefore + nzero + nafter + 1 + edigits +
+		 (sign != S_NONE ? 1 : 0));
 
-  if (dtp->u.p.g0_no_blanks)
-    {
-      w -= nblanks;
-      nblanks = 0;
-    }
-
   /* Create the ouput buffer.  */
-  out = write_block (dtp, w);
+  if ((ft == FMT_G && dtp->u.p.g0_no_blanks)
+      || (ft == FMT_F && f->u.real.w == 0))
+    out = write_block (dtp, w);
+  else
+    out = write_block (dtp, f->u.real.w);
+
   if (out == NULL)
     return FAILURE;
 
   /* Check the value fits in the specified field width.  */
-  if (nblanks < 0 || edigits == -1 || w == 1 || (w == 2 && sign != S_NONE))
+  if (lblanks < 0 || tblanks < 0 || edigits == -1 || w == 1
+      || (w == 2 && sign != S_NONE))
     {
       if (unlikely (is_char4_unit (dtp)))
 	{
 	  gfc_char4_t *out4 = (gfc_char4_t *) out;
-	  memset4 (out4, '*', w);
+	  memset4 (out4, '*', f->u.real.w);
 	  return FAILURE;
 	}
-      star_fill (out, w);
+      star_fill (out, f->u.real.w);
       return FAILURE;
     }
 
   /* See if we have space for a zero before the decimal point.  */
-  if (nbefore == 0 && nblanks > 0)
+  if (nbefore == 0 && lblanks > 0)
     {
       leadzero = 1;
-      nblanks--;
+      lblanks--;
     }
   else
     leadzero = 0;
@@ -461,10 +573,10 @@ output_float (st_parameter_dt *dtp, const fnode *f
       gfc_char4_t *out4 = (gfc_char4_t *) out;
       /* Pad to full field width.  */
 
-      if ( ( nblanks > 0 ) && !dtp->u.p.no_leading_blank)
+      if ((lblanks > 0) && !dtp->u.p.no_leading_blank)
 	{
-	  memset4 (out4, ' ', nblanks);
-	  out4 += nblanks;
+	  memset4 (out4, ' ', lblanks);
+	  out4 += lblanks;
 	}
 
       /* Output the initial sign (if any).  */
@@ -539,21 +651,17 @@ output_float (st_parameter_dt *dtp, const fnode *f
 	  memcpy4 (out4, buffer, edigits);
 	}
 
-      if (dtp->u.p.no_leading_blank)
-	{
-	  out4 += edigits;
-	  memset4 (out4, ' ' , nblanks);
-	  dtp->u.p.no_leading_blank = 0;
-	}
+      /* Emit trailing blanks if needed.  */
+      memset4 (out4, ' ' , tblanks);
       return SUCCESS;
     } /* End of character(kind=4) internal unit code.  */
 
-  /* Pad to full field width.  */
+  /* Pad leading to full field width.  */
 
-  if ( ( nblanks > 0 ) && !dtp->u.p.no_leading_blank)
+  if (lblanks > 0 && !dtp->u.p.no_leading_blank)
     {
-      memset (out, ' ', nblanks);
-      out += nblanks;
+      memset (out, ' ', lblanks);
+      out += lblanks;
     }
 
   /* Output the initial sign (if any).  */
@@ -627,13 +735,10 @@ output_float (st_parameter_dt *dtp, const fnode *f
       memcpy (out, buffer, edigits);
     }
 
-  if (dtp->u.p.no_leading_blank)
-    {
-      out += edigits;
-      memset( out , ' ' , nblanks );
-      dtp->u.p.no_leading_blank = 0;
-    }
-
+  /* Emit trailing blanks if needed.  */
+  if (tblanks)
+    memset( out , ' ' , tblanks );
+  
 #undef STR
 #undef STR1
 #undef MIN_FIELD_WIDTH
@@ -770,185 +875,6 @@ write_infnan (st_parameter_dt *dtp, const fnode *f
     }
 }
 
-
-/* Returns the value of 10**d.  */
-
-#define CALCULATE_EXP(x) \
-inline static GFC_REAL_ ## x \
-calculate_exp_ ## x  (int d)\
-{\
-  int i;\
-  GFC_REAL_ ## x r = 1.0;\
-  for (i = 0; i< (d >= 0 ? d : -d); i++)\
-    r *= 10;\
-  r = (d >= 0) ? r : 1.0 / r;\
-  return r;\
-}
-
-CALCULATE_EXP(4)
-
-CALCULATE_EXP(8)
-
-#ifdef HAVE_GFC_REAL_10
-CALCULATE_EXP(10)
-#endif
-
-#ifdef HAVE_GFC_REAL_16
-CALCULATE_EXP(16)
-#endif
-#undef CALCULATE_EXP
-
-/* Generate corresponding I/O format for FMT_G and output.
-   The rules to translate FMT_G to FMT_E or FMT_F from DEC fortran
-   LRM (table 11-2, Chapter 11, "I/O Formatting", P11-25) is:
-
-   Data Magnitude                              Equivalent Conversion
-   0< m < 0.1-0.5*10**(-d-1)                   Ew.d[Ee]
-   m = 0                                       F(w-n).(d-1), n' '
-   0.1-0.5*10**(-d-1)<= m < 1-0.5*10**(-d)     F(w-n).d, n' '
-   1-0.5*10**(-d)<= m < 10-0.5*10**(-d+1)      F(w-n).(d-1), n' '
-   10-0.5*10**(-d+1)<= m < 100-0.5*10**(-d+2)  F(w-n).(d-2), n' '
-   ................                           ..........
-   10**(d-1)-0.5*10**(-1)<= m <10**d-0.5       F(w-n).0,n(' ')
-   m >= 10**d-0.5                              Ew.d[Ee]
-
-   notes: for Gw.d ,  n' ' means 4 blanks
-	  for Gw.dEe, n' ' means e+2 blanks
-	  for rounding modes adjustment, r, See Fortran F2008 10.7.5.2.2
-	  the asm volatile is required for 32-bit x86 platforms.  */
-
-#define OUTPUT_FLOAT_FMT_G(x) \
-static void \
-output_float_FMT_G_ ## x (st_parameter_dt *dtp, const fnode *f, \
-		      GFC_REAL_ ## x m, char *buffer, size_t size, \
-		      int sign_bit, bool zero_flag, int ndigits, \
-                      int edigits, int comp_d) \
-{ \
-  int e = f->u.real.e;\
-  int d = f->u.real.d;\
-  int w = f->u.real.w;\
-  fnode *newf;\
-  GFC_REAL_ ## x rexp_d, r = 0.5;\
-  int low, high, mid;\
-  int ubound, lbound;\
-  char *p, pad = ' ';\
-  int save_scale_factor, nb = 0;\
-  try result;\
-\
-  save_scale_factor = dtp->u.p.scale_factor;\
-  newf = (fnode *) get_mem (sizeof (fnode));\
-\
-  switch (dtp->u.p.current_unit->round_status)\
-    {\
-      case ROUND_ZERO:\
-	r = sign_bit ? 1.0 : 0.0;\
-	break;\
-      case ROUND_UP:\
-	r = 1.0;\
-	break;\
-      case ROUND_DOWN:\
-	r = 0.0;\
-	break;\
-      default:\
-	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)))\
-      || ((m == 0.0) && !(compile_options.allow_std\
-			  & (GFC_STD_F2003 | GFC_STD_F2008))))\
-    { \
-      newf->format = FMT_E;\
-      newf->u.real.w = w;\
-      newf->u.real.d = d - comp_d;\
-      newf->u.real.e = e;\
-      nb = 0;\
-      goto finish;\
-    }\
-\
-  mid = 0;\
-  low = 0;\
-  high = d + 1;\
-  lbound = 0;\
-  ubound = d + 1;\
-\
-  while (low <= high)\
-    { \
-      volatile GFC_REAL_ ## x temp;\
-      mid = (low + high) / 2;\
-\
-      temp = (calculate_exp_ ## x (mid - 1) * (1 - r * rexp_d));\
-\
-      if (m < temp)\
-        { \
-          ubound = mid;\
-          if (ubound == lbound + 1)\
-            break;\
-          high = mid - 1;\
-        }\
-      else if (m > temp)\
-        { \
-          lbound = mid;\
-          if (ubound == lbound + 1)\
-            { \
-              mid ++;\
-              break;\
-            }\
-          low = mid + 1;\
-        }\
-      else\
-	{\
-	  mid++;\
-	  break;\
-	}\
-    }\
-\
-  nb = e <= 0 ? 4 : e + 2;\
-  nb = nb >= w ? w - 1 : nb;\
-  newf->format = FMT_F;\
-  newf->u.real.w = w - nb;\
-  newf->u.real.d = m == 0.0 ? d - 1 : -(mid - d - 1) ;\
-  dtp->u.p.scale_factor = 0;\
-\
- finish:\
-  result = output_float (dtp, newf, buffer, size, sign_bit, zero_flag, \
-			 ndigits, edigits);\
-  dtp->u.p.scale_factor = save_scale_factor;\
-\
-  free (newf);\
-\
-  if (nb > 0 && !dtp->u.p.g0_no_blanks)\
-    {\
-      p = write_block (dtp, nb);\
-      if (p == NULL)\
-	return;\
-      if (result == FAILURE)\
-        pad = '*';\
-      if (unlikely (is_char4_unit (dtp)))\
-	{\
-	  gfc_char4_t *p4 = (gfc_char4_t *) p;\
-	  memset4 (p4, pad, nb);\
-	}\
-      else \
-	memset (p, pad, nb);\
-    }\
-}\
-
-OUTPUT_FLOAT_FMT_G(4)
-
-OUTPUT_FLOAT_FMT_G(8)
-
-#ifdef HAVE_GFC_REAL_10
-OUTPUT_FLOAT_FMT_G(10)
-#endif
-
-#ifdef HAVE_GFC_REAL_16
-OUTPUT_FLOAT_FMT_G(16)
-#endif
-
-#undef OUTPUT_FLOAT_FMT_G
-
-
 /* Define a macro to build code for write_float.  */
 
   /* Note: Before output_float is called, snprintf is used to print to buffer the
@@ -1003,19 +929,15 @@ __qmath_(quadmath_snprintf) (buffer, sizeof buffer
 \
 	DTOA ## y\
 \
-	if (f->format != FMT_G)\
-	  output_float (dtp, f, buffer, size, sign_bit, zero_flag, ndigits, \
-			edigits);\
-	else \
-	  output_float_FMT_G_ ## x (dtp, f, tmp, buffer, size, sign_bit, \
-				    zero_flag, ndigits, edigits, comp_d);\
+	output_float (dtp, f, buffer, size, sign_bit, zero_flag, ndigits, \
+		      edigits);\
 }\
 
 /* Output a real number according to its format.  */
 
 static void
 write_float (st_parameter_dt *dtp, const fnode *f, const char *source, \
-            int len, int comp_d)
+            int len)
 {
 
 #if defined(HAVE_GFC_REAL_16) || __LDBL_DIG__ > 18

^ permalink raw reply	[flat|nested] 8+ messages in thread

* Re: [patch, libgfortran] PR48906 Wrong rounding results with -m32
  2011-06-03 12:51 [patch, libgfortran] Wrong rounding results with -m32 jerry DeLisle
@ 2011-06-10 17:39 ` jerry DeLisle
  2011-06-11  8:24   ` Janne Blomqvist
  0 siblings, 1 reply; 8+ messages in thread
From: jerry DeLisle @ 2011-06-10 17:39 UTC (permalink / raw)
  To: jerry DeLisle; +Cc: gfortran, gcc patches, Thomas Henlich

[-- Attachment #1: Type: text/plain, Size: 1077 bytes --]

On 06/03/2011 05:51 AM, jerry DeLisle wrote:
> Hi,
>
> The attached patch, which includes test cases, fixes this bug by eliminating the
> code which used floating point instructions to determine the 'r' value as
> outlined in the Fortran standard under G formatting.
>
> Essentially, the code now examines the d and e values to determine the number of
> digits before and after the decimal point and whether or not to display the 'E'
> exponent symbol. Adjustments are made for various corner cases, including when
> rounding has resulted in a carry. (see PR for details of the trials)
>
> This patch is intrusive. It results in a minor performance improvement. It
> eliminates a bit of code.
>
> Regression tested on x86-64.
>
> OK for trunk? then later back port to 4.6 after some proving time?

Attached is updated patch based on comments from Thomas Henlich.  See my comment 
#33 and subsequent comments in the PR48906 explaining the issue with test case 
fmt_g0_6.f08 which is revised by the updated patch.

Regression tested on X86-64.

OK for trunk. (please ;) )

Jerry

[-- Attachment #2: pr48906-final-2.diff --]
[-- Type: text/x-patch, Size: 19578 bytes --]

Index: gcc/testsuite/gfortran.dg/char4_iunit_1.f03
===================================================================
--- gcc/testsuite/gfortran.dg/char4_iunit_1.f03	(revision 174673)
+++ gcc/testsuite/gfortran.dg/char4_iunit_1.f03	(working copy)
@@ -24,7 +24,7 @@ program char4_iunit_1
   write(string, *) .true., .false. , .true.
   if (string .ne. 4_" T F T                                    ") call abort
   write(string, *) 1.2345e-06, 4.2846e+10_8
-  if (string .ne. 4_"   1.23450002E-06   42846000000.000000      ") call abort
+  if (string .ne. 4_"  1.234500019E-06   42846000000.000000      ") call abort
   write(string, *) nan, inf
   if (string .ne. 4_"              NaN         Infinity    ") call abort
   write(string, '(10x,f3.1,3x,f9.1)') nan, inf
Index: gcc/testsuite/gfortran.dg/fmt_g0_6.f08
===================================================================
--- gcc/testsuite/gfortran.dg/fmt_g0_6.f08	(revision 174673)
+++ gcc/testsuite/gfortran.dg/fmt_g0_6.f08	(working copy)
@@ -9,7 +9,7 @@ program test_g0fr
     
     call check_all(0.0_RT, 15, 2, 0)
     call check_all(0.991_RT, 15, 2, 0)
-    call check_all(0.995_RT, 15, 2, 0)
+    call check_all(0.99500000001_RT, 15, 2, 0)
     call check_all(0.996_RT, 15, 2, 0)
     call check_all(0.999_RT, 15, 2, 0)
 contains
Index: libgfortran/io/write.c
===================================================================
--- libgfortran/io/write.c	(revision 174673)
+++ libgfortran/io/write.c	(working copy)
@@ -1155,35 +1155,35 @@ write_z (st_parameter_dt *dtp, const fnode *f, con
 void
 write_d (st_parameter_dt *dtp, const fnode *f, const char *p, int len)
 {
-  write_float (dtp, f, p, len, 0);
+  write_float (dtp, f, p, len);
 }
 
 
 void
 write_e (st_parameter_dt *dtp, const fnode *f, const char *p, int len)
 {
-  write_float (dtp, f, p, len, 0);
+  write_float (dtp, f, p, len);
 }
 
 
 void
 write_f (st_parameter_dt *dtp, const fnode *f, const char *p, int len)
 {
-  write_float (dtp, f, p, len, 0);
+  write_float (dtp, f, p, len);
 }
 
 
 void
 write_en (st_parameter_dt *dtp, const fnode *f, const char *p, int len)
 {
-  write_float (dtp, f, p, len, 0);
+  write_float (dtp, f, p, len);
 }
 
 
 void
 write_es (st_parameter_dt *dtp, const fnode *f, const char *p, int len)
 {
-  write_float (dtp, f, p, len, 0);
+  write_float (dtp, f, p, len);
 }
 
 
@@ -1476,7 +1476,7 @@ write_real (st_parameter_dt *dtp, const char *sour
   int org_scale = dtp->u.p.scale_factor;
   dtp->u.p.scale_factor = 1;
   set_fnode_default (dtp, &f, length);
-  write_float (dtp, &f, source , length, 1);
+  write_float (dtp, &f, source , length);
   dtp->u.p.scale_factor = org_scale;
 }
 
@@ -1487,19 +1487,11 @@ void
 write_real_g0 (st_parameter_dt *dtp, const char *source, int length, int d)
 {
   fnode f;
-  int comp_d; 
   set_fnode_default (dtp, &f, length);
   if (d > 0)
     f.u.real.d = d;
-
-  /* Compensate for extra digits when using scale factor, d is not
-     specified, and the magnitude is such that E editing is used.  */
-  if (dtp->u.p.scale_factor > 0 && d == 0)
-    comp_d = 1;
-  else
-    comp_d = 0;
   dtp->u.p.g0_no_blanks = 1;
-  write_float (dtp, &f, source , length, comp_d);
+  write_float (dtp, &f, source , length);
   dtp->u.p.g0_no_blanks = 0;
 }
 
Index: libgfortran/io/write_float.def
===================================================================
--- libgfortran/io/write_float.def	(revision 174673)
+++ libgfortran/io/write_float.def	(working copy)
@@ -59,8 +59,22 @@ calculate_sign (st_parameter_dt *dtp, int negative
 }
 
 
-/* Output a real number according to its format which is FMT_G free.  */
+/* Output a real number according to its format.
 
+   What this does:
+
+	- Calculates the number of digits to emit before and after the
+	  decimal point.
+	- Performs the required rounding.
+	- Calculates the number of leading blanks to emit.
+	- Calculates the number of trailing blanks to emit.
+	- Calculates the exponent format if any.
+	- Determines if the sign should be displayed.
+	- Allocates the write buffer according to the width.
+	- Emits the formatted number according to the above.
+
+   buffer contains the digit string when we enter this function.  */
+
 static try
 output_float (st_parameter_dt *dtp, const fnode *f, char *buffer, size_t size, 
 	      int sign_bit, bool zero_flag, int ndigits, int edigits)
@@ -78,8 +92,10 @@ output_float (st_parameter_dt *dtp, const fnode *f
   int nafter;
   /* Number of zeros after the decimal point, whatever the precision.  */
   int nzero_real;
+  /* Flag for optional leading zero, if there is room.  */
   int leadzero;
-  int nblanks;
+  /* Number of leading and trailing blanks.  */
+  int lblanks, tblanks; 
   sign_t sign;
 
   ft = f->format;
@@ -87,7 +103,6 @@ output_float (st_parameter_dt *dtp, const fnode *f
   d = f->u.real.d;
   p = dtp->u.p.scale_factor;
 
-  rchar = '5';
   nzero_real = -1;
 
   /* We should always know the field width and precision.  */
@@ -142,19 +157,50 @@ output_float (st_parameter_dt *dtp, const fnode *f
       expchar = 0;
       break;
 
+    case FMT_G:
+      nbefore = 0;
+      nzero = 0;
+      nafter = d;
+      expchar = 0;
+
+      if (e == d)
+	{
+	  if (e > 0)
+	    {
+	      nbefore = e;
+	      nafter = d > nbefore ? d - nbefore : 0;
+	    }
+	  else
+	    expchar = 'E';
+	  break;
+	}
+      else if (e < d && e >= 0)
+	{
+	  nbefore = e;
+	  nafter = d > nbefore ? d - nbefore : 0;
+	  break;
+	}
+      else if (e >= -d && e < 0 && f->u.real.e == -1)
+	{
+	  nafter = d;
+	  expchar = 'E';
+	}
+
+    /* Fall through.  */
+
     case FMT_E:
     case FMT_D:
-      i = dtp->u.p.scale_factor;
       if (d <= 0 && p == 0)
 	{
 	  generate_error (&dtp->common, LIBERROR_FORMAT, "Precision not "
-			  "greater than zero in format specifier 'E' or 'D'");
+			  "greater than zero in format specifier 'D', 'E', "
+			  "or 'G'");
 	  return FAILURE;
 	}
       if (p <= -d || p >= d + 2)
 	{
 	  generate_error (&dtp->common, LIBERROR_FORMAT, "Scale factor "
-			  "out of range in format specifier 'E' or 'D'");
+			  "out of range in format specifier 'E', 'D', or 'G'");
 	  return FAILURE;
 	}
 
@@ -179,10 +225,11 @@ output_float (st_parameter_dt *dtp, const fnode *f
 	  nafter = d;
 	}
 
-      if (ft == FMT_E)
+      if (ft == FMT_E || ft == FMT_G)
 	expchar = 'E';
       else
 	expchar = 'D';
+      
       break;
 
     case FMT_EN:
@@ -224,6 +271,7 @@ output_float (st_parameter_dt *dtp, const fnode *f
 
   /* Round the value.  The value being rounded is an unsigned magnitude.
      The ROUND_COMPATIBLE is rounding away from zero when there is a tie.  */
+  rchar = '5';
   switch (dtp->u.p.current_unit->round_status)
     {
       case ROUND_ZERO: /* Do nothing and truncation occurs.  */
@@ -275,6 +323,7 @@ output_float (st_parameter_dt *dtp, const fnode *f
   rchar = '0';
   if (w > 0 && d == 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++)
     {
@@ -284,7 +333,7 @@ output_float (st_parameter_dt *dtp, const fnode *f
   goto skip;
     
   do_rnd:
- 
+
   if (nbefore + nafter == 0)
     {
       ndigits = 0;
@@ -321,6 +370,7 @@ output_float (st_parameter_dt *dtp, const fnode *f
 	         zero.  */
 	      digits--;
 	      digits[0] = '1';
+
 	      if (ft == FMT_F)
 		{
 		  if (nzero > 0)
@@ -331,6 +381,27 @@ output_float (st_parameter_dt *dtp, const fnode *f
 		  else
 		    nbefore++;
 		}
+	      else if (ft == FMT_G)
+		{
+		  if (e < d && e >= 0)
+		    {
+		      nbefore++;
+		      nafter--;
+		    }
+		  else if (e == d)
+		    {
+		      nbefore -= d;
+		      nafter += d;
+		      e++;
+		      expchar = 'E';
+		    }
+		  else
+		    {
+		      e++;
+		      if (e == 0)
+		        expchar = 0;
+		    }
+		}
 	      else if (ft == FMT_EN)
 		{
 		  nbefore++;
@@ -346,11 +417,68 @@ output_float (st_parameter_dt *dtp, const fnode *f
 	}
     }
 
+  /* Scan the digits string and count the number of zeros.  If we make it
+     all the way through the loop, we know the value is zero after the
+     rounding completed above. To format properly, we need to know if the
+     rounded result is zero and if so, we set the zero_flag which may have
+     been already set for actual zero.  */
+  for (i = 0; i < ndigits; i++)
+    {
+      if (digits[i] != '0')
+	break;
+    }
+  if (i == ndigits)
+    {
+      zero_flag = true;
+      /* The output is zero, so set the sign according to the sign bit unless
+	 -fno-sign-zero was specified.  */
+      if (compile_options.sign_zero == 1)
+        sign = calculate_sign (dtp, sign_bit);
+      else
+	sign = calculate_sign (dtp, 0);
+    }
+
   skip:
 
+  /* Pick a field size if none was specified, taking into account small
+     values that may have been rounded to zero.  */
+  if (w <= 0)
+    {
+      if (zero_flag)
+	w = d + (sign != S_NONE ? 2 : 1) + (d == 0 ? 1 : 0);
+      else
+	{
+	  w = nbefore + nzero + nafter + (sign != S_NONE ? 2 : 1);
+	  w = w == 1 ? 2 : w;
+	}
+    }
+  
+  /* Figure out the leading and trailing blanks.  */
+  lblanks = tblanks = 0;
+
+  /* G formatting is a special case of several things.  */
+  if (ft == FMT_G)
+    {
+      if (zero_flag && nbefore == 0 && nafter > 0)
+	{
+	  nbefore++;
+	  nafter--;
+	}
+      tblanks = edigits = f->u.real.e == -1 ? 4 : f->u.real.e + 2;
+      lblanks = w - (nbefore + nzero + nafter + 1 + tblanks +
+		     (sign != S_NONE ? 1 : 0));
+    }
+  else
+    {
+      lblanks = w - (nbefore + nzero + nafter + 1 + e);
+      if (sign != S_NONE)
+	lblanks--;
+    }
+
   /* Calculate the format of the exponent field.  */
   if (expchar)
     {
+      tblanks = 0;
       edigits = 1;
       for (i = abs (e); i >= 10; i /= 10)
 	edigits++;
@@ -379,76 +507,57 @@ output_float (st_parameter_dt *dtp, const fnode *f
   else
     edigits = 0;
 
-  /* Scan the digits string and count the number of zeros.  If we make it
-     all the way through the loop, we know the value is zero after the
-     rounding completed above.  */
-  for (i = 0; i < ndigits; i++)
+  if (ft == FMT_G)
     {
-      if (digits[i] != '0')
-	break;
-    }
-
-  /* To format properly, we need to know if the rounded result is zero and if
-     so, we set the zero_flag which may have been already set for
-     actual zero.  */
-  if (i == ndigits)
-    {
-      zero_flag = true;
-      /* The output is zero, so set the sign according to the sign bit unless
-	 -fno-sign-zero was specified.  */
-      if (compile_options.sign_zero == 1)
-        sign = calculate_sign (dtp, sign_bit);
-      else
-	sign = calculate_sign (dtp, 0);
-    }
-
-  /* Pick a field size if none was specified, taking into account small
-     values that may have been rounded to zero.  */
-  if (w <= 0)
-    {
-      if (zero_flag)
-	w = d + (sign != S_NONE ? 2 : 1) + (d == 0 ? 1 : 0);
-      else
+      if (dtp->u.p.g0_no_blanks)
 	{
-	  w = nbefore + nzero + nafter + (sign != S_NONE ? 2 : 1);
-	  w = w == 1 ? 2 : w;
+	  w = w - tblanks - lblanks;
+	  lblanks = tblanks = 0;
 	}
+      else if (expchar == 0)
+	{
+	  lblanks = w - (nbefore + nzero + nafter + 1 + tblanks +
+			 (sign != S_NONE ? 1 : 0));
+	  if (dtp->u.p.no_leading_blank)
+	    {
+	      tblanks += lblanks;
+	      lblanks = 0; 
+	    }
+	}
     }
-  
-  /* Work out how much padding is needed.  */
-  nblanks = w - (nbefore + nzero + nafter + edigits + 1);
-  if (sign != S_NONE)
-    nblanks--;
+  else
+    lblanks = w - (nbefore + nzero + nafter + 1 + edigits +
+		 (sign != S_NONE ? 1 : 0));
 
-  if (dtp->u.p.g0_no_blanks)
-    {
-      w -= nblanks;
-      nblanks = 0;
-    }
-
   /* Create the ouput buffer.  */
-  out = write_block (dtp, w);
+  if ((ft == FMT_G && dtp->u.p.g0_no_blanks)
+      || (ft == FMT_F && f->u.real.w == 0))
+    out = write_block (dtp, w);
+  else
+    out = write_block (dtp, f->u.real.w);
+
   if (out == NULL)
     return FAILURE;
 
   /* Check the value fits in the specified field width.  */
-  if (nblanks < 0 || edigits == -1 || w == 1 || (w == 2 && sign != S_NONE))
+  if (lblanks < 0 || tblanks < 0 || edigits == -1 || w == 1
+      || (w == 2 && sign != S_NONE))
     {
       if (unlikely (is_char4_unit (dtp)))
 	{
 	  gfc_char4_t *out4 = (gfc_char4_t *) out;
-	  memset4 (out4, '*', w);
+	  memset4 (out4, '*', f->u.real.w);
 	  return FAILURE;
 	}
-      star_fill (out, w);
+      star_fill (out, f->u.real.w);
       return FAILURE;
     }
 
   /* See if we have space for a zero before the decimal point.  */
-  if (nbefore == 0 && nblanks > 0)
+  if (nbefore == 0 && lblanks > 0)
     {
       leadzero = 1;
-      nblanks--;
+      lblanks--;
     }
   else
     leadzero = 0;
@@ -461,10 +570,10 @@ output_float (st_parameter_dt *dtp, const fnode *f
       gfc_char4_t *out4 = (gfc_char4_t *) out;
       /* Pad to full field width.  */
 
-      if ( ( nblanks > 0 ) && !dtp->u.p.no_leading_blank)
+      if ((lblanks > 0) && !dtp->u.p.no_leading_blank)
 	{
-	  memset4 (out4, ' ', nblanks);
-	  out4 += nblanks;
+	  memset4 (out4, ' ', lblanks);
+	  out4 += lblanks;
 	}
 
       /* Output the initial sign (if any).  */
@@ -539,21 +648,17 @@ output_float (st_parameter_dt *dtp, const fnode *f
 	  memcpy4 (out4, buffer, edigits);
 	}
 
-      if (dtp->u.p.no_leading_blank)
-	{
-	  out4 += edigits;
-	  memset4 (out4, ' ' , nblanks);
-	  dtp->u.p.no_leading_blank = 0;
-	}
+      /* Emit trailing blanks if needed.  */
+      memset4 (out4, ' ' , tblanks);
       return SUCCESS;
     } /* End of character(kind=4) internal unit code.  */
 
-  /* Pad to full field width.  */
+  /* Pad leading to full field width.  */
 
-  if ( ( nblanks > 0 ) && !dtp->u.p.no_leading_blank)
+  if (lblanks > 0 && !dtp->u.p.no_leading_blank)
     {
-      memset (out, ' ', nblanks);
-      out += nblanks;
+      memset (out, ' ', lblanks);
+      out += lblanks;
     }
 
   /* Output the initial sign (if any).  */
@@ -627,13 +732,10 @@ output_float (st_parameter_dt *dtp, const fnode *f
       memcpy (out, buffer, edigits);
     }
 
-  if (dtp->u.p.no_leading_blank)
-    {
-      out += edigits;
-      memset( out , ' ' , nblanks );
-      dtp->u.p.no_leading_blank = 0;
-    }
-
+  /* Emit trailing blanks if needed.  */
+  if (tblanks)
+    memset( out , ' ' , tblanks );
+  
 #undef STR
 #undef STR1
 #undef MIN_FIELD_WIDTH
@@ -770,185 +872,6 @@ write_infnan (st_parameter_dt *dtp, const fnode *f
     }
 }
 
-
-/* Returns the value of 10**d.  */
-
-#define CALCULATE_EXP(x) \
-inline static GFC_REAL_ ## x \
-calculate_exp_ ## x  (int d)\
-{\
-  int i;\
-  GFC_REAL_ ## x r = 1.0;\
-  for (i = 0; i< (d >= 0 ? d : -d); i++)\
-    r *= 10;\
-  r = (d >= 0) ? r : 1.0 / r;\
-  return r;\
-}
-
-CALCULATE_EXP(4)
-
-CALCULATE_EXP(8)
-
-#ifdef HAVE_GFC_REAL_10
-CALCULATE_EXP(10)
-#endif
-
-#ifdef HAVE_GFC_REAL_16
-CALCULATE_EXP(16)
-#endif
-#undef CALCULATE_EXP
-
-/* Generate corresponding I/O format for FMT_G and output.
-   The rules to translate FMT_G to FMT_E or FMT_F from DEC fortran
-   LRM (table 11-2, Chapter 11, "I/O Formatting", P11-25) is:
-
-   Data Magnitude                              Equivalent Conversion
-   0< m < 0.1-0.5*10**(-d-1)                   Ew.d[Ee]
-   m = 0                                       F(w-n).(d-1), n' '
-   0.1-0.5*10**(-d-1)<= m < 1-0.5*10**(-d)     F(w-n).d, n' '
-   1-0.5*10**(-d)<= m < 10-0.5*10**(-d+1)      F(w-n).(d-1), n' '
-   10-0.5*10**(-d+1)<= m < 100-0.5*10**(-d+2)  F(w-n).(d-2), n' '
-   ................                           ..........
-   10**(d-1)-0.5*10**(-1)<= m <10**d-0.5       F(w-n).0,n(' ')
-   m >= 10**d-0.5                              Ew.d[Ee]
-
-   notes: for Gw.d ,  n' ' means 4 blanks
-	  for Gw.dEe, n' ' means e+2 blanks
-	  for rounding modes adjustment, r, See Fortran F2008 10.7.5.2.2
-	  the asm volatile is required for 32-bit x86 platforms.  */
-
-#define OUTPUT_FLOAT_FMT_G(x) \
-static void \
-output_float_FMT_G_ ## x (st_parameter_dt *dtp, const fnode *f, \
-		      GFC_REAL_ ## x m, char *buffer, size_t size, \
-		      int sign_bit, bool zero_flag, int ndigits, \
-                      int edigits, int comp_d) \
-{ \
-  int e = f->u.real.e;\
-  int d = f->u.real.d;\
-  int w = f->u.real.w;\
-  fnode *newf;\
-  GFC_REAL_ ## x rexp_d, r = 0.5;\
-  int low, high, mid;\
-  int ubound, lbound;\
-  char *p, pad = ' ';\
-  int save_scale_factor, nb = 0;\
-  try result;\
-\
-  save_scale_factor = dtp->u.p.scale_factor;\
-  newf = (fnode *) get_mem (sizeof (fnode));\
-\
-  switch (dtp->u.p.current_unit->round_status)\
-    {\
-      case ROUND_ZERO:\
-	r = sign_bit ? 1.0 : 0.0;\
-	break;\
-      case ROUND_UP:\
-	r = 1.0;\
-	break;\
-      case ROUND_DOWN:\
-	r = 0.0;\
-	break;\
-      default:\
-	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)))\
-      || ((m == 0.0) && !(compile_options.allow_std\
-			  & (GFC_STD_F2003 | GFC_STD_F2008))))\
-    { \
-      newf->format = FMT_E;\
-      newf->u.real.w = w;\
-      newf->u.real.d = d - comp_d;\
-      newf->u.real.e = e;\
-      nb = 0;\
-      goto finish;\
-    }\
-\
-  mid = 0;\
-  low = 0;\
-  high = d + 1;\
-  lbound = 0;\
-  ubound = d + 1;\
-\
-  while (low <= high)\
-    { \
-      volatile GFC_REAL_ ## x temp;\
-      mid = (low + high) / 2;\
-\
-      temp = (calculate_exp_ ## x (mid - 1) * (1 - r * rexp_d));\
-\
-      if (m < temp)\
-        { \
-          ubound = mid;\
-          if (ubound == lbound + 1)\
-            break;\
-          high = mid - 1;\
-        }\
-      else if (m > temp)\
-        { \
-          lbound = mid;\
-          if (ubound == lbound + 1)\
-            { \
-              mid ++;\
-              break;\
-            }\
-          low = mid + 1;\
-        }\
-      else\
-	{\
-	  mid++;\
-	  break;\
-	}\
-    }\
-\
-  nb = e <= 0 ? 4 : e + 2;\
-  nb = nb >= w ? w - 1 : nb;\
-  newf->format = FMT_F;\
-  newf->u.real.w = w - nb;\
-  newf->u.real.d = m == 0.0 ? d - 1 : -(mid - d - 1) ;\
-  dtp->u.p.scale_factor = 0;\
-\
- finish:\
-  result = output_float (dtp, newf, buffer, size, sign_bit, zero_flag, \
-			 ndigits, edigits);\
-  dtp->u.p.scale_factor = save_scale_factor;\
-\
-  free (newf);\
-\
-  if (nb > 0 && !dtp->u.p.g0_no_blanks)\
-    {\
-      p = write_block (dtp, nb);\
-      if (p == NULL)\
-	return;\
-      if (result == FAILURE)\
-        pad = '*';\
-      if (unlikely (is_char4_unit (dtp)))\
-	{\
-	  gfc_char4_t *p4 = (gfc_char4_t *) p;\
-	  memset4 (p4, pad, nb);\
-	}\
-      else \
-	memset (p, pad, nb);\
-    }\
-}\
-
-OUTPUT_FLOAT_FMT_G(4)
-
-OUTPUT_FLOAT_FMT_G(8)
-
-#ifdef HAVE_GFC_REAL_10
-OUTPUT_FLOAT_FMT_G(10)
-#endif
-
-#ifdef HAVE_GFC_REAL_16
-OUTPUT_FLOAT_FMT_G(16)
-#endif
-
-#undef OUTPUT_FLOAT_FMT_G
-
-
 /* Define a macro to build code for write_float.  */
 
   /* Note: Before output_float is called, snprintf is used to print to buffer the
@@ -1003,19 +926,15 @@ __qmath_(quadmath_snprintf) (buffer, sizeof buffer
 \
 	DTOA ## y\
 \
-	if (f->format != FMT_G)\
-	  output_float (dtp, f, buffer, size, sign_bit, zero_flag, ndigits, \
-			edigits);\
-	else \
-	  output_float_FMT_G_ ## x (dtp, f, tmp, buffer, size, sign_bit, \
-				    zero_flag, ndigits, edigits, comp_d);\
+	output_float (dtp, f, buffer, size, sign_bit, zero_flag, ndigits, \
+		      edigits);\
 }\
 
 /* Output a real number according to its format.  */
 
 static void
 write_float (st_parameter_dt *dtp, const fnode *f, const char *source, \
-            int len, int comp_d)
+            int len)
 {
 
 #if defined(HAVE_GFC_REAL_16) || __LDBL_DIG__ > 18

^ permalink raw reply	[flat|nested] 8+ messages in thread

* Re: [patch, libgfortran] PR48906 Wrong rounding results with -m32
  2011-06-10 17:39 ` [patch, libgfortran] PR48906 " jerry DeLisle
@ 2011-06-11  8:24   ` Janne Blomqvist
  2011-06-11  9:13     ` Thomas Henlich
  0 siblings, 1 reply; 8+ messages in thread
From: Janne Blomqvist @ 2011-06-11  8:24 UTC (permalink / raw)
  To: jerry DeLisle; +Cc: gfortran, gcc patches, Thomas Henlich

On Fri, Jun 10, 2011 at 20:27, jerry DeLisle <jvdelisle@charter.net> wrote:
> On 06/03/2011 05:51 AM, jerry DeLisle wrote:
>>
>> Hi,
>>
>> The attached patch, which includes test cases, fixes this bug by
>> eliminating the
>> code which used floating point instructions to determine the 'r' value as
>> outlined in the Fortran standard under G formatting.
>>
>> Essentially, the code now examines the d and e values to determine the
>> number of
>> digits before and after the decimal point and whether or not to display
>> the 'E'
>> exponent symbol. Adjustments are made for various corner cases, including
>> when
>> rounding has resulted in a carry. (see PR for details of the trials)
>>
>> This patch is intrusive. It results in a minor performance improvement. It
>> eliminates a bit of code.
>>
>> Regression tested on x86-64.
>>
>> OK for trunk? then later back port to 4.6 after some proving time?
>
> Attached is updated patch based on comments from Thomas Henlich.  See my
> comment #33 and subsequent comments in the PR48906 explaining the issue with
> test case fmt_g0_6.f08 which is revised by the updated patch.
>
> Regression tested on X86-64.
>
> OK for trunk. (please ;) )

Index: gcc/testsuite/gfortran.dg/char4_iunit_1.f03
===================================================================
--- gcc/testsuite/gfortran.dg/char4_iunit_1.f03	(revision 174673)
+++ gcc/testsuite/gfortran.dg/char4_iunit_1.f03	(working copy)
@@ -24,7 +24,7 @@ program char4_iunit_1
   write(string, *) .true., .false. , .true.
   if (string .ne. 4_" T F T                                    ") call abort
   write(string, *) 1.2345e-06, 4.2846e+10_8
-  if (string .ne. 4_"   1.23450002E-06   42846000000.000000      ") call abort
+  if (string .ne. 4_"  1.234500019E-06   42846000000.000000      ") call abort
   write(string, *) nan, inf


I don't agree with this; with the patch we now output 10 significant
digits, whereas 9 is sufficient for a binary->ascii->binary roundtrip.
So please retain the "reduce d by one when E editing is used" thing
for list format and G0. This is just a side effect of using 1PGw.d
format for list format and G0 in order to avoid duplicating code, but
we don't need to follow this particular craziness that is due to how
the scale factor is specified in the standard.

Otherwise it looks nice. Good job!

FWIW, as this is quite tricky code and not a regression, IMHO we
should not backport it.

-- 
Janne Blomqvist

^ permalink raw reply	[flat|nested] 8+ messages in thread

* Re: [patch, libgfortran] PR48906 Wrong rounding results with -m32
  2011-06-11  8:24   ` Janne Blomqvist
@ 2011-06-11  9:13     ` Thomas Henlich
  2011-06-11 14:19       ` jerry DeLisle
  2011-06-14  6:14       ` jerry DeLisle
  0 siblings, 2 replies; 8+ messages in thread
From: Thomas Henlich @ 2011-06-11  9:13 UTC (permalink / raw)
  To: Janne Blomqvist; +Cc: jerry DeLisle, gfortran, gcc patches

> I don't agree with this; with the patch we now output 10 significant
> digits, whereas 9 is sufficient for a binary->ascii->binary roundtrip.
> So please retain the "reduce d by one when E editing is used" thing
> for list format and G0. This is just a side effect of using 1PGw.d
> format for list format and G0 in order to avoid duplicating code, but
> we don't need to follow this particular craziness that is due to how
> the scale factor is specified in the standard.

I hadn't noticed this, but I agree with Janne.

It should be easy to implement:

After the switch between F and E editing, we just need to shift the
decimal point and decrement the exponent. No new rounding is required,
because we keep the number of significant digits.

^ permalink raw reply	[flat|nested] 8+ messages in thread

* Re: [patch, libgfortran] PR48906 Wrong rounding results with -m32
  2011-06-11  9:13     ` Thomas Henlich
@ 2011-06-11 14:19       ` jerry DeLisle
  2011-06-11 16:05         ` Thomas Henlich
  2011-06-14  6:14       ` jerry DeLisle
  1 sibling, 1 reply; 8+ messages in thread
From: jerry DeLisle @ 2011-06-11 14:19 UTC (permalink / raw)
  To: Thomas Henlich; +Cc: Janne Blomqvist, gfortran, gcc patches

On 06/11/2011 12:23 AM, Thomas Henlich wrote:
>> I don't agree with this; with the patch we now output 10 significant
>> digits, whereas 9 is sufficient for a binary->ascii->binary roundtrip.
>> So please retain the "reduce d by one when E editing is used" thing
>> for list format and G0. This is just a side effect of using 1PGw.d
>> format for list format and G0 in order to avoid duplicating code, but
>> we don't need to follow this particular craziness that is due to how
>> the scale factor is specified in the standard.
>
> I hadn't noticed this, but I agree with Janne.
>
> It should be easy to implement:
>
> After the switch between F and E editing, we just need to shift the
> decimal point and decrement the exponent. No new rounding is required,
> because we keep the number of significant digits.
>

Our default formats for kind=4 are:

static void
set_fnode_default (st_parameter_dt *dtp, fnode *f, int length)
{
   f->format = FMT_G;
   switch (length)
     {
     case 4:
       f->u.real.w = 16;
       f->u.real.d = 9;
       f->u.real.e = 2;
       break;

This was established as solution to PR48488 where we had two choices for 
selecting the significant digits. Nine significant digits was established as a 
requirement to guarantee round trip in all cases. The char4_iunit_1.f03 test 
case was revised because after we corrected the formatting in PR48906, it 
started to fail and I observed the test case was looking for the wrong number of 
significant digits.

Based on this, I would suggest we leave it as I have it, which is correct.

Jerry

^ permalink raw reply	[flat|nested] 8+ messages in thread

* Re: [patch, libgfortran] PR48906 Wrong rounding results with -m32
  2011-06-11 14:19       ` jerry DeLisle
@ 2011-06-11 16:05         ` Thomas Henlich
  0 siblings, 0 replies; 8+ messages in thread
From: Thomas Henlich @ 2011-06-11 16:05 UTC (permalink / raw)
  To: jerry DeLisle; +Cc: Janne Blomqvist, gfortran, gcc patches

On Sat, Jun 11, 2011 at 14:41, jerry DeLisle <jvdelisle@charter.net> wrote:
> This was established as solution to PR48488 where we had two choices for
> selecting the significant digits. Nine significant digits was established as
> a requirement to guarantee round trip in all cases. The char4_iunit_1.f03
> test case was revised because after we corrected the formatting in PR48906,
> it started to fail and I observed the test case was looking for the wrong
> number of significant digits.
>
> Based on this, I would suggest we leave it as I have it, which is correct.

I'm afraid it's not.

1.23450002E-06 has nine significant digits. That's how it should be.

We don't want 1PG16.9E2 editing for list-directed and G0,
but G16.9E2 for the F editing range and 1PE16.8E2 editing for the E range.

This is to make sure the result always has nine significant digits,
whether in the F or E range.

^ permalink raw reply	[flat|nested] 8+ messages in thread

* Re: [patch, libgfortran] PR48906 Wrong rounding results with -m32
  2011-06-11  9:13     ` Thomas Henlich
  2011-06-11 14:19       ` jerry DeLisle
@ 2011-06-14  6:14       ` jerry DeLisle
  2011-06-14  8:30         ` Thomas Henlich
  1 sibling, 1 reply; 8+ messages in thread
From: jerry DeLisle @ 2011-06-14  6:14 UTC (permalink / raw)
  To: gfortran; +Cc: Thomas Henlich, Janne Blomqvist, gcc patches

[-- Attachment #1: Type: text/plain, Size: 1077 bytes --]

On 06/11/2011 12:23 AM, Thomas Henlich wrote:
>> I don't agree with this; with the patch we now output 10 significant
>> digits, whereas 9 is sufficient for a binary->ascii->binary roundtrip.
>> So please retain the "reduce d by one when E editing is used" thing
>> for list format and G0. This is just a side effect of using 1PGw.d
>> format for list format and G0 in order to avoid duplicating code, but
>> we don't need to follow this particular craziness that is due to how
>> the scale factor is specified in the standard.
>
> I hadn't noticed this, but I agree with Janne.
>
> It should be easy to implement:
>
> After the switch between F and E editing, we just need to shift the
> decimal point and decrement the exponent. No new rounding is required,
> because we keep the number of significant digits.
>

OK, after a little bit of experimentation, I have arrived at the updated patch 
attached.

This has been regression tested and passes all test cases I am aware of.  I also 
have included a new test case gcc/testsuite/gfortran.dg/fmt_g.f90.

OK for trunk?

Jerry

[-- Attachment #2: pr48906-final-3.diff --]
[-- Type: text/x-patch, Size: 22667 bytes --]

Index: gcc/testsuite/gfortran.dg/fmt_g.f90
===================================================================
--- gcc/testsuite/gfortran.dg/fmt_g.f90	(revision 0)
+++ gcc/testsuite/gfortran.dg/fmt_g.f90	(revision 0)
@@ -0,0 +1,88 @@
+! { dg-do run }
+! PR48906 Wrong rounding results with -m32
+! This test case presents various corner cases bumped into while reworking
+! handling of G formatting of reals without using floating point operations.
+  implicit none
+  integer, parameter :: RT = 8
+  character(28) :: buf
+  write(buf, "(rc,f11.2,4x,'<RC')") .995_RT
+  if (buf.ne."       0.99    <RC") call abort
+  write(buf, "(rc,g15.2,'<RC')") .995_RT
+  if (buf.ne."       0.99    <RC") call abort
+  write(buf, "(rd,f11.2,4x,'<RD')") .995_RT
+  if (buf.ne."       0.99    <RD") call abort
+  write(buf, "(rd,g15.2,'<RD')") .995_RT
+  if (buf.ne."       0.99    <RD") call abort
+  write(buf, "(ru,f11.1,4x,'<RU')") .995_RT
+  if (buf.ne."        1.0    <RU") call abort
+  write(buf, "(ru,g15.2,'<RU')") .995_RT
+  if (buf.ne."        1.0    <RU") call abort
+  write(buf,"(ru,g15.2)")  .099d0
+  if (buf.ne."       0.10    ") call abort
+  write(buf,"(rc,g15.1)")  .095d0
+  if (buf.ne."        0.1    ") call abort
+  write(buf,"(rc,g15.2)")  .0995d0
+  if (buf.ne."       0.10    ") call abort
+  write(buf,"(ru,g15.3)")  .0999d0
+  if (buf.ne."      0.100    ") call abort
+  write(buf,"(rc,g0.4,'<')")  .1234d-3
+  if (buf.ne.".1234E-003<") call abort
+  write(buf,"(rc,g0.4,'<')")  .1234d-2
+  if (buf.ne.".1234E-002<") call abort
+  write(buf,"(rc,g0.4,'<')")  .1234d-1
+  if (buf.ne.".1234E-001<") call abort
+  write(buf,"(rc,g0.4,'<')")  .1234d0
+  if (buf.ne.".1234<") call abort
+  write(buf,"(rc,g0.4,'<')")  .1234d1
+  if (buf.ne."1.234<") call abort
+  write(buf,"(rc,g0.4,'<')")  .1234d2
+  if (buf.ne."12.34<") call abort
+  write(buf,"(rc,g0.4,'<')")  .1234d3
+  if (buf.ne."123.4<") call abort
+  write(buf,"(rc,g0.4,'<')")  .1234d4
+  if (buf.ne."1234.<") call abort
+  write(buf,"(rc,g0.4,'<')")  .1234d5
+  if (buf.ne.".1234E+005<") call abort
+  write(buf,"(rc,g0.4,'<')")  .1234d6
+  if (buf.ne.".1234E+006<") call abort
+  write(buf,"(rc,g0.4,'<')")  .1234d7
+  if (buf.ne.".1234E+007<") call abort
+  write(buf,"(rc,1p,g15.4e2,'<')")  0.1234
+  if (buf.ne."     0.1234    <") call abort
+  write(buf,"(rc,g15.4e2,'<')")  0.1234
+  if (buf.ne."     0.1234    <") call abort
+  write(buf,"(rc,-1p,g15.4e2,'<')")  0.1234 
+  if (buf.ne."     0.1234    <") call abort
+  write(buf,"(rc,g10.2,'<')")  99.5
+  if (buf.ne."  0.10E+03<") call abort
+  write(buf,"(rc,g10.2,'<')")  995.
+  if (buf.ne."  0.10E+04<") call abort
+  write(buf,"(rc,g10.3,'<')")  999.5
+  if (buf.ne." 0.100E+04<") call abort
+  write(buf,"(rc,g10.3,'<')")  9995.
+  if (buf.ne." 0.100E+05<") call abort
+  write(buf,"(rc,e10.3,'<')")  9995.
+  if (buf.ne." 0.100E+05<") call abort
+  write(buf,"(g10.3,a)")  -0.0012346, "z"
+  if (buf.ne."-0.123E-02z") call abort
+  write(buf,"(g0.5,'<')")  -10000.
+  if (buf.ne."-10000.<") call abort
+  write(buf,"(g14.5e5,'<')")  -10000.
+  if (buf.ne."-10000.       <") call abort
+  write(buf,"(g15.5e5,'<')")  -10000.
+  if (buf.ne." -10000.       <") call abort
+  write(buf,"(g16.5e5,'<')")  -10000.
+  if (buf.ne."  -10000.       <") call abort
+  write(buf,"(rc,g11.1)")  2.
+  if (buf.ne."     2.    ") call abort
+  write(buf,"(rc,g11.2)")  20.
+  if (buf.ne."    20.    ") call abort
+  write(buf,"(rc,g11.3)")  200.
+  if (buf.ne."   200.    ") call abort
+  write(buf,"(rc,g11.4)")  2000.
+  if (buf.ne."  2000.    ") call abort
+  write(buf,"(rc,g11.2)")  .04
+  if (buf.ne."   0.40E-01") call abort
+  write(buf,"(rc,g11.3)")  .04
+  if (buf.ne."  0.400E-01") call abort
+end
Index: gcc/testsuite/gfortran.dg/fmt_g0_6.f08
===================================================================
--- gcc/testsuite/gfortran.dg/fmt_g0_6.f08	(revision 174673)
+++ gcc/testsuite/gfortran.dg/fmt_g0_6.f08	(working copy)
@@ -9,7 +9,7 @@ program test_g0fr
     
     call check_all(0.0_RT, 15, 2, 0)
     call check_all(0.991_RT, 15, 2, 0)
-    call check_all(0.995_RT, 15, 2, 0)
+    call check_all(0.99500000001_RT, 15, 2, 0)
     call check_all(0.996_RT, 15, 2, 0)
     call check_all(0.999_RT, 15, 2, 0)
 contains
Index: libgfortran/io/write.c
===================================================================
--- libgfortran/io/write.c	(revision 174673)
+++ libgfortran/io/write.c	(working copy)
@@ -1155,35 +1155,35 @@ write_z (st_parameter_dt *dtp, const fnode *f, con
 void
 write_d (st_parameter_dt *dtp, const fnode *f, const char *p, int len)
 {
-  write_float (dtp, f, p, len, 0);
+  write_float (dtp, f, p, len);
 }
 
 
 void
 write_e (st_parameter_dt *dtp, const fnode *f, const char *p, int len)
 {
-  write_float (dtp, f, p, len, 0);
+  write_float (dtp, f, p, len);
 }
 
 
 void
 write_f (st_parameter_dt *dtp, const fnode *f, const char *p, int len)
 {
-  write_float (dtp, f, p, len, 0);
+  write_float (dtp, f, p, len);
 }
 
 
 void
 write_en (st_parameter_dt *dtp, const fnode *f, const char *p, int len)
 {
-  write_float (dtp, f, p, len, 0);
+  write_float (dtp, f, p, len);
 }
 
 
 void
 write_es (st_parameter_dt *dtp, const fnode *f, const char *p, int len)
 {
-  write_float (dtp, f, p, len, 0);
+  write_float (dtp, f, p, len);
 }
 
 
@@ -1476,7 +1476,7 @@ write_real (st_parameter_dt *dtp, const char *sour
   int org_scale = dtp->u.p.scale_factor;
   dtp->u.p.scale_factor = 1;
   set_fnode_default (dtp, &f, length);
-  write_float (dtp, &f, source , length, 1);
+  write_float (dtp, &f, source , length);
   dtp->u.p.scale_factor = org_scale;
 }
 
@@ -1487,19 +1487,11 @@ void
 write_real_g0 (st_parameter_dt *dtp, const char *source, int length, int d)
 {
   fnode f;
-  int comp_d; 
   set_fnode_default (dtp, &f, length);
   if (d > 0)
     f.u.real.d = d;
-
-  /* Compensate for extra digits when using scale factor, d is not
-     specified, and the magnitude is such that E editing is used.  */
-  if (dtp->u.p.scale_factor > 0 && d == 0)
-    comp_d = 1;
-  else
-    comp_d = 0;
   dtp->u.p.g0_no_blanks = 1;
-  write_float (dtp, &f, source , length, comp_d);
+  write_float (dtp, &f, source , length);
   dtp->u.p.g0_no_blanks = 0;
 }
 
Index: libgfortran/io/write_float.def
===================================================================
--- libgfortran/io/write_float.def	(revision 174673)
+++ libgfortran/io/write_float.def	(working copy)
@@ -59,8 +59,22 @@ calculate_sign (st_parameter_dt *dtp, int negative
 }
 
 
-/* Output a real number according to its format which is FMT_G free.  */
+/* Output a real number according to its format.
 
+   What this does:
+
+	- Calculates the number of digits to emit before and after the
+	  decimal point.
+	- Performs the required rounding.
+	- Calculates the number of leading blanks to emit.
+	- Calculates the number of trailing blanks to emit.
+	- Calculates the exponent format if any.
+	- Determines if the sign should be displayed.
+	- Allocates the write buffer according to the width.
+	- Emits the formatted number according to the above.
+
+   buffer contains the digit string when we enter this function.  */
+
 static try
 output_float (st_parameter_dt *dtp, const fnode *f, char *buffer, size_t size, 
 	      int sign_bit, bool zero_flag, int ndigits, int edigits)
@@ -78,8 +92,10 @@ output_float (st_parameter_dt *dtp, const fnode *f
   int nafter;
   /* Number of zeros after the decimal point, whatever the precision.  */
   int nzero_real;
+  /* Flag for optional leading zero, if there is room.  */
   int leadzero;
-  int nblanks;
+  /* Number of leading and trailing blanks.  */
+  int lblanks, tblanks; 
   sign_t sign;
 
   ft = f->format;
@@ -87,7 +103,6 @@ output_float (st_parameter_dt *dtp, const fnode *f
   d = f->u.real.d;
   p = dtp->u.p.scale_factor;
 
-  rchar = '5';
   nzero_real = -1;
 
   /* We should always know the field width and precision.  */
@@ -142,19 +157,57 @@ output_float (st_parameter_dt *dtp, const fnode *f
       expchar = 0;
       break;
 
+    case FMT_G:
+      nbefore = 0;
+      nzero = 0;
+      nafter = d;
+      expchar = 0;
+
+      if (e == d)
+	{
+	  if (e > 0)
+	    {
+	      nbefore = e;
+	      nafter = d > nbefore ? d - nbefore : 0;
+	    }
+	  else
+	    expchar = 'E';
+	  break;
+	}
+      else if (e < d && e >= 0)
+	{
+	  nbefore = e;
+	  nafter = d > nbefore ? d - nbefore : 0;
+	  break;
+	}
+      else if (e >= -d && e < 0)
+	{
+	  if (f->u.real.e == -1)
+	    nafter = d;
+	  else if (p == 1)
+	    {
+	      p = 0;
+	      e--;
+	      nbefore++;
+	      nafter--;
+	      expchar = 'E';
+	      break;
+	    }
+	}
+    /* Fall through.  */
     case FMT_E:
     case FMT_D:
-      i = dtp->u.p.scale_factor;
       if (d <= 0 && p == 0)
 	{
 	  generate_error (&dtp->common, LIBERROR_FORMAT, "Precision not "
-			  "greater than zero in format specifier 'E' or 'D'");
+			  "greater than zero in format specifier 'D', 'E', "
+			  "or 'G'");
 	  return FAILURE;
 	}
       if (p <= -d || p >= d + 2)
 	{
 	  generate_error (&dtp->common, LIBERROR_FORMAT, "Scale factor "
-			  "out of range in format specifier 'E' or 'D'");
+			  "out of range in format specifier 'E', 'D', or 'G'");
 	  return FAILURE;
 	}
 
@@ -179,10 +232,11 @@ output_float (st_parameter_dt *dtp, const fnode *f
 	  nafter = d;
 	}
 
-      if (ft == FMT_E)
+      if (ft == FMT_E || ft == FMT_G)
 	expchar = 'E';
       else
 	expchar = 'D';
+      
       break;
 
     case FMT_EN:
@@ -224,6 +278,7 @@ output_float (st_parameter_dt *dtp, const fnode *f
 
   /* Round the value.  The value being rounded is an unsigned magnitude.
      The ROUND_COMPATIBLE is rounding away from zero when there is a tie.  */
+  rchar = '5';
   switch (dtp->u.p.current_unit->round_status)
     {
       case ROUND_ZERO: /* Do nothing and truncation occurs.  */
@@ -275,6 +330,7 @@ output_float (st_parameter_dt *dtp, const fnode *f
   rchar = '0';
   if (w > 0 && d == 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++)
     {
@@ -284,7 +340,7 @@ output_float (st_parameter_dt *dtp, const fnode *f
   goto skip;
     
   do_rnd:
- 
+
   if (nbefore + nafter == 0)
     {
       ndigits = 0;
@@ -321,6 +377,7 @@ output_float (st_parameter_dt *dtp, const fnode *f
 	         zero.  */
 	      digits--;
 	      digits[0] = '1';
+
 	      if (ft == FMT_F)
 		{
 		  if (nzero > 0)
@@ -331,6 +388,27 @@ output_float (st_parameter_dt *dtp, const fnode *f
 		  else
 		    nbefore++;
 		}
+	      else if (ft == FMT_G)
+		{
+		  if (e < d && e >= 0)
+		    {
+		      nbefore++;
+		      nafter--;
+		    }
+		  else if (e == d)
+		    {
+		      nbefore -= d;
+		      nafter += d;
+		      e++;
+		      expchar = 'E';
+		    }
+		  else
+		    {
+		      e++;
+		      if (e == 0)
+		        expchar = 0;
+		    }
+		}
 	      else if (ft == FMT_EN)
 		{
 		  nbefore++;
@@ -346,11 +424,68 @@ output_float (st_parameter_dt *dtp, const fnode *f
 	}
     }
 
+  /* Scan the digits string and count the number of zeros.  If we make it
+     all the way through the loop, we know the value is zero after the
+     rounding completed above. To format properly, we need to know if the
+     rounded result is zero and if so, we set the zero_flag which may have
+     been already set for actual zero.  */
+  for (i = 0; i < ndigits; i++)
+    {
+      if (digits[i] != '0')
+	break;
+    }
+  if (i == ndigits)
+    {
+      zero_flag = true;
+      /* The output is zero, so set the sign according to the sign bit unless
+	 -fno-sign-zero was specified.  */
+      if (compile_options.sign_zero == 1)
+        sign = calculate_sign (dtp, sign_bit);
+      else
+	sign = calculate_sign (dtp, 0);
+    }
+
   skip:
 
+  /* Pick a field size if none was specified, taking into account small
+     values that may have been rounded to zero.  */
+  if (w <= 0)
+    {
+      if (zero_flag)
+	w = d + (sign != S_NONE ? 2 : 1) + (d == 0 ? 1 : 0);
+      else
+	{
+	  w = nbefore + nzero + nafter + (sign != S_NONE ? 2 : 1);
+	  w = w == 1 ? 2 : w;
+	}
+    }
+  
+  /* Figure out the leading and trailing blanks.  */
+  lblanks = tblanks = 0;
+
+  /* G formatting is a special case of several things.  */
+  if (ft == FMT_G)
+    {
+      if (zero_flag && nbefore == 0 && nafter > 0)
+	{
+	  nbefore++;
+	  nafter--;
+	}
+      tblanks = edigits = f->u.real.e == -1 ? 4 : f->u.real.e + 2;
+      lblanks = w - (nbefore + nzero + nafter + 1 + tblanks +
+		     (sign != S_NONE ? 1 : 0));
+    }
+  else
+    {
+      lblanks = w - (nbefore + nzero + nafter + 1 + e);
+      if (sign != S_NONE)
+	lblanks--;
+    }
+
   /* Calculate the format of the exponent field.  */
   if (expchar)
     {
+      tblanks = 0;
       edigits = 1;
       for (i = abs (e); i >= 10; i /= 10)
 	edigits++;
@@ -379,76 +514,57 @@ output_float (st_parameter_dt *dtp, const fnode *f
   else
     edigits = 0;
 
-  /* Scan the digits string and count the number of zeros.  If we make it
-     all the way through the loop, we know the value is zero after the
-     rounding completed above.  */
-  for (i = 0; i < ndigits; i++)
+  if (ft == FMT_G)
     {
-      if (digits[i] != '0')
-	break;
-    }
-
-  /* To format properly, we need to know if the rounded result is zero and if
-     so, we set the zero_flag which may have been already set for
-     actual zero.  */
-  if (i == ndigits)
-    {
-      zero_flag = true;
-      /* The output is zero, so set the sign according to the sign bit unless
-	 -fno-sign-zero was specified.  */
-      if (compile_options.sign_zero == 1)
-        sign = calculate_sign (dtp, sign_bit);
-      else
-	sign = calculate_sign (dtp, 0);
-    }
-
-  /* Pick a field size if none was specified, taking into account small
-     values that may have been rounded to zero.  */
-  if (w <= 0)
-    {
-      if (zero_flag)
-	w = d + (sign != S_NONE ? 2 : 1) + (d == 0 ? 1 : 0);
-      else
+      if (dtp->u.p.g0_no_blanks)
 	{
-	  w = nbefore + nzero + nafter + (sign != S_NONE ? 2 : 1);
-	  w = w == 1 ? 2 : w;
+	  w = w - tblanks - lblanks;
+	  lblanks = tblanks = 0;
 	}
+      else if (expchar == 0)
+	{
+	  lblanks = w - (nbefore + nzero + nafter + 1 + tblanks +
+			 (sign != S_NONE ? 1 : 0));
+	  if (dtp->u.p.no_leading_blank)
+	    {
+	      tblanks += lblanks;
+	      lblanks = 0; 
+	    }
+	}
     }
-  
-  /* Work out how much padding is needed.  */
-  nblanks = w - (nbefore + nzero + nafter + edigits + 1);
-  if (sign != S_NONE)
-    nblanks--;
+  else
+    lblanks = w - (nbefore + nzero + nafter + 1 + edigits +
+		 (sign != S_NONE ? 1 : 0));
 
-  if (dtp->u.p.g0_no_blanks)
-    {
-      w -= nblanks;
-      nblanks = 0;
-    }
-
   /* Create the ouput buffer.  */
-  out = write_block (dtp, w);
+  if ((ft == FMT_G && dtp->u.p.g0_no_blanks)
+      || (ft == FMT_F && f->u.real.w == 0))
+    out = write_block (dtp, w);
+  else
+    out = write_block (dtp, f->u.real.w);
+
   if (out == NULL)
     return FAILURE;
 
   /* Check the value fits in the specified field width.  */
-  if (nblanks < 0 || edigits == -1 || w == 1 || (w == 2 && sign != S_NONE))
+  if (lblanks < 0 || tblanks < 0 || edigits == -1 || w == 1
+      || (w == 2 && sign != S_NONE))
     {
       if (unlikely (is_char4_unit (dtp)))
 	{
 	  gfc_char4_t *out4 = (gfc_char4_t *) out;
-	  memset4 (out4, '*', w);
+	  memset4 (out4, '*', f->u.real.w);
 	  return FAILURE;
 	}
-      star_fill (out, w);
+      star_fill (out, f->u.real.w);
       return FAILURE;
     }
 
   /* See if we have space for a zero before the decimal point.  */
-  if (nbefore == 0 && nblanks > 0)
+  if (nbefore == 0 && lblanks > 0)
     {
       leadzero = 1;
-      nblanks--;
+      lblanks--;
     }
   else
     leadzero = 0;
@@ -461,10 +577,10 @@ output_float (st_parameter_dt *dtp, const fnode *f
       gfc_char4_t *out4 = (gfc_char4_t *) out;
       /* Pad to full field width.  */
 
-      if ( ( nblanks > 0 ) && !dtp->u.p.no_leading_blank)
+      if ((lblanks > 0) && !dtp->u.p.no_leading_blank)
 	{
-	  memset4 (out4, ' ', nblanks);
-	  out4 += nblanks;
+	  memset4 (out4, ' ', lblanks);
+	  out4 += lblanks;
 	}
 
       /* Output the initial sign (if any).  */
@@ -539,21 +655,17 @@ output_float (st_parameter_dt *dtp, const fnode *f
 	  memcpy4 (out4, buffer, edigits);
 	}
 
-      if (dtp->u.p.no_leading_blank)
-	{
-	  out4 += edigits;
-	  memset4 (out4, ' ' , nblanks);
-	  dtp->u.p.no_leading_blank = 0;
-	}
+      /* Emit trailing blanks if needed.  */
+      memset4 (out4, ' ' , tblanks);
       return SUCCESS;
     } /* End of character(kind=4) internal unit code.  */
 
-  /* Pad to full field width.  */
+  /* Pad leading to full field width.  */
 
-  if ( ( nblanks > 0 ) && !dtp->u.p.no_leading_blank)
+  if (lblanks > 0 && !dtp->u.p.no_leading_blank)
     {
-      memset (out, ' ', nblanks);
-      out += nblanks;
+      memset (out, ' ', lblanks);
+      out += lblanks;
     }
 
   /* Output the initial sign (if any).  */
@@ -627,13 +739,10 @@ output_float (st_parameter_dt *dtp, const fnode *f
       memcpy (out, buffer, edigits);
     }
 
-  if (dtp->u.p.no_leading_blank)
-    {
-      out += edigits;
-      memset( out , ' ' , nblanks );
-      dtp->u.p.no_leading_blank = 0;
-    }
-
+  /* Emit trailing blanks if needed.  */
+  if (tblanks)
+    memset( out , ' ' , tblanks );
+  
 #undef STR
 #undef STR1
 #undef MIN_FIELD_WIDTH
@@ -770,185 +879,6 @@ write_infnan (st_parameter_dt *dtp, const fnode *f
     }
 }
 
-
-/* Returns the value of 10**d.  */
-
-#define CALCULATE_EXP(x) \
-inline static GFC_REAL_ ## x \
-calculate_exp_ ## x  (int d)\
-{\
-  int i;\
-  GFC_REAL_ ## x r = 1.0;\
-  for (i = 0; i< (d >= 0 ? d : -d); i++)\
-    r *= 10;\
-  r = (d >= 0) ? r : 1.0 / r;\
-  return r;\
-}
-
-CALCULATE_EXP(4)
-
-CALCULATE_EXP(8)
-
-#ifdef HAVE_GFC_REAL_10
-CALCULATE_EXP(10)
-#endif
-
-#ifdef HAVE_GFC_REAL_16
-CALCULATE_EXP(16)
-#endif
-#undef CALCULATE_EXP
-
-/* Generate corresponding I/O format for FMT_G and output.
-   The rules to translate FMT_G to FMT_E or FMT_F from DEC fortran
-   LRM (table 11-2, Chapter 11, "I/O Formatting", P11-25) is:
-
-   Data Magnitude                              Equivalent Conversion
-   0< m < 0.1-0.5*10**(-d-1)                   Ew.d[Ee]
-   m = 0                                       F(w-n).(d-1), n' '
-   0.1-0.5*10**(-d-1)<= m < 1-0.5*10**(-d)     F(w-n).d, n' '
-   1-0.5*10**(-d)<= m < 10-0.5*10**(-d+1)      F(w-n).(d-1), n' '
-   10-0.5*10**(-d+1)<= m < 100-0.5*10**(-d+2)  F(w-n).(d-2), n' '
-   ................                           ..........
-   10**(d-1)-0.5*10**(-1)<= m <10**d-0.5       F(w-n).0,n(' ')
-   m >= 10**d-0.5                              Ew.d[Ee]
-
-   notes: for Gw.d ,  n' ' means 4 blanks
-	  for Gw.dEe, n' ' means e+2 blanks
-	  for rounding modes adjustment, r, See Fortran F2008 10.7.5.2.2
-	  the asm volatile is required for 32-bit x86 platforms.  */
-
-#define OUTPUT_FLOAT_FMT_G(x) \
-static void \
-output_float_FMT_G_ ## x (st_parameter_dt *dtp, const fnode *f, \
-		      GFC_REAL_ ## x m, char *buffer, size_t size, \
-		      int sign_bit, bool zero_flag, int ndigits, \
-                      int edigits, int comp_d) \
-{ \
-  int e = f->u.real.e;\
-  int d = f->u.real.d;\
-  int w = f->u.real.w;\
-  fnode *newf;\
-  GFC_REAL_ ## x rexp_d, r = 0.5;\
-  int low, high, mid;\
-  int ubound, lbound;\
-  char *p, pad = ' ';\
-  int save_scale_factor, nb = 0;\
-  try result;\
-\
-  save_scale_factor = dtp->u.p.scale_factor;\
-  newf = (fnode *) get_mem (sizeof (fnode));\
-\
-  switch (dtp->u.p.current_unit->round_status)\
-    {\
-      case ROUND_ZERO:\
-	r = sign_bit ? 1.0 : 0.0;\
-	break;\
-      case ROUND_UP:\
-	r = 1.0;\
-	break;\
-      case ROUND_DOWN:\
-	r = 0.0;\
-	break;\
-      default:\
-	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)))\
-      || ((m == 0.0) && !(compile_options.allow_std\
-			  & (GFC_STD_F2003 | GFC_STD_F2008))))\
-    { \
-      newf->format = FMT_E;\
-      newf->u.real.w = w;\
-      newf->u.real.d = d - comp_d;\
-      newf->u.real.e = e;\
-      nb = 0;\
-      goto finish;\
-    }\
-\
-  mid = 0;\
-  low = 0;\
-  high = d + 1;\
-  lbound = 0;\
-  ubound = d + 1;\
-\
-  while (low <= high)\
-    { \
-      volatile GFC_REAL_ ## x temp;\
-      mid = (low + high) / 2;\
-\
-      temp = (calculate_exp_ ## x (mid - 1) * (1 - r * rexp_d));\
-\
-      if (m < temp)\
-        { \
-          ubound = mid;\
-          if (ubound == lbound + 1)\
-            break;\
-          high = mid - 1;\
-        }\
-      else if (m > temp)\
-        { \
-          lbound = mid;\
-          if (ubound == lbound + 1)\
-            { \
-              mid ++;\
-              break;\
-            }\
-          low = mid + 1;\
-        }\
-      else\
-	{\
-	  mid++;\
-	  break;\
-	}\
-    }\
-\
-  nb = e <= 0 ? 4 : e + 2;\
-  nb = nb >= w ? w - 1 : nb;\
-  newf->format = FMT_F;\
-  newf->u.real.w = w - nb;\
-  newf->u.real.d = m == 0.0 ? d - 1 : -(mid - d - 1) ;\
-  dtp->u.p.scale_factor = 0;\
-\
- finish:\
-  result = output_float (dtp, newf, buffer, size, sign_bit, zero_flag, \
-			 ndigits, edigits);\
-  dtp->u.p.scale_factor = save_scale_factor;\
-\
-  free (newf);\
-\
-  if (nb > 0 && !dtp->u.p.g0_no_blanks)\
-    {\
-      p = write_block (dtp, nb);\
-      if (p == NULL)\
-	return;\
-      if (result == FAILURE)\
-        pad = '*';\
-      if (unlikely (is_char4_unit (dtp)))\
-	{\
-	  gfc_char4_t *p4 = (gfc_char4_t *) p;\
-	  memset4 (p4, pad, nb);\
-	}\
-      else \
-	memset (p, pad, nb);\
-    }\
-}\
-
-OUTPUT_FLOAT_FMT_G(4)
-
-OUTPUT_FLOAT_FMT_G(8)
-
-#ifdef HAVE_GFC_REAL_10
-OUTPUT_FLOAT_FMT_G(10)
-#endif
-
-#ifdef HAVE_GFC_REAL_16
-OUTPUT_FLOAT_FMT_G(16)
-#endif
-
-#undef OUTPUT_FLOAT_FMT_G
-
-
 /* Define a macro to build code for write_float.  */
 
   /* Note: Before output_float is called, snprintf is used to print to buffer the
@@ -1003,19 +933,15 @@ __qmath_(quadmath_snprintf) (buffer, sizeof buffer
 \
 	DTOA ## y\
 \
-	if (f->format != FMT_G)\
-	  output_float (dtp, f, buffer, size, sign_bit, zero_flag, ndigits, \
-			edigits);\
-	else \
-	  output_float_FMT_G_ ## x (dtp, f, tmp, buffer, size, sign_bit, \
-				    zero_flag, ndigits, edigits, comp_d);\
+	output_float (dtp, f, buffer, size, sign_bit, zero_flag, ndigits, \
+		      edigits);\
 }\
 
 /* Output a real number according to its format.  */
 
 static void
 write_float (st_parameter_dt *dtp, const fnode *f, const char *source, \
-            int len, int comp_d)
+            int len)
 {
 
 #if defined(HAVE_GFC_REAL_16) || __LDBL_DIG__ > 18

^ permalink raw reply	[flat|nested] 8+ messages in thread

* Re: [patch, libgfortran] PR48906 Wrong rounding results with -m32
  2011-06-14  6:14       ` jerry DeLisle
@ 2011-06-14  8:30         ` Thomas Henlich
  0 siblings, 0 replies; 8+ messages in thread
From: Thomas Henlich @ 2011-06-14  8:30 UTC (permalink / raw)
  To: jerry DeLisle; +Cc: gfortran, Janne Blomqvist, gcc patches

On Tue, Jun 14, 2011 at 06:51, jerry DeLisle <jvdelisle@charter.net> wrote:
>> It should be easy to implement:
>>
>> After the switch between F and E editing, we just need to shift the
>> decimal point and decrement the exponent. No new rounding is required,
>> because we keep the number of significant digits.
>>
>
> OK, after a little bit of experimentation, I have arrived at the updated
> patch attached.
>
> This has been regression tested and passes all test cases I am aware of.  I
> also have included a new test case gcc/testsuite/gfortran.dg/fmt_g.f90.
>
> OK for trunk?

I have reviewed your patch, and I noticed that you placed the
digit-shifting code quite at the top of output_float(), where the
final value of e is not even known. Due to rounding, e can be modified
after this point, so your code will generate invalid output in some
cases, for example:

print "(-2PG0)", nearest(0.1d0, -1.0d0) ! 1.0000000000000000E+001
expected .0099999999999999992E+001

Please put the code where at belongs, after the switch between F and E
editing (based on the final value of e).

The same applies to the scale factor in general, e.g.

print "(-2pg12.3)", 0.096    ! 1.00E+01 expected 0.001E+02
print "(-1pg12.3)", 0.0996   ! 1.00E+00 expected 0.010E+01
print "(-2pg12.3)", 0.09996  ! 1.00E+01 expected 0.100
print "(-1pg12.3)", 0.09996  ! 1.00E+00 expected 0.100
print "(1pg12.3)",  0.099996 ! 1.000E-01 expected 0.100
print "(2pg12.3)",  0.099996 ! 10.00E-02 expected 0.100
print "(-2pg12.3)",  999.6  ! 0.100E+04 expected 0.001E+06
print "(-1pg12.3)",  999.6  ! 0.100E+04 expected 0.010E+05
print "(1pg12.3)",  999.6  ! 0.100E+04 expected 9.996E+02
print "(2pg12.3)",  999.6  ! 0.100E+04 expected 99.96E+01

Please revise your code to fix this. A working approach I have outlined in
http://gcc.gnu.org/bugzilla/show_bug.cgi?id=48906#c28
and an (alpha) implementation is here:
http://gcc.gnu.org/bugzilla/show_bug.cgi?id=48906#c31

Thomas

^ permalink raw reply	[flat|nested] 8+ messages in thread

end of thread, other threads:[~2011-06-14  7:19 UTC | newest]

Thread overview: 8+ messages (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
2011-06-03 12:51 [patch, libgfortran] Wrong rounding results with -m32 jerry DeLisle
2011-06-10 17:39 ` [patch, libgfortran] PR48906 " jerry DeLisle
2011-06-11  8:24   ` Janne Blomqvist
2011-06-11  9:13     ` Thomas Henlich
2011-06-11 14:19       ` jerry DeLisle
2011-06-11 16:05         ` Thomas Henlich
2011-06-14  6:14       ` jerry DeLisle
2011-06-14  8:30         ` Thomas Henlich

This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).