This is the mail archive of the gcc-patches@gcc.gnu.org mailing list for the GCC project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

[gfortran] PR18218 - loss of floating point precision.


The formatted real input routines used a very poor exponentiation algorithm 
(successive multiplication by 10). This can introduce large rounding errors.
This wasn't so obvious on x86 as the excess FP precision was working in out 
favour.  I haven't yet checked if PR18218 is actually fixed, but this was 
definitely a contributory factor.

The patch below changes the IO library to use strto{f,d} to do the hard work.

Tested on i686-linux and powerpc-linux.
Applied to mainline.

Paul

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.
gcc/testsuite/
 * gfortran.dg/list_read.c: New test.

Index: config.h.in
===================================================================
RCS file: /var/cvsroot/gcc-cvs/gcc/libgfortran/config.h.in,v
retrieving revision 1.10
diff -u -p -r1.10 config.h.in
--- config.h.in 30 Oct 2004 16:23:22 -0000 1.10
+++ config.h.in 10 Nov 2004 00:49:29 -0000
@@ -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
 
Index: configure
===================================================================
RCS file: /var/cvsroot/gcc-cvs/gcc/libgfortran/configure,v
retrieving revision 1.18
diff -u -p -r1.18 configure
--- configure 5 Nov 2004 12:50:48 -0000 1.18
+++ configure 10 Nov 2004 01:44:10 -0000
@@ -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
Index: configure.ac
===================================================================
RCS file: /var/cvsroot/gcc-cvs/gcc/libgfortran/configure.ac,v
retrieving revision 1.13
diff -u -p -r1.13 configure.ac
--- configure.ac 5 Nov 2004 12:50:53 -0000 1.13
+++ configure.ac 10 Nov 2004 00:49:13 -0000
@@ -159,7 +159,7 @@ AC_CHECK_HEADER([complex.h],[AC_DEFINE([
 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])])
Index: io/read.c
===================================================================
RCS file: /var/cvsroot/gcc-cvs/gcc/libgfortran/io/read.c,v
retrieving revision 1.3
diff -u -p -r1.3 read.c
--- io/read.c 5 Aug 2004 08:37:42 -0000 1.3
+++ io/read.c 10 Nov 2004 01:52:25 -0000
@@ -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
- * rounding. */
+ * infinities.  */
 
 int
 convert_real (void *dest, const char *buffer, int length)
@@ -101,13 +101,18 @@ convert_real (void *dest, const char *bu
   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 *bu
   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
- * 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
- * the input. */
+   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,
+   and the implicit decimal point if a decimal point is not present in
+   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 lengt
       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 lengt
   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 lengt
  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)
-     *n++ = '0';
-   p++;
+   ndigits++;
+   *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
- * the d parameter before explict conversion takes place. */
+   of the exponent in order to take account of the scale factor and
+   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.
-   * The only thing that can go wrong at this point is overflow or
-   * underflow. */
+  if (exponent > 0)
+    {
+      edigits = 2;
+      i = exponent;
+    }
+  else
+    {
+      edigits = 3;
+      i = -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);
 
-  convert_precision_real (dest, val_sign, buffer, length, exponent);
+  /* Do the actual conversion.  */
+  string_to_real (dest, buffer, length);
 
   if (buffer != scratch)
      free_mem (buffer);

Attachment: read_float_1.f90
Description: Text document


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]