[Bug c/51441] New: Incorrect FP rounding on addition of doubles

[Bug c/51441] New: Incorrect FP rounding on addition of doubles

[exponent-bias at yandex dot ru]

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

             Bug #: 51441
           Summary: Incorrect FP rounding on addition of doubles
    Classification: Unclassified
           Product: gcc
           Version: 4.4.3
            Status: UNCONFIRMED
          Severity: normal
          Priority: P3
         Component: c
        AssignedTo: unassigned@gcc.gnu.org
        ReportedBy: exponent-bias@yandex.ru


A well-known sentence of IEEE 754 reads: "Each of the computational operations
that return a numeric result specified by this standard shall be performed as
if it first produced an intermediate result correct to infinite precision and
with unbounded range, and then rounded that intermediate result, if necessary,
to fit in the destination's format..."

Given two operands in IEEE's binary32 format:
    0x1.002p0  (C99 hexadecimal notation)
and
    1.6331239353195372e+16
I believe the result of addition of the operands under the to-nearest rounding
mode shall be:
    1.6331239353195374e+16
whereas GCC yields:
    1.6331239353195372e+16
(note the difference in the last sign of the mantissa).

The analysis is this: the difference of exponents of the operands is 53 bits
which is also the width of mantissa in the binary32 format, including the
implied most significant raised bit. Thus, mantissa of the intermediate result
of the addition should be result of concatenation of mantissas of the operands:
mantissa of 1.6331239353195372e+16 followed by mantissa of 0x1.002p0. That is,
the most significant mantissa's bit of 0x1.002p0 is the first bit that shall be
rounded away to fit the destination's format. We know that first bit is raised
and we know there is one more bit raised in the mantissa, so under the
to-nearest rounding mode the intermediate result shall be rounded to the next
representable number that is not less than the intermediate result, which is
1.6331239353195374e+16.

The test code:
----
#include <stdlib.h>
#include <stdio.h>

double a, b, r;

int main(void)
{
    a = 0x1.002p0;
    b = 1.6331239353195372e+16;

    r = a + b;

    printf("%+.16e + %+.16e = %+.16e\n", a, b, r);

    return EXIT_SUCCESS;
}
----

Compilation flags:
gcc --float-store -g -O0 -o test test.c

Actual output:
+1.0004882812500000e+00 + +1.6331239353195372e+16 = +1.6331239353195372e+16

Expected output:
+1.0004882812500000e+00 + +1.6331239353195372e+16 = +1.6331239353195374e+16

Compiler version:
gcc-4.4.real (Ubuntu 4.4.3-4ubuntu5) 4.4.3

[Bug c/51441] Incorrect FP rounding on addition of doubles

[exponent-bias at yandex dot ru]

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

--- Comment #1 from Alex Makov <exponent-bias at yandex dot ru> 2011-12-06 20:05:39 UTC ---
The references to binary32 should be read as binary64, of course--just a typo.

[Bug target/51441] Incorrect FP rounding on addition of doubles

[pinskia at gcc dot gnu.org]

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

Andrew Pinski <pinskia at gcc dot gnu.org> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
             Target|                            |i?86-*-*
          Component|c                           |target

--- Comment #2 from Andrew Pinski <pinskia at gcc dot gnu.org> 2011-12-06 20:15:55 UTC ---
This only effects x87.

[Bug target/51441] Incorrect FP rounding on addition of doubles

[joseph at codesourcery dot com]

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

--- Comment #3 from joseph at codesourcery dot com <joseph at codesourcery dot com> 2011-12-06 20:17:36 UTC ---
Please identify the target for which the compiler is configured.  If 
32-bit x86, then note that in accordance with FLT_EVAL_METHOD you have 
double rounding (two double values are converted to long double, added as 
long doubles and then the long double sum rounded to double).

[Bug c/51441] Incorrect FP rounding on addition of doubles

[exponent-bias at yandex dot ru]

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

Alex Makov <exponent-bias at yandex dot ru> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
          Component|target                      |c

--- Comment #4 from Alex Makov <exponent-bias at yandex dot ru> 2011-12-07 07:20:29 UTC ---
(In reply to comment #3).

> Please identify the target for which the compiler is configured.  If 
> 32-bit x86, then note that in accordance with FLT_EVAL_METHOD you have 
> double rounding (two double values are converted to long double, added as 
> long doubles and then the long double sum rounded to double).

The target is 32-bit x86. Given the "intermediate result" is first rounded to
long double and then rounded once more to double, then the result should match
with what GCC yields indeed. Moreover, after adding the "-mpc64" option to the
command line the test code prints +1.6331239353195374e+16 which is the expect
output. So, this whole thing looks to be false alarm.

Anyway, thanks a lot for your feedback.

[Bug c/51441] Incorrect FP rounding on addition of doubles

[jakub at gcc dot gnu.org]

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

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

           What    |Removed                     |Added
----------------------------------------------------------------------------
             Status|UNCONFIRMED                 |RESOLVED
                 CC|                            |jakub at gcc dot gnu.org
         Resolution|                            |INVALID

--- Comment #5 from Jakub Jelinek <jakub at gcc dot gnu.org> 2011-12-07 09:35:16 UTC ---
Invalid.