Floating point problems (accuracy).

[James E Wilson]

Phil Prentice wrote:
>   We are having a few problems regarding floating point inaccuracies under
>   Linux.  Please see the code below.

This is PR 323 in our bug database.  One of our most commonly reported 
problems.  We call this the excess precision problem.
	http://gcc.gnu.org/bugzilla/show_bug.cgi?id=323

This problem has been known for about 15 years now, and hasn't been 
fixed yet.  My suspicion is that it never will be.

The problem here is partly due to the design of the x87 FP unit, which 
supports full precision 80-bit operations, but not 32-bit or 64-bit 
operations.  The problem is also partly due to the design of the FSF gcc 
i386 backend, which lies, and says we do have 32-bit and 64-bit 
operations.  This results in a mismatch between what the compiler does 
and what the hardware does, which results in various FP accuracy problems.

Unfortunately, there are no good fixes for this.  If we fix the gcc i386 
backend, then we will likely have a significant FP performance loss, 
which may annoy more people than the excess precision problem does.  And 
we will still be left with double rounding problems, which means we 
still won't have completely correct FP problems.  The double rounding 
problem can't be fixed in gcc, it requires a hardware solution.  A new 
x87 FP optimization pass may help performance, but would require a lot 
of work, and it isn't clear that all of the work is worth the benefit.

The only fool proof solution is to avoid use of the x87 FP unit.  The 
good news is that if you have an AMD AMD64 or Intel EM64T system, and 
run x86_64-linux on it, then you will get correct FP code by default. 
So we do have a solution, of sorts, though it requires the whole world 
to migrate from x86-linux to x86_64-linux.  x86_64-linux does not have a 
problem because it uses the SSE registers instead of the x87 registers 
by default for 64-bit code.

Meanwhile, there are a few workarounds
1) You can try changing the hardware FP rounding mode to double 
precision.  You lose long double support, but you now get correct double 
precision FP.  You still have excess precision problems for float.  This 
may break libraries (e.g. glibc math library) that depend on long double 
support.
2) Use -ffloat-store.  This forces variables to be allocated to memory 
instead of FP registers, so that they will be corrected rounded.  This 
will eliminate the excess precision problem for some programs, but not 
all of them, as you may still have excess precision problems for 
temporaries.  This will reduce FP performance.
3) Use -fpmath=sse.  Only works for recent compilers, and recent processors.

If you are ambitious, you could try adding an option to the x86 backend 
that disables all of the 32-bit and 64-bit x87 FP instruction patterns. 
  This should result in slow but more correct FP code which may be 
useful to some people.  This is a probably a significant amount of work 
though, and it isn't clear how useful the result will be.
-- 
Jim Wilson, GNU Tools Support, http://www.SpecifixInc.com