[Bug libfortran/48047] New: Incorrect output rounding of double precision numbers

[Bug libfortran/48047] New: Incorrect output rounding of double precision numbers

[thenlich at users dot sourceforge.net]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=48047

           Summary: Incorrect output rounding of double precision numbers
           Product: gcc
           Version: 4.6.0
            Status: UNCONFIRMED
          Severity: normal
          Priority: P3
         Component: libfortran
        AssignedTo: unassigned@gcc.gnu.org
        ReportedBy: thenlich@users.sourceforge.net


Created attachment 23603
  --> http://gcc.gnu.org/bugzilla/attachment.cgi?id=23603
Test case

The Fortran library does not round real(8) numbers correctly on output if 39
decimal digits are requested and real(16) is supported. This violates IEEE Std
754-2008 which demands that the minimum of supported significant digits for
correct rounding of all supported binary formats (H) is at least H=39 if the
binary128 format is supported.

Thus, GCC misses this requirement by 1 digit.

The attached program fails because the exact value, as given by
quadmath_snprintf(..., (__float128)0.14285714285714285) is
0.142857142857142849212692681248881854116916656494...

IEEE Std 754-2008 says:
===
5.12.2 External decimal character sequences representing finite numbers   
...
For the purposes of discussing the limits on correctly rounded conversion,
define the following quantities:
...
- for binary128, Pmin (binary128) = 36
...
- M = max(Pmin(bf)) for all supported binary formats bf
...

There might be an implementation-defined limit on the number of significant
digits that can be converted with correct rounding to and from supported binary
formats. That limit, H, shall be such that H >= M + 3 and it should be that H
is unbounded.

For all supported binary formats the conversion operations shall support
correctly rounded conversions to or from external character sequences for all
significant digit counts from 1 through H (that is, for all
expressible counts if H is unbounded).
...
===

[Bug libfortran/48047] Incorrect output rounding of double precision numbers

