[Bug c/19837] New: Floating-point numerical regression with and without optimization

["wirawan0 at softhome dot net" <gcc-bugzilla at gcc dot gnu dot org>]

VERSION KNOWN TO BE AFFECTED: 3.4.3, 3.2.2

My gcc gives the following version/compilation info:

~/Projects/BUGS $ g++ -v
Reading specs from /usr/local/gcc-3.4.3/lib/gcc/i686-pc-linux-gnu/3.4.3/specs
Configured with: ../gcc-3.4.3/configure --enable-threads=posix
--prefix=/usr/local/gcc-3.4.3 --enable-languages=c,c++,f77 --with-system-zlib
--enable-shared --enable-__cxa_atexit
Thread model: posix
gcc version 3.4.3


DESCRIPTION

OK. Here is the origin of this bug: I have a formula:

alpha = 1.0 + (x * x) - x * sqrt(x * x + 2.0);

where alpha and x are double-precision numbers. Now, given a very large
value of x (for example, x >= 5000.0), then g++ (and also gcc) would
evaluate wrongly the formula above. Somehow, the outcome of the value
depends on the optimization level.

Let's take x = 20000.0 as a severe testcase. If you do Taylor expansion on
alpha for large x values (thus, expanding in terms of 1/x^2), the leading
term of alpha is

0.5 / (x * x)

So, the leading term is 1.25e-09, which is essentially a an accurate
answer. Precise numerical evaluation using Waterloo Maple (version 9, in
this case) yields

1.24999999687500000976562496582032e-09   ("exact result", 33-digits precision)

What happens if I evaluate the function above using gcc?

gcc command                result (alpha)
g++ testcase.cpp -O0       2.526212483644485e-08
g++ testcase.cpp -O1       1.251464709639549e-09
g++ testcase.cpp -O2       1.251464709639549e-09
g++ testcase.cpp -O3       0

Use the snippet below as the testcase.

What is this? The -O0 compiled code is exceedingly wrong by predicting a
number 20 times bigger.
The -O1 and -O2 results are essentially correct. Why they disagree with
the exact answer is understandable, since the subtraction is between two
VERY LARGE numbers:

alpha = 400000001 - 20000*sqrt(400000002)

If the sqrt thing can be accurate, then we should be left with a small
number of residual precision, such as predicted by gcc-compiled code above
with -O1 and -O2.

NOTE that the same formula is used in another program (my own code, not
included here). There, the -O3 executable code gives the same number as
-O1 and -O2 above. So, there is no consistency here, although the gist of
the problem remains.

I'm sorry if this bug is probably a result of pushing numerics to its
edge. I should evaluate the expression using a series expansion instead of
naively using the formula like above. But I report this problem anyway, at
least to document the problem to the community.

Although this bug is (originally) rather specific to a formula, I think
that other like this one, which computes the difference between two large
numbers, could suffer from a similar problem.


REPRODUCING THE ERROR

Use the stock g++ 3.4.3. Run on an x86 machine. I'm not sure if AMD chips
would yield the same error. I bet it would, because of the executable
binary compatibility. But you have to test it.
This may have been a specific problem to the Intel x86 machines only (I
have tested this problem using Pentium M and Pentium 4 chips, see below).
I have NOT tested using AMD chips. This may not apply to other machines.

I have tried the `-mieee-fp' switch, it does not fix the problem.

This error is also reproducible if you compile the testcase using
`gcc -x c' (the C compiler) instead of g++ (the C++ compiler).


OTHER VERSIONS

g++ version 3.2.2 running on an Intel Pentium 4 also reproduces the same
results above, except that the -O3 binary gives the same answer as -O2
above, instead of zero.
The g++'s identification string is:

[wirawan@snowflake wirawan]$ g++ -v
Reading specs from /usr/lib/gcc-lib/i386-redhat-linux/3.2.2/specs
Configured with: ../configure --prefix=/usr --mandir=/usr/share/man
--infodir=/usr/share/info --enable-shared --enable-threads=posix
--disable-checking --with-system-zlib --enable-__cxa_atexit
--host=i386-redhat-linux
Thread model: posix
gcc version 3.2.2 20030222 (Red Hat Linux 3.2.2-5)

g++ version 3.1 on Alpha ev67 gives zero regardless of the optimization
level (-O0, -O1, -O2, -O3).

g++ version 3.3 on an UltraSPARC III box also gives zero regardless of the
optimization level (-O0, -O1, -O2, -O3).

SUN C++ (version "5.4 Patch 111715-08 2003/06/08") also gives zero
regardless of the optimization level (no switch, or -xO3).

Intel C++ version 7.1 (Build 20030617Z) gives 1.251464709639549e-09
regardless of the optimization level.

-- 
           Summary: Floating-point numerical regression with and without
                    optimization
           Product: gcc
           Version: 3.4.3
            Status: UNCONFIRMED
          Severity: normal
          Priority: P2
         Component: c
        AssignedTo: unassigned at gcc dot gnu dot org
        ReportedBy: wirawan0 at softhome dot net
                CC: gcc-bugs at gcc dot gnu dot org
 GCC build triplet: i686-pc-linux-gnu
  GCC host triplet: i686-pc-linux-gnu
GCC target triplet: i686-pc-linux-gnu


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