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
