[Patch, Fortran] Rework of floating point number reading

[Daniel Kraft <d at domob dot eu>]

Hi,

here's a full patch for my proposed changes to libgfortran's floating 
point reading.  With the attached test case, I see a speed-up from 31s 
to 18s (the 31 being with unpatched gfortran 4.3, but also 4.4 had a 
similar speed before I did the patch).  In libgfortran/io/read.c:read_f, 
currently the number is "parsed" but only to put it back as string, 
handling special things like BLANK_NULL/BLANK_ZERO; this string is then 
again parsed using strtod & friends.  With the patch, strtod is never 
called but instead the number directly calculated during this first pass.

I see no regressions on GNU/Linux-x86-32.  Is this patch ok in general? 
 I do however not want to commit it as-is; my goal is to get rid of 
io/read.c:convert_real, but I saw it is used two times in list_read.c, 
also.  There seems to be a similar practice of saving to a string and 
re-parsing with convert_real; if my patch is ok, I would like to work 
out a common system for handling all of those cases in one and do only 
one pass there, too.

As I'm no expert on libgfortran and IO, I would very welcome any 
comments!  Also, as I fear my patch could introduce some rounding-errors 
to parsing, I would like to get test results from some other platforms 
as well.

And finally, profiling showed that the test-program attached spends most 
time in real_f after the patch, but still around 6% in pow; I guess this 
is because of the "powl (10.0, exponent)" at the end; by using a 
lookup-table, we could get rid of this call.  Do you think this is worth it?

Thanks,
Daniel

--
Done:  Arc-Bar-Cav-Rog-Sam-Tou-Val-Wiz
To go: Hea-Kni-Mon-Pri-Ran

