[Bug libc/3479] Incorrect rounding in strtod()

["john at thesalmons dot org" <sourceware-bugzilla at sourceware dot org>]

------- Additional Comments From john at thesalmons dot org  2007-02-04 19:37 -------
I just encountered something very similar.  I started looking at C++
operator>>, which led me to scanf, which led me to strtof.  It appears
that glibc's strtof doesn't adhere to the "Recommended Practice" of
section 6.4.4.2 and 7.20.1.3 of C99.  I see numerous FAILs in the
attached code with strtof.  I haven't seen any problems with strtod,
but I haven't done exhaustive tests.  Platforms:

  Opteron, Centos, glibc-2.3.4-2, gcc-4.1.1
  Pentium 4, Fedora, glibc-2.3.3-27.1, gcc-3.3.3-7

I have not tested against glibc's CVS, but even if it passes, I would
suggest that something like the tests here belong in the test suite.

The recommended practice for literal floating constants in 6.4.4.2 is:

     The translation-time conversion of floating constants should
     match the execution-time conversion of character strings by library
     functions, such as strtod, given matching inputs suitable for both
     conversions, the same result format, and default execution-time
     rounding.63) 

     63) The specification for the library functions
     recommends more accurate conversion than required for floating
     constants (see 7.20.1.3).

The attached code checks this condition.  It's *much* easier to check
this than to check the more detailed requirement of 7.20.1.3.  As a
result, the code doesn't actually prove that glibc is in error.  It's
possible that gcc's handling of literal floats is in error.  


The recommended practice for strto* in 7.20.1.3 is:

     If the subject sequence has the decimal form and at most
     DECIMAL_DIG (defined in <float.h>) significant digits, the result
     should be correctly rounded. If the subject sequence D has the
     decimal form and more than DECIMAL_DIG significant digits,
     consider the two bounding, adjacent decimal strings L and U, both
     having DECIMAL_DIG significant digits, such that the values of L,
     D, and U satisfy L <= D <= U. The result should be one of the
     (equal or adjacent) values that would be obtained by correctly
     rounding L and U according to the current rounding direction,
     with the extra stipulation that the error with respect to D
     should have a correct sign for the current rounding direction.

The decimal strings in the attached test case are all less than or
equal to DECIMAL_DIG in length.

A close look at the output suggests that glibc is in error, at least
if one believes the output of printf("%.30g").  For example, for the
first FAIL:

The input string is:
"1.00000005960464477550"
The possible neighboring floats are:
"1.00000011920928955078125000000"  (printf("%.30g", 1 + FLT_EPSILON))
and
"1.0"

According to dc, the distance from the input to 1+FLT_EPSILON is:
.00000005960464477528125000000
and the distance from the input to 1.0 is
.00000005960464477550

so 1+FLT_EPSILON should be the "correct" properly rounded result.
strtof returns 1.0.

---------------------------------------

[jsalmon@river lexical_cast]$ cat badstrtod2.c
#define _ISOC99_SOURCE_
#include <stdio.h>
#include <stdlib.h>
#include <limits.h>
#include <math.h>
#include <float.h>

void checkf(const char *s, float f){
    float cf = strtof(s, NULL);
    long double cld = strtold(s, NULL);

    if( cf != f ){
	printf("FAIL strtof (\"%s\") %a is not the same as literal %a\n",
	       s, cf, f);
    }else{
	printf("OK   strtof (\"%s\") %a is the same as literal\n",
	       s, cf);
    }
}

void checkd(const char *s, double d){
    double cd = strtod(s, NULL);
    if( cd != d ){
	printf("FAIL strtod (\"%s\") %a is not the same as literal %a\n",
	       s, cd, d);
    }else{
	printf("OK   strtod (\"%s\") %a is the same as literal\n",
	       s, cd);
    }
}

int main(int argc, char **argv){
    long double ld;

    /* Put some context in the output.  These are not used for testing */
    ld = 1.L + (long double)FLT_EPSILON;
    printf("1 + FLT_EPSILON: %La %#.30Lg\n", ld, ld);
    ld = 1.L + ((long double)FLT_EPSILON)/2.;
    printf("1 + FLT_EPSILON/2.: %La %#.30Lg\n", ld, ld);
    ld = 1.L + ((long double)FLT_EPSILON)/2. + LDBL_EPSILON;
    printf("1 + FLT_EPSILON/2. + LDBL_EPSILON: %La %#.30Lg\n", ld, ld);

#define CHKF(s) checkf(#s, s##f);
    /* 1 + FLT_EPSILON/2. + LDBL_EPSILON */
    CHKF(1.00000005960464477550);
    CHKF(1.0000000596046447755);
    CHKF(1.000000059604644776);
    CHKF(1.000000059604644775);
    CHKF(1.00000005960464478);
    CHKF(1.0000000596046448);
    CHKF(1.000000059604645);
    CHKF(1.00000005960464);
    CHKF(1.0000000596046);
    CHKF(1.000000059605);
    CHKF(1.00000005960);
    CHKF(1.0000000596);
    CHKF(1.000000060);
    CHKF(1.00000006);
    CHKF(1.0000001);
    CHKF(1.000000);

    ld = 1.L + (long double)DBL_EPSILON;
    printf("1 + DBL_EPSILON: %La %#.30Lg\n", ld, ld);
    ld = 1.L + ((long double)DBL_EPSILON)/2.;
    printf("1 + DBL_EPSILON/2.: %La %#.30Lg\n", ld, ld);
    ld = 1.L + ((long double)DBL_EPSILON)/2. + LDBL_EPSILON;
    printf("1 + DBL_EPSILON/2. + LDBL_EPSILON: %La %#.30Lg\n", ld, ld);

#define CHKD(s) checkd(#s, s);
    CHKD(1.00000000000000011113);
    CHKD(1.00000000000000011103);
    CHKD(1.00000000000000011102);
    CHKD(1.00000000000000011101);
    CHKD(1.0000000000000001111);
    CHKD(1.000000000000000111);
    CHKD(1.00000000000000011);
    CHKD(1.0000000000000001);

    return 0;
}
[jsalmon@river lexical_cast]$ 
[jsalmon@river lexical_cast]$ gcc -g -std=c99  badstrtod2.c -lm -o badstrtod2 &&
badstrtod2 
1 + FLT_EPSILON: 0x8.00001p-3 1.00000011920928955078125000000
1 + FLT_EPSILON/2.: 0x8.000008p-3 1.00000005960464477539062500000
1 + FLT_EPSILON/2. + LDBL_EPSILON: 0x8.000008000000001p-3
1.00000005960464477549904521725
FAIL strtof ("1.00000005960464477550") 0x1p+0 is not the same as literal
0x1.000002p+0
FAIL strtof ("1.0000000596046447755") 0x1p+0 is not the same as literal
0x1.000002p+0
FAIL strtof ("1.000000059604644776") 0x1p+0 is not the same as literal 0x1.000002p+0
OK   strtof ("1.000000059604644775") 0x1p+0 is the same as literal
FAIL strtof ("1.00000005960464478") 0x1p+0 is not the same as literal 0x1.000002p+0
FAIL strtof ("1.0000000596046448") 0x1p+0 is not the same as literal 0x1.000002p+0
FAIL strtof ("1.000000059604645") 0x1p+0 is not the same as literal 0x1.000002p+0
OK   strtof ("1.00000005960464") 0x1p+0 is the same as literal
OK   strtof ("1.0000000596046") 0x1p+0 is the same as literal
FAIL strtof ("1.000000059605") 0x1p+0 is not the same as literal 0x1.000002p+0
OK   strtof ("1.00000005960") 0x1p+0 is the same as literal
OK   strtof ("1.0000000596") 0x1p+0 is the same as literal
OK   strtof ("1.000000060") 0x1.000002p+0 is the same as literal
OK   strtof ("1.00000006") 0x1.000002p+0 is the same as literal
OK   strtof ("1.0000001") 0x1.000002p+0 is the same as literal
OK   strtof ("1.000000") 0x1p+0 is the same as literal
1 + DBL_EPSILON: 0x8.0000000000008p-3 1.00000000000000022204460492503
1 + DBL_EPSILON/2.: 0x8.0000000000004p-3 1.00000000000000011102230246252
1 + DBL_EPSILON/2. + LDBL_EPSILON: 0x8.000000000000401p-3
1.00000000000000011113072267976
OK   strtod ("1.00000000000000011113") 0x1.0000000000001p+0 is the same as literal
OK   strtod ("1.00000000000000011103") 0x1.0000000000001p+0 is the same as literal
OK   strtod ("1.00000000000000011102") 0x1p+0 is the same as literal
OK   strtod ("1.00000000000000011101") 0x1p+0 is the same as literal
OK   strtod ("1.0000000000000001111") 0x1.0000000000001p+0 is the same as literal
OK   strtod ("1.000000000000000111") 0x1p+0 is the same as literal
OK   strtod ("1.00000000000000011") 0x1p+0 is the same as literal
OK   strtod ("1.0000000000000001") 0x1p+0 is the same as literal
[jsalmon@river lexical_cast]$ 


-- 


http://sourceware.org/bugzilla/show_bug.cgi?id=3479

------- You are receiving this mail because: -------
You are on the CC list for the bug, or are watching someone who is.