[kargl at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=48047

kargl at gcc dot gnu.org changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |kargl at gcc dot gnu.org

--- Comment #1 from kargl at gcc dot gnu.org 2011-03-09 18:01:42 UTC ---
gfortran does not claim conformance to IEEE 754, and
gfortran currently does not include even a crude 
implementation of the IEEE 754 intrinsic modules.
gfortran doesn't even guarantee 36 decimal digits in
its REAL(16) type

program z
print '(I0,1X,I0)', digits(1.e0_8),  digits(1.e0_16)
print '(I0,1X,I0)', precision(1.e0_8),  precision(1.e0_16)
end program z
troutmask:sgk[213] gfc4x -o z d.f90 && ./z
53 113
15 33

It is unclear what you want to happen here.

[Bug libfortran/48047] Incorrect output rounding of double precision numbers

[burnus at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=48047

Tobias Burnus <burnus at gcc dot gnu.org> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |burnus at gcc dot gnu.org

--- Comment #2 from Tobias Burnus <burnus at gcc dot gnu.org> 2011-03-09 19:06:12 UTC ---
For what it is worth: On x86-64-linux (openSUSE Factory, GLIBC 2.11.3), I get:

Got:      0.142857142857142849212692681248881854120
Expected: 0.142857142857142849212692681248881854116
Got:      0.142857142857142849212692681248881854120
Expected: 0.142857142857142849212692681248881854117

Using ifort 11.1 on the same machine:

Got:      0.142857142857142849212692681248881854100
Expected: 0.142857142857142849212692681248881854116
Got:      0.142857142857142849212692681248881854200
Expected: 0.142857142857142849212692681248881854117

Using crayftn 7.2.4, I get on another x86-64-linux (SLES 10/Cray):

Got:      0.142857142857142849000000000000000000000
Expected: 0.142857142857142849212692681248881854116
Got:      +.142857142857142849000000000000000000000
Expected: 0.142857142857142849212692681248881854117

My impression is that all twelve numbers are numerically identical and lead to
the same binary representation. Hence, I believe the three compilers print the
correct result.

[Bug libfortran/48047] Incorrect output rounding of double precision numbers

[jakub at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=48047

Jakub Jelinek <jakub at gcc dot gnu.org> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |jakub at gcc dot gnu.org

--- Comment #3 from Jakub Jelinek <jakub at gcc dot gnu.org> 2011-03-09 20:06:52 UTC ---
38 comes from:
1056      if (ndigits > MIN_FIELD_WIDTH - 4 - edigits)
1057        ndigits = MIN_FIELD_WIDTH - 4 - edigits;
(46 - 4 - 4).  snprintf/quadmath_snprintf returns only what it is asked for.

[Bug libfortran/48047] Incorrect output rounding of double precision numbers

[thenlich at users dot sourceforge.net]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=48047

--- Comment #4 from Thomas Henlich <thenlich at users dot sourceforge.net> 2011-03-10 09:36:52 UTC ---
Created attachment 23610
  --> http://gcc.gnu.org/bugzilla/attachment.cgi?id=23610
Patch to fix rounding issue

Proposing this patch (untested) for field width.

Trading 2 extra bytes for (partial) IEEE 754-2008 compliance seems like a good
thing.

It facilitates writing programs portable to other (IEEE 754-2008 compliant)
compilers.

[Bug libfortran/48047] Incorrect output rounding of double precision numbers

[thenlich at users dot sourceforge.net]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=48047

--- Comment #5 from Thomas Henlich <thenlich at users dot sourceforge.net> 2011-03-10 09:39:03 UTC ---
Created attachment 23611
  --> http://gcc.gnu.org/bugzilla/attachment.cgi?id=23611
Comprehensive test for IEEE 754-2008 clause 5.12.2 compliant output rounding

[Bug libfortran/48047] Incorrect output rounding of double precision numbers

[thenlich at users dot sourceforge.net]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=48047

Thomas Henlich <thenlich at users dot sourceforge.net> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
  Attachment #23611|0                           |1
        is obsolete|                            |

--- Comment #6 from Thomas Henlich <thenlich at users dot sourceforge.net> 2011-03-10 09:41:40 UTC ---
Created attachment 23612
  --> http://gcc.gnu.org/bugzilla/attachment.cgi?id=23612
Comprehensive test for IEEE 754-2008 clause 5.12.2 compliant output rounding
(corrected)

[Bug libfortran/48047] Incorrect output rounding of double precision numbers

[jvdelisle at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=48047

--- Comment #7 from Jerry DeLisle <jvdelisle at gcc dot gnu.org> 2011-03-11 00:59:04 UTC ---
I am OK with this simple patch, but will have to run the testsuite to make sure
its all OK on our other testsuite cases. (Not very many at the moment)

[Bug libfortran/48047] Incorrect output rounding of double precision numbers

[jvdelisle at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=48047

--- Comment #8 from Jerry DeLisle <jvdelisle at gcc dot gnu.org> 2011-03-11 18:49:30 UTC ---
Regression tested OK.  Can I get a release manager OK for this? If not, it can
wait for 4.7.

[Bug libfortran/48047] Incorrect output rounding of double precision numbers

[jakub at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=48047

--- Comment #9 from Jakub Jelinek <jakub at gcc dot gnu.org> 2011-03-11 19:12:36 UTC ---
(In reply to comment #8)
> Regression tested OK.  Can I get a release manager OK for this? If not, it can
> wait for 4.7.

This is fine for 4.6.

[Bug libfortran/48047] Incorrect output rounding of double precision numbers

[jvdelisle at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=48047

--- Comment #10 from Jerry DeLisle <jvdelisle at gcc dot gnu.org> 2011-03-12 00:49:55 UTC ---
Author: jvdelisle
Date: Sat Mar 12 00:49:51 2011
New Revision: 170895

URL: http://gcc.gnu.org/viewcvs?root=gcc&view=rev&rev=170895
Log:
2011-03-11  Jerry DeLisle  <jvdelisle@gcc.gnu.org>

    PR libgfortran/48047
    * io/write_float.def (write_float): Change MIN_FIELD_WIDTH to 48.

Modified:
    trunk/libgfortran/ChangeLog
    trunk/libgfortran/io/write_float.def

[Bug libfortran/48047] Incorrect output rounding of double precision numbers

[thenlich at users dot sourceforge.net]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=48047

Thomas Henlich <thenlich at users dot sourceforge.net> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
             Status|UNCONFIRMED                 |RESOLVED
         Resolution|                            |FIXED

--- Comment #11 from Thomas Henlich <thenlich at users dot sourceforge.net> 2011-05-02 11:28:08 UTC ---
It is now compliant with IEEE 754/2008.