Index: libgfortran/io/read.c
===================================================================
--- libgfortran/io/read.c	(revision 142940)
+++ libgfortran/io/read.c	(working copy)
@@ -786,15 +786,23 @@ read_f (st_parameter_dt *dtp, const fnod
   size_t wu;
   int w, seen_dp, exponent;
   int exponent_sign, val_sign;
-  int ndigits;
-  int edigits;
-  int i;
-  char *p, *buffer;
-  char *digits;
-  char scratch[SCRATCH_SIZE];
+  char *p;
+  int dec_exp; /* Additional exponent from placement of decimal point.  */
+
+  /* Take maximum precision result variable.  */
+#ifdef HAVE_GFC_REAL_16
+  GFC_REAL_16 result;
+#elif defined (HAVE_GFC_REAL_10)
+  GFC_REAL_10 result;
+#else
+  GFC_REAL_8 result;
+#endif
 
   val_sign = 1;
   seen_dp = 0;
+  result = 0.0;
+  dec_exp = 0;
+  exponent = 0;
   wu = f->u.w;
 
   p = gfc_alloca (wu);
@@ -806,7 +814,10 @@ read_f (st_parameter_dt *dtp, const fnod
 
   p = eat_leading_spaces (&w, p);
   if (w == 0)
-    goto zero;
+    {
+      result = 0.0;
+      goto done;
+    }
 
   /* Optional sign */
 
@@ -821,7 +832,10 @@ read_f (st_parameter_dt *dtp, const fnod
   exponent_sign = 1;
   p = eat_leading_spaces (&w, p);
   if (w == 0)
-    goto zero;
+    {
+      result = 0.0;
+      goto done;
+    }
 
   /* A digit, a '.' or a exponent character ('e', 'E', 'd' or 'D')
      is required at this point */
@@ -830,11 +844,8 @@ read_f (st_parameter_dt *dtp, const fnod
       && *p != 'e' && *p != 'E')
     goto bad_float;
 
-  /* Remember the position of the first digit.  */
-  digits = p;
-  ndigits = 0;
-
-  /* Scan through the string to find the exponent.  */
+  /* Scan through the string to find the exponent.  Calculate value of
+     mantissa on the way.  */
   while (w > 0)
     {
       switch (*p)
@@ -850,8 +861,22 @@ read_f (st_parameter_dt *dtp, const fnod
 	  if (seen_dp)
 	    goto bad_float;
 	  seen_dp = 1;
-	  /* Fall through */
+	  p++;
+	  w--;
+          break;
 
+        case ' ':
+          if (dtp->u.p.blank_status == BLANK_ZERO)
+            *p = '0';
+          else if (dtp->u.p.blank_status == BLANK_NULL)
+            {
+              ++p;
+              --w;
+              break;
+            }
+          else
+            goto done;
+          /* Fall through.  */
 	case '0':
 	case '1':
 	case '2':
@@ -862,10 +887,12 @@ read_f (st_parameter_dt *dtp, const fnod
 	case '7':
 	case '8':
 	case '9':
-	case ' ':
-	  ndigits++;
-	  p++;
-	  w--;
+          result *= 10;
+          result += (*p - '0');
+          if (seen_dp)
+            --dec_exp;
+          ++p;
+          --w;
 	  break;
 
 	case '-':
@@ -900,35 +927,6 @@ read_f (st_parameter_dt *dtp, const fnod
   next_record (dtp, 1);
   return;
 
-  /* The value read is zero */
- zero:
-  switch (length)
-    {
-      case 4:
-	*((GFC_REAL_4 *) dest) = 0;
-	break;
-
-      case 8:
-	*((GFC_REAL_8 *) dest) = 0;
-	break;
-
-#ifdef HAVE_GFC_REAL_10
-      case 10:
-	*((GFC_REAL_10 *) dest) = 0;
-	break;
-#endif
-
-#ifdef HAVE_GFC_REAL_16
-      case 16:
-	*((GFC_REAL_16 *) dest) = 0;
-	break;
-#endif
-
-      default:
-	internal_error (&dtp->common, "Unsupported real kind during IO");
-    }
-  return;
-
   /* At this point the start of an exponent has been found */
  exp1:
   while (w > 0 && *p == ' ')
@@ -972,7 +970,7 @@ read_f (st_parameter_dt *dtp, const fnod
       while (w > 0)
         {
           if (*p != ' ')
-  	  goto bad_float;
+            goto bad_float;
           p++;
           w--;
         }
@@ -1008,61 +1006,32 @@ read_f (st_parameter_dt *dtp, const fnod
   if (!seen_dp)
     exponent -= f->u.real.d;
 
-  if (exponent > 0)
-    {
-      edigits = 2;
-      i = exponent;
-    }
-  else
-    {
-      edigits = 3;
-      i = -exponent;
-    }
-
-  while (i >= 10)
+  result *= val_sign * powl (10.0, exponent + dec_exp);
+  switch (length)
     {
-      i /= 10;
-      edigits++;
-    }
+    case 4:
+      *((GFC_REAL_4*) dest) = (GFC_REAL_4) result;
+      break;
 
-  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 decimal point in place.  */
-  p = buffer;
-  if (val_sign < 0)
-    *(p++) = '-';
-  for (; ndigits > 0; ndigits--)
-    {
-      if (*digits == ' ')
-        {
-	  if (dtp->u.p.blank_status == BLANK_ZERO) *digits = '0';
-	  if (dtp->u.p.blank_status == BLANK_NULL)
-            {
-              digits++;
-              continue;
-            } 
-        }
-      *p = *digits;
-      p++;
-      digits++;
-    }
-  *(p++) = 'e';
-  sprintf (p, "%d", exponent);
+    case 8:
+      *((GFC_REAL_8*) dest) = (GFC_REAL_8) result;
+      break;
 
-  /* Do the actual conversion.  */
-  convert_real (dtp, dest, buffer, length);
+#if defined(HAVE_GFC_REAL_10)
+    case 10:
+      *((GFC_REAL_10*) dest) = (GFC_REAL_10) result;
+      break;
+#endif
 
-  if (buffer != scratch)
-     free_mem (buffer);
+#if defined(HAVE_GFC_REAL_16)
+    case 16:
+      *((GFC_REAL_16*) dest) = (GFC_REAL_16) result;
+      break;
+#endif
 
+    default:
+      internal_error (&dtp->common, "Unsupported real kind during IO");
+    }
 }
 
 

! Time floating-point IO.

PROGRAM floatIO
  IMPLICIT NONE

  ! Floating-point numbers to read/write
  INTEGER, PARAMETER :: wp = 8
  REAL(KIND=wp), PARAMETER :: nums(6) = (/ 1.0_wp, .023_wp, -100.005_wp, &
                                           -.5e7_wp, 1e10_wp, &
                                           123456789000.00001_wp /)

  ! Number of repetitions.
  INTEGER, PARAMETER :: rep = 4000000

  CHARACTER(len=255) :: buffer
  REAL(KIND=8) :: readNums(SIZE (nums))
  INTEGER :: i

100 FORMAT (6F24.6)
  WRITE (buffer, 100) nums
  PRINT *, buffer

  DO i = 1, rep
    READ (buffer, 100) readNums
    IF (i == 1 .AND. ANY (readNums /= nums)) THEN
      PRINT *, nums
      PRINT *, readNums
      CALL abort ()
    END IF
  END DO
END PROGRAM floatIO