From 5e7b82bb9e3d9f11c49ed360c650aee50b472e97 Mon Sep 17 00:00:00 2001 From: pbrook Date: Wed, 10 Nov 2004 02:19:27 +0000 Subject: [PATCH] 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. git-svn-id: svn+ssh://gcc.gnu.org/svn/gcc/trunk@90382 138bc75d-0d04-0410-961f-82ee72b054a4 --- gcc/testsuite/ChangeLog | 5 + gcc/testsuite/gfortran.dg/read_float_1.f90 | 12 ++ libgfortran/ChangeLog | 10 ++ libgfortran/config.h.in | 3 + libgfortran/configure | 3 +- libgfortran/configure.ac | 2 +- libgfortran/io/read.c | 239 ++++++++++------------------- 7 files changed, 118 insertions(+), 156 deletions(-) create mode 100644 gcc/testsuite/gfortran.dg/read_float_1.f90 diff --git a/gcc/testsuite/ChangeLog b/gcc/testsuite/ChangeLog index df1eaa3b8ed..ce5a2b13aef 100644 --- a/gcc/testsuite/ChangeLog +++ b/gcc/testsuite/ChangeLog @@ -1,3 +1,8 @@ +2004-11-10 Paul Brook + + PR fortran/18218 + * gfortran.dg/list_read.c: New test. + 2004-11-09 Joseph S. Myers PR c/18322 diff --git a/gcc/testsuite/gfortran.dg/read_float_1.f90 b/gcc/testsuite/gfortran.dg/read_float_1.f90 new file mode 100644 index 00000000000..e38036f8b92 --- /dev/null +++ b/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 + diff --git a/libgfortran/ChangeLog b/libgfortran/ChangeLog index 3028e07d98b..000d49656d0 100644 --- a/libgfortran/ChangeLog +++ b/libgfortran/ChangeLog @@ -1,3 +1,13 @@ +2004-11-10 Paul Brook + + 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 * configure.ac: Use AC_PROG_FC, FC and FCFLAGS instead of diff --git a/libgfortran/config.h.in b/libgfortran/config.h.in index ecb9a6c896f..d31c21b84c9 100644 --- a/libgfortran/config.h.in +++ b/libgfortran/config.h.in @@ -156,6 +156,9 @@ /* Define to 1 if you have the 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 header file. */ #undef HAVE_SYS_MMAN_H diff --git a/libgfortran/configure b/libgfortran/configure index 3cf4dbee2be..e486b4e9ec8 100755 --- a/libgfortran/configure +++ b/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 diff --git a/libgfortran/configure.ac b/libgfortran/configure.ac index d99eded4428..3adcfd61c3d 100644 --- a/libgfortran/configure.ac +++ b/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])]) diff --git a/libgfortran/io/read.c b/libgfortran/io/read.c index 260a3dca5c6..6999158c13a 100644 --- a/libgfortran/io/read.c +++ b/libgfortran/io/read.c @@ -24,6 +24,7 @@ Boston, MA 02111-1307, USA. */ #include #include #include +#include #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 *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 - * 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 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) - *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); -- 2.11.4.GIT