re PR fortran/18218 (Miscompare in sixtrack benchmark caused by loss of precision)

PR fortran/18218
	* configure.ac: Check for strtof.
	* configure: Regenerate.
	* config.h.in: Regenerate.
	* io/read.c (convert_real): Use strtof if available.
	(convert_precision_real): Remove.
	(read_f): Avoid poor exponentiation algorithm.
gcc/testsuite/
	* gfortran.dg/list_read.c: New test.

From-SVN: r90382

gcc/testsuite/ChangeLog

@@ -1,3 +1,8 @@
+2004-11-10  Paul Brook  <paul@codesourcery.com>
+
+	PR fortran/18218
+	* gfortran.dg/list_read.c: New test.
+
 2004-11-09  Joseph S. Myers  <joseph@codesourcery.com>
 
 	PR c/18322

gcc/testsuite/gfortran.dg/read_float_1.f90

@@ -0,0 +1,12 @@
+! { dg-do run }
+! PR18218
+! The IO library has an algorithm that involved repeated multiplication by 10,
+! resulting in introducing large cumulative floating point errors.
+program foo
+  character*20 s
+  real*8 d
+  s = "-.18774312893273    "
+  read(unit=s, fmt='(g20.14)') d
+  if (d + 0.18774312893273d0 .gt. 1d-13) call abort
+end program
+

libgfortran/ChangeLog

@@ -1,3 +1,13 @@
+2004-11-10  Paul Brook  <paul@codesourcery.com>
+
+	PR fortran/18218
+	* configure.ac: Check for strtof.
+	* configure: Regenerate.
+	* config.h.in: Regenerate.
+	* io/read.c (convert_real): Use strtof if available.
+	(convert_precision_real): Remove.
+	(read_f): Avoid poor exponentiation algorithm.
+
 2004-11-05  Andreas Schwab  <schwab@suse.de>
 
 	* configure.ac: Use AC_PROG_FC, FC and FCFLAGS instead of

libgfortran/config.h.in

@@ -156,6 +156,9 @@
 /* Define to 1 if you have the <string.h> header file. */
 #undef HAVE_STRING_H
 
+/* Define to 1 if you have the `strtof' function. */
+#undef HAVE_STRTOF
+
 /* Define to 1 if you have the <sys/mman.h> header file. */
 #undef HAVE_SYS_MMAN_H
 

libgfortran/configure

@@ -6821,7 +6821,8 @@ fi
 
 
 
-for ac_func in getrusage times mkstemp
+
+for ac_func in getrusage times mkstemp strtof
 do
 as_ac_var=`echo "ac_cv_func_$ac_func" | $as_tr_sh`
 echo "$as_me:$LINENO: checking for $ac_func" >&5

libgfortran/configure.ac

@@ -159,7 +159,7 @@ AC_CHECK_HEADER([complex.h],[AC_DEFINE([HAVE_COMPLEX_H], [1], [complex.h exists]
 AC_CHECK_LIB([m],[csin],[need_math="no"],[need_math="yes"])
 
 # Check for library functions.
-AC_CHECK_FUNCS(getrusage times mkstemp)
+AC_CHECK_FUNCS(getrusage times mkstemp strtof)
 
 # Check libc for getgid, getpid, getuid
 AC_CHECK_LIB([c],[getgid],[AC_DEFINE([HAVE_GETGID],[1],[libc includes getgid])])

libgfortran/io/read.c

@@ -24,6 +24,7 @@ Boston, MA 02111-1307, USA.  */
 #include <errno.h>
 #include <ctype.h>
 #include <stdlib.h>
+#include <stdio.h>
 #include "libgfortran.h"
 #include "io.h"
 
@@ -89,8 +90,7 @@ max_value (int length, int signed_flag)
 /* convert_real()-- Convert a character representation of a floating
  * point number to the machine number.  Returns nonzero if there is a
  * range problem during conversion.  TODO: handle not-a-numbers and
- * infinities.  Handling of kind 4 is probably wrong because of double
+ * infinities.  */
- * rounding. */
 
 int
 convert_real (void *dest, const char *buffer, int length)
@@ -101,13 +101,18 @@ convert_real (void *dest, const char *buffer, int length)
   switch (length)
     {
     case 4:
-      *((float *) dest) = (float) strtod (buffer, NULL);
+      *((float *) dest) =
+#if defined(HAVE_STRTOF)
+	strtof (buffer, NULL);
+#else
+	(float) strtod (buffer, NULL);
+#endif
       break;
     case 8:
       *((double *) dest) = strtod (buffer, NULL);
       break;
     default:
-      internal_error ("Bad real number kind");
+      internal_error ("Unsupported real kind during IO");
     }
 
   if (errno != 0)
@@ -120,114 +125,6 @@ convert_real (void *dest, const char *buffer, int length)
   return 0;
 }
 
-static int
-convert_precision_real (void *dest, int sign,
-                       char *buffer, int length, int exponent)
-{
-  int w, new_dp_pos, i, slen, k, dp;
-  char * p, c;
-  double fval;
-  float tf;
-
-  fval =0.0;
-  tf = 0.0;
-  dp = 0;
-  new_dp_pos = 0;
-
-  slen = strlen (buffer);
-  w = slen;
-  p = buffer;
-
-/*  for (i = w - 1; i > 0; i --)
-    {
-       if (buffer[i] == '0' || buffer[i] == 0)
-         buffer[i] = 0;
-       else
-         break;
-    }
-*/
-  for (i = 0; i < w; i++)
-    {
-       if (buffer[i] == '.')
-         break;
-    }
-
-  new_dp_pos = i;
-  new_dp_pos += exponent;
-
-  while (w > 0)
-    {
-      c = *p;
-      switch (c)
-        {
-        case '0':
-        case '1':
-        case '2':
-        case '3':
-        case '4':
-        case '5':
-        case '6':
-        case '7':
-        case '8':
-        case '9':
-          fval = fval * 10.0 + c - '0';
-          p++;
-          w--;
-          break;
-
-        case '.':
-          dp = 1;
-          p++;
-          w--;
-          break;
-
-       default:
-          p++;
-          w--;
-          break;
-     }
-  }
-
-  if (sign)
-    fval = - fval;
-
-  i = new_dp_pos - slen + dp;
-  k = abs(i);
-  tf = 1.0;
-
-  while (k > 0)
-    {
-       tf *= 10.0 ;
-       k -- ;
-    }
-
-  if (fval != 0.0)
-    {
-       if (i < 0)
-         {
-           fval = fval / tf;
-         }
-        else
-         {
-           fval = fval * tf;
-         }
-    }
-
-  switch (length)
-    {
-    case 4:
-      *((float *) dest) = (float)fval;
-      break;
-    case 8:
-      *((double *) dest) = fval;
-      break;
-    default:
-      internal_error ("Bad real number kind");
-    }
-
-  return 0;
-}
-
 
 /* read_l()-- Read a logical value */
 
@@ -576,19 +473,23 @@ overflow:
 
 
 /* read_f()-- Read a floating point number with F-style editing, which
- * is what all of the other floating point descriptors behave as.  The
+   is what all of the other floating point descriptors behave as.  The
- * tricky part is that optional spaces are allowed after an E or D,
+   tricky part is that optional spaces are allowed after an E or D,
- * and the implicit decimal point if a decimal point is not present in
+   and the implicit decimal point if a decimal point is not present in
- * the input. */
+   the input.  */
 
 void
 read_f (fnode * f, char *dest, int length)
 {
   int w, seen_dp, exponent;
   int exponent_sign, val_sign;
-  char *p, *buffer, *n;
+  int ndigits;
+  int edigits;
+  int i;
+  char *p, *buffer;
+  char *digits;
 
-  val_sign = 0;
+  val_sign = 1;
   seen_dp = 0;
   w = f->u.w;
   p = read_block (&w);
@@ -601,32 +502,26 @@ read_f (fnode * f, char *dest, int length)
       switch (length)
 	{
 	case 4:
-	  *((float *) dest) = 0.0;
+	  *((float *) dest) = 0.0f;
 	  break;
 
 	case 8:
 	  *((double *) dest) = 0.0;
 	  break;
+
+	default:
+	  internal_error ("Unsupported real kind during IO");
 	}
 
       return;
     }
 
-  if (w + 2 < SCRATCH_SIZE)
-    buffer = scratch;
-  else
-    buffer = get_mem (w + 2);
-
-  memset(buffer, 0, w + 2);
-
-  n = buffer;
-
   /* Optional sign */
 
   if (*p == '-' || *p == '+')
     {
       if (*p == '-')
-        val_sign = 1;
+        val_sign = -1;
       p++;
 
       if (--w == 0)
@@ -640,10 +535,21 @@ read_f (fnode * f, char *dest, int length)
   if (!isdigit (*p) && *p != '.')
     goto bad_float;
 
+  /* Remember the position of the first digit.  */
+  digits = p;
+  ndigits = 0;
+
+  /* Scan through the string to find the exponent.  */
   while (w > 0)
     {
       switch (*p)
 	{
+	case '.':
+	  if (seen_dp)
+	    goto bad_float;
+	  seen_dp = 1;
+	  /* Fall through */
+
 	case '0':
 	case '1':
 	case '2':
@@ -654,23 +560,9 @@ read_f (fnode * f, char *dest, int length)
 	case '7':
 	case '8':
 	case '9':
-	  *n++ = *p++;
-	  w--;
-	  break;
-
-	case '.':
-	  if (seen_dp)
-	    goto bad_float;
-	  seen_dp = 1;
-
-	  *n++ = *p++;
-	  w--;
-	  break;
-
 	case ' ':
-	  if (g.blank_status == BLANK_ZERO)
+	  ndigits++;
-	    *n++ = '0';
+	  *p++;
-	  p++;
 	  w--;
 	  break;
 
@@ -732,8 +624,8 @@ exp1:
     goto bad_float;
 
 /* At this point a digit string is required.  We calculate the value
- * of the exponent in order to take account of the scale factor and
+   of the exponent in order to take account of the scale factor and
- * the d parameter before explict conversion takes place. */
+   the d parameter before explict conversion takes place. */
 
 exp2:
   if (!isdigit (*p))
@@ -746,9 +638,6 @@ exp2:
   while (w > 0 && isdigit (*p))
     {
       exponent = 10 * exponent + *p - '0';
-      if (exponent > 999999)
-	goto bad_float;
-
       p++;
       w--;
     }
@@ -766,14 +655,56 @@ exp2:
   exponent = exponent * exponent_sign;
 
 done:
+  /* Use the precision specified in the format if no decimal point has been
+     seen.  */
   if (!seen_dp)
     exponent -= f->u.real.d;
 
-  /* The number is syntactically correct and ready for conversion.
+  if (exponent > 0)
-   * The only thing that can go wrong at this point is overflow or
+    {
-   * underflow. */
+      edigits = 2;
+      i = exponent;
+    }
+  else
+    {
+      edigits = 3;
+      i = -exponent;
+    }
 
-  convert_precision_real (dest, val_sign, buffer, length, exponent);
+  while (i >= 10)
+    {
+      i /= 10;
+      edigits++;
+    }
+
+  i = ndigits + edigits + 1;
+  if (val_sign < 0)
+    i++;
+
+  if (i < SCRATCH_SIZE) 
+    buffer = scratch;
+  else
+    buffer = get_mem (i);
+
+  /* Reformat the string into a temporary buffer.  As we're using atof it's
+     easiest to just leave the dcimal point in place.  */
+  p = buffer;
+  if (val_sign < 0)
+    *(p++) = '-';
+  for (; ndigits > 0; ndigits--)
+    {
+      if (*digits == ' ' && g.blank_status == BLANK_ZERO)
+	*p = '0';
+      else
+	*p = *digits;
+      p++;
+      digits++;
+    }
+  *(p++) = 'e';
+  sprintf (p, "%d", exponent);
+
+  /* Do the actual conversion.  */
+  string_to_real (dest, buffer, length);
 
   if (buffer != scratch)
      free_mem (buffer);