This is GCC Bugzilla
This is GCC Bugzilla Version 2.20+
View Bug Activity | Format For Printing | Clone This Bug
A simple program gives incorrect floating point results when compiled with the -O flag. The complete source code is in the How-To-Repeat section below. This program prints "error" when compiled with -O; clearly wrong. The unoptimized code correctly prints no output. Also, if I try to print out the value of y2 before the comparison with y, then the correct value is output, and there is no "error" message, even when the code is optimized. In other words, interposing a print statement between the declaration (and initialization) of y2 and the comparison changes the program behavior. The error disappears for certain different values of the numerical constants. The error also goes away when an "l" suffix is appended to the two 1.0 constants. However, it is also possible to create the same type of error using only variables and no numerical constants in the assignments to y and y2; in this situation there is no "l" suffix workaround. The gzipped .ii file is attached. Here is the output of gcc: > /usr/local/bin/gcc -v --save-temps -O bug.C Reading specs from /usr/local/lib/gcc-lib/i686-pc-linux-gnu/2.95.2/specs gcc version 2.95.2 19991024 (release) /usr/local/lib/gcc-lib/i686-pc-linux-gnu/2.95.2/cpp -lang-c++ -v -D__GNUC__=2 -D__GNUG__=2 -D__GNUC_MINOR__=95 -D__cplusplus -D__ELF__ -Dunix -D__i386__ -Dlinux -D__ELF__ -D__unix__ -D__i386__ -D__linux__ -D__unix -D__linux -Asystem(posix) -D__EXCEPTIONS -D__OPTIMIZE__ -Acpu(i386) -Amachine(i386) -Di386 -D__i386 -D__i386__ -Di686 -Dpentiumpro -D__i686 -D__i686__ -D__pentiumpro -D__pentiumpro__ bug.C bug.ii GNU CPP version 2.95.2 19991024 (release) (i386 Linux/ELF) #include "..." search starts here: #include <...> search starts here: /usr/local/lib/gcc-lib/i686-pc-linux-gnu/2.95.2/../../../../include/g++-3 /usr/local/include /usr/local/lib/gcc-lib/i686-pc-linux-gnu/2.95.2/../../../../i686-pc-linux-gnu/include /usr/local/lib/gcc-lib/i686-pc-linux-gnu/2.95.2/include /usr/include End of search list. The following default directories have been omitted from the search path: End of omitted list. /usr/local/lib/gcc-lib/i686-pc-linux-gnu/2.95.2/cc1plus bug.ii -quiet -dumpbase bug.cc -O -version -o bug.s GNU C++ version 2.95.2 19991024 (release) (i686-pc-linux-gnu) compiled by GNU C version 2.95.2 19991024 (release). as -V -Qy -o bug.o bug.s GNU assembler version 2.9.1 (i386-redhat-linux), using BFD version 2.9.1.0.24 /usr/local/lib/gcc-lib/i686-pc-linux-gnu/2.95.2/collect2 -m elf_i386 -dynamic-linker /lib/ld-linux.so.2 /usr/lib/crt1.o /usr/lib/crti.o /usr/local/lib/gcc-lib/i686-pc-linux-gnu/2.95.2/crtbegin.o -L/usr/local/lib/gcc-lib/i686-pc-linux-gnu/2.95.2 -L/usr/local/lib bug.o -lgcc -lc -lgcc /usr/local/lib/gcc-lib/i686-pc-linux-gnu/2.95.2/crtend.o /usr/lib/crtn.o Release: gcc version 2.95.2 19991024 (release) Environment: (Red Hat) Linux 2.2.12-20, i686 pc How-To-Repeat: #include <stdio.h> void test(double x, double y) { const double y2 = x + 1.0; if (y != y2) printf("error\n"); } void main() { const double x = .012; const double y = x + 1.0; test(x, y); }
Fix: This bug is elusive and only occurs for certain floating point values. It goes away by tweaking the program in various ways (e.g. printing out the value of y2 before comparing it) or by compiling unoptimized.
State-Changed-From-To: open->closed State-Changed-Why: See any faq on numerical analysis that mentions the x86. You are seeing the results of excess precision in the FPU. Either change the rounding precision in the FPCR, or work around the problem with -ffloat-store.
From: rth@gcc.gnu.org To: gcc-gnats@gcc.gnu.org, mirtich@merl.com, nobody@gcc.gnu.org Cc: Subject: Re: optimization/323 Date: 16 Jan 2001 14:30:04 -0000 Synopsis: optimized code gives strange floating point results State-Changed-From-To: open->closed State-Changed-By: rth State-Changed-When: Tue Jan 16 06:30:04 2001 State-Changed-Why: See any faq on numerical analysis that mentions the x86. You are seeing the results of excess precision in the FPU. Either change the rounding precision in the FPCR, or work around the problem with -ffloat-store. http://gcc.gnu.org/cgi-bin/gnatsweb.pl?cmd=view&pr=323&database=gcc
*** Bug 11394 has been marked as a duplicate of this bug. ***
*** Bug 629 has been marked as a duplicate of this bug. ***
*** Bug 4251 has been marked as a duplicate of this bug. ***
*** Bug 6900 has been marked as a duplicate of this bug. ***
*** Bug 10417 has been marked as a duplicate of this bug. ***
Reopening so I can mark it as ...
Invalid and so that bugzilla can pick up it for the dup count.
*** Bug 10644 has been marked as a duplicate of this bug. ***
*** Bug 2582 has been marked as a duplicate of this bug. ***
*** Bug 8606 has been marked as a duplicate of this bug. ***
*** Bug 6824 has been marked as a duplicate of this bug. ***
*** Bug 11518 has been marked as a duplicate of this bug. ***
*** Bug 11618 has been marked as a duplicate of this bug. ***
*** Bug 11655 has been marked as a duplicate of this bug. ***
*** Bug 11767 has been marked as a duplicate of this bug. ***
*** Bug 11892 has been marked as a duplicate of this bug. ***
*** Bug 11914 has been marked as a duplicate of this bug. ***
*** Bug 11934 has been marked as a duplicate of this bug. ***
*** Bug 12129 has been marked as a duplicate of this bug. ***
*** Bug 12331 has been marked as a duplicate of this bug. ***
*** Bug 13265 has been marked as a duplicate of this bug. ***
*** Bug 13890 has been marked as a duplicate of this bug. ***
*** Bug 14367 has been marked as a duplicate of this bug. ***
*** Bug 14384 has been marked as a duplicate of this bug. ***
*** Bug 14652 has been marked as a duplicate of this bug. ***
*** Bug 14989 has been marked as a duplicate of this bug. ***
*** Bug 15134 has been marked as a duplicate of this bug. ***
*** Bug 15437 has been marked as a duplicate of this bug. ***
*** Bug 15852 has been marked as a duplicate of this bug. ***
*** Bug 16223 has been marked as a duplicate of this bug. ***
*** Bug 16607 has been marked as a duplicate of this bug. ***
*** Bug 15386 has been marked as a duplicate of this bug. ***
*** Bug 16686 has been marked as a duplicate of this bug. ***
*** Bug 17802 has been marked as a duplicate of this bug. ***
Excessive precision should be discarded before performing any comparison. It is done e.g. by the Intel fortran compiler. -ffloat-store rounds only the values stored in a variable, and in some situations may trigger instead of suppressing the bug. Existing numerical software packages (e.g. LAPACK and SLATEC fortran libraries) fail their tests with gcc just due to this issue. So it is a bug, even if almost all other compilers reproduce it. If you refer me to a FAQ on numerical analysis, please note that the FAQ states things as they are, and not as they should be.
Some musings based on my encounter with this bug... Comparison of doubles is not just restricted to numerical analysis. Bug report 17802 (aka #323) arose out of a binary search of a timer queue. Is it wrong to assume that if "int secs = 2" and "int usecs = 100000" then "secs + 1e-6 * usecs" will return the same result every time it is evaluated as an inline function? Calling a trivial function twice with an identical parameter may return a value less than, equal to or greater than itself (depending on register allocation in the context of the calls). This is a really nasty piece of nondeterministic behaviour. Weren't compilers and "high level" languages invented to get around exactly this kind of hardware dependency? Compiling with -ffloat-store did not help in my case. Compiling without -O isn't always possible either. Since the function that causes the problem is the trivial calculation above, inlined for performance, I can't guarrantee that every user of the header file it appears in will know to compile without -O (or be willing to). The outcome will be users reporting a bug in my code. This bug is being reported regularly, and probably affecting a lot more people that don't report it.
*** Bug 18567 has been marked as a duplicate of this bug. ***
*** Bug 18756 has been marked as a duplicate of this bug. ***
To summarize, this defect effectively states that: assert( (x/y) == (x/y) ) may cause an assertion if compiled with optimization. While I understand why it happens, that doesn't mean it isn't a defect. This makes it impossible to turn on the optimizer with any code using floating point and still expect to get a correct result. Perhaps in some situations this is okay, but in general this is not. This would also mean the following are also invalid code -- which I'm fairly certain the C/C++ standards would state otherwise: a = (x/y); assert( a == x/y ) //may Abort if( a == x/y ) assert( a == x/y ) //may Abort
*** Bug 18784 has been marked as a duplicate of this bug. ***
*** Bug 19011 has been marked as a duplicate of this bug. ***
*** Bug 19469 has been marked as a duplicate of this bug. ***
*** Bug 19675 has been marked as a duplicate of this bug. ***
*** Bug 19837 has been marked as a duplicate of this bug. ***
*** Bug 20026 has been marked as a duplicate of this bug. ***
*** Bug 20674 has been marked as a duplicate of this bug. ***
*** Bug 7719 has been marked as a duplicate of this bug. ***
(In reply to comment #2) > State-Changed-From-To: open->closed > State-Changed-Why: See any faq on numerical analysis that mentions the x86. > You are seeing the results of excess precision in the FPU. > Either change the rounding precision in the FPCR, or work > around the problem with -ffloat-store. > I had this bug on x86 architecture, with no optimization of the code (no -OX) and with float-store on. My workaround was to store the return of the double function in a auxliar double variable before comparison. Have you an other suggestion ?
(In reply to comment #59) > > I had this bug on x86 architecture, with no optimization of the code (no -OX) > and with float-store on. My workaround was to store the return of the double > function in a auxliar double variable before comparison. > Have you an other suggestion ? > The way I've "fixed" (more like avoided) this problem is to have: #include <fpu_control.h> void set_math_double_precision() { fpu_control_t fpu_control = 0x027f ; _FPU_SETCW(fpu_control); } and make sure this function is called before doing any FP operations. It only needs to be called once. Richard.
*** Bug 21809 has been marked as a duplicate of this bug. ***
Reopening...
...to end this pointless discussion. Some people call this a bug in the x87 series. Other call it a bug in gcc. See these mails at least for the reason why this could be considered a bug in gcc: http://gcc.gnu.org/ml/gcc/2003-08/msg01195.html http://gcc.gnu.org/ml/gcc/2003-08/msg01234.html http://gcc.gnu.org/ml/gcc/2003-08/msg01257.html Regardless of where one wishes to put the blame, this problem will _not_ be fixed. Period.
*** Bug 1098 has been marked as a duplicate of this bug. ***
*** Bug 22394 has been marked as a duplicate of this bug. ***
*** Bug 23298 has been marked as a duplicate of this bug. ***
*** Bug 23407 has been marked as a duplicate of this bug. ***
*** Bug 23318 has been marked as a duplicate of this bug. ***
*** Bug 24014 has been marked as a duplicate of this bug. ***
*** Bug 24387 has been marked as a duplicate of this bug. ***
*** Bug 24479 has been marked as a duplicate of this bug. ***
*** Bug 7935 has been marked as a duplicate of this bug. ***
*** Bug 25303 has been marked as a duplicate of this bug. ***
*** Bug 26000 has been marked as a duplicate of this bug. ***
*** Bug 28191 has been marked as a duplicate of this bug. ***
*** Bug 16825 has been marked as a duplicate of this bug. ***
*** Bug 29597 has been marked as a duplicate of this bug. ***
*** Bug 30255 has been marked as a duplicate of this bug. ***
The option -ffloat-store, recommended by Richard Henderson, has the effect of decreasing the performance of floating-point operations for the entire compilation unit. If you want a minimal fix that does not affect other functions in the same compilation unit, you can use 'volatile double' instead of 'double'. It's like a one-shot -ffloat-store. Example: #include <stdio.h> void test(double x, double y) { const volatile double y2 = x + 1.0; if (y != y2) printf("error\n"); } void main() { const double x = .012; const double y = x + 1.0; test(x, y); }
> const volatile double y2 = x + 1.0; I'd also favor this "selective" approach, because the instability is harmless in most cases. But it is dangerous sometimes, as mentioned in the binary search or when sorting, where the faulty cmpfn() could turn the sort function to infinite loop. So, could you please advise a body of hypothetical double discard_extended_precision(double a); function (possibly some 387 FP insn?) that reliably truncates the argument to 64bit? Would inline double discard_extended_precision(volatile double a) { return a; } do the trick?
*** Bug 30752 has been marked as a duplicate of this bug. ***
*** Bug 31008 has been marked as a duplicate of this bug. ***
*** Bug 31114 has been marked as a duplicate of this bug. ***
I'd like to welcome the newest members of the bug 323 community, where all x87 floating point errors in gcc come to die! All floating point errors that use the x87 are welcome, despite the fact that many of them are easily fixable, and many are not! We're all one happy family, making the egregious mistake of wanting accuracy out of the most accurate general purpose FPU on the market! Cheers, Clint
I think that Uros' patch to add a -mpc switch for precision control would "fix" this. The real fix would be to automatically insert fldcw instructions before float/double operations, in order to limit the precision of the operations. But I think that it would kill speed even more than -ffloat-store.
> I think that Uros' patch to add a -mpc switch for precision control would > "fix" this. > The real fix would be to automatically insert fldcw instructions before > float/double operations, in order to limit the precision of the operations. > But I think that it would kill speed even more than -ffloat-store. Unfortunately, it is not that simple. The -mpc switch and the fldcw instructions control the size of the significant, but they don't control the range of the exponent. So it will solve the issue with the first testcase of this bug-report, but you could still build examples where two execution paths that perform the same floating-point computations produce completely different results.
The following paper explains how this kind of behaviour occurs, why it is "correct", why it is difficult to fix but how it can be partly fixed, and how this breaks many testing and proving techniques: http://hal.archives-ouvertes.fr/hal-00128124
*** Bug 32391 has been marked as a duplicate of this bug. ***
*** Bug 32976 has been marked as a duplicate of this bug. ***
*** Bug 33611 has been marked as a duplicate of this bug. ***
*** Bug 34616 has been marked as a duplicate of this bug. ***
*** Bug 34951 has been marked as a duplicate of this bug. ***
*** Bug 34992 has been marked as a duplicate of this bug. ***
*** Bug 35389 has been marked as a duplicate of this bug. ***
Subject: RE: optimized code gives strange floating point results Not sure this is the same issues as 323. All three numbers, 8, 1 and 65, should be able to represented in IEEE 754 floating-point format exactly without any rounding or approximation. That is 8 = 1* 2^3 1 = 1* 2^0 65 = (1 + 1/64) * 2^6 -----Original Message----- From: pinskia at gcc dot gnu dot org [mailto:gcc-bugzilla@gcc.gnu.org] Sent: Wednesday, February 27, 2008 12:39 PM To: Wei, Yongbin Subject: [Bug rtl-optimization/323] optimized code gives strange floating point results ------- Comment #103 from pinskia at gcc dot gnu dot org 2008-02-27 20:38 ------- *** Bug 35389 has been marked as a duplicate of this bug. *** -- pinskia at gcc dot gnu dot org changed: What |Removed |Added ------------------------------------------------------------------------ ---- CC| |ywei at qualcomm dot com http://gcc.gnu.org/bugzilla/show_bug.cgi?id=323 ------- You are receiving this mail because: ------- You are on the CC list for the bug, or are watching someone who is.
*** Bug 35488 has been marked as a duplicate of this bug. ***
*** Bug 35489 has been marked as a duplicate of this bug. ***
*** Bug 35585 has been marked as a duplicate of this bug. ***
I also encountered such problems and was going to report it as a bug in GCC... But in the GCC bug (not) reporting guide, there is fortunately a link to this page and here (comment #96) is a link to David Monniaux's paper about floating-point computations. This explains it closely but it is maybe too long. I have almost read it and hope I have understood it properly. So I'll give a brief explanation (for those who don't know it yet) of the reasons of such a strange behaviour. Then I'll assess where the bug actually is (in GCC or CPU). Then I'll write the solution (!) and finally a few recommendations to the GCC team. EXPLANATION The x87 FPU was originally designed in (or before) 1980. I think that's why it is quite simple: it has only one unit for all FP data types. Of course, the precision must be of the widest type, which is the 80-bit long double. Consider you have a program, where all the FP variables are of the type double. You perform some FP operations and one of them is e.g. 1e-300/1e300, which results in 1e-600. Despite this value cannot be held by a "double", it is stored in an 80-bit FPU register as the result. Consider you use the variable "x" to hold that result. If the program has been compiled with optimization, the value need not be stored in RAM. So, say, it is still in the register. Consider you need x to be nonzero, so you perform the test x != 0. Since 1e-600 is not zero, the test yields true. While you perform some other computations, the value is moved to RAM and converted to 0 because x is of type "double". Now you want to use your certainly nonzero x... Hard luck :-( Note that if the result doesn't have its corresponding variable and you perform the test directly on an expression, the problem can come to light even without optimization. It could seem that performing all FP operations in extended precision can bring benefits only. But it introduces a serious pitfall: moving a value may change the value!!! WHERE'S THE BUG This is really not a GCC bug. The bug is actually in the x87 FPU because it doesn't obey the IEEE standard. SOLUTION The x87 FPU is still present in contemporary processors (including AMD) due to compatibility. I think most of PC software still uses it. But new processors have also another FPU, called SSE, and this do obey the IEEE. GCC in 32-bit mode compiles for x87 by default but it is able to compile for the SSE, too. So the solution is to add these options to the compilation command: -march=* -msse -mfpmath=sse Yes, this definitely resolves the problem - but not for all processors. The * can be one of the following: pentium3, pentium3m, pentium-m, pentium4, pentium4m, prescott, nocona, athlon-4, athlon-xp, athlon-mp, k8, opteron, athlon64, athlon-fx and c3-2 (I'm unsure about athlon and athlon-tbird). Beside -msse, you can also add some of -mmmx, -msse2, -msse3 and -m3dnow, if the CPU supports them (see GCC doc or CPU doc). If you wish to compile for processors which don't have SSE, you have a few possibilities: (1) A very simple solution: Use long double everywhere. (But be careful when transfering binary data in long double format between computers because this format is not standardized and so the concrete bit representations vary between different CPU architectures.) (2) A partial but simple solution: Do comparisons on volatile variables only. (3) A similar solution: Try to implement a "discard_extended_precision" function suggested by Egon in comment #88. (4) A complex solution: Before doing any mathematical operation or comparison, put the operands into variables and put also the result to a variable (i.e. don't use complex expressions). For example, instead of { c = 2*(a+b); } , write { double s = a+b; c = 2*s; } . I'm unsure about arrays but I think they should be OK. When you have modified your code in this manner, then compile it either without optimization or, when using optimization, use -ffloat-store. In order to avoid double rounding (i.e. rounding twice), it is also good to decrease the FPU precision by changing its control word in the beginning of your program (see comment #60). Then you should also apply -frounding-math. (5) A radical solution: Find a job/hobby where computers are not used at all. RECOMMENDATIONS I think this problem is really serious and general. Therefore, programmers should be warned soon enough. This recommendation should be addressed especially to authors of programming coursebooks. But I think there could also be a paragraph about it in the GCC documentation (I haven't read it wholly but it doesn't seem there's any warning against x87). And, of course, there should be a warning in the bug reporting guide (http://gcc.gnu.org/bugs.html). It's fine there's a link to this page (Bug 323) but the example with (int)(a/b) is insufficient. It only demonstrates that real numbers are often not represented exactly in the computer. It doesn't demonstrate the x87 pitfall. Hence there should be an example such as the initial code of this "GCC bug 323 report". Because when one sees the example with (int)(a/b), he can say "It's trivial" and not click the link (as I did the first time). EPILOGUE I hope my effort of writing this "comment #109" will be helpful for many people. If you want more info, read the David Monniaux's work or something else about FPUs. Thanks to David Monniaux.
I used an old version of GCC documentation so I omitted some new processors with SSE: core2, k8-sse3, opteron-sse3, athlon64-sse3, amdfam10 and barcelona. I think you can use -march=pentium3 for all Intel's CPUs (of course, starting with P3). I'm unsure about AMD. (Maybe you know it better.)
(In reply to comment #109) > WHERE'S THE BUG > This is really not a GCC bug. The bug is actually in the x87 FPU because it > doesn't obey the IEEE standard. Concerning the standards: The x87 FPU does obey the IEEE754-1985 standard, which *allows* extended precision, and double precision is *available*. In fact, one could say that GCC even obeys the IEEE standard (which doesn't define bindings: the definition of "destination" page 4 of the IEEE754-1985 standard is rather vague and lets the language to define it exactly), but it doesn't obey the ISO C99 standard on some point. Concerning the x87 FPU: One can say however that the x87 is a badly designed because it is not possible to statically specify the precision. Nevertheless the OS/language implementations should take care of this problem. Note: the solution chosen by some OS'es (*BSD, MS-Windows...) is to configure the processor to the IEEE double precision by default (thus "long double" is also in double precision, but this is OK as far as the C language is concerned, there's still a problem with "float", but in practice, nobody cares AFAIK). > If you wish to compile for processors which don't have SSE, you have a few > possibilities: > (1) A very simple solution: Use long double everywhere. This avoids the bug, but this is not possible for software that requires double precision exactly, e.g. XML tools that use XPath. See other examples here: http://www.vinc17.org/research/extended.en.html Also this makes maintenance of software more difficult because long double can be much slower on some platforms, which support this type in software to provide more precision (e.g. PowerPC Linux and Mac OS X implement a double-double arithmetic, Solaris and HPUX implement quadruple precision). > (But be careful when transfering binary data in long double format between > computers because this format is not standardized and so the concrete bit > representations vary between different CPU architectures.) Well, this is not specific to long double anyway: there exist 3 possible endianess for the double format (x86, PowerPC, ARM). > (2) A partial but simple solution: Do comparisons on volatile variables only. Yes (but this is also a problem concerning the maintenance of portable programs). > (4) A complex solution: [...] Yes, this is the workaround I use in practice. > RECOMMENDATIONS > I think this problem is really serious and general. Therefore, programmers > should be warned soon enough. Yes, but note that this is not the only problem with compilers. See e.g. http://gcc.gnu.org/bugzilla/show_bug.cgi?id=36578 for a bug related to casts to long double on x86_64 and ia64. This one is now tested by: http://www.vinc17.org/software/tst-ieee754.c (which has also tested bug 323 for a long time).
(In reply to comment #111) > Concerning the standards: The x87 FPU does obey the IEEE754-1985 standard, > which *allows* extended precision, and double precision is *available*. It's true that double *precision* is available on x87. But not the *IEEE-754 "double precision" type*. Beside the precision of mantissa, this includes also the range of exponent. On the x87, it is possible to set the precision of mantissa but not the range of exponent. That's why I believe it doesn't obey the IEEE. (I haven't ever seen the IEEE-754 standard but I base on the work of David Monniaux.) > Note: the solution chosen by some OS'es (*BSD, MS-Windows...) is to configure > the processor to the IEEE double precision by default (thus "long double" is > also in double precision, but this is OK as far as the C language is concerned, > there's still a problem with "float", but in practice, nobody cares AFAIK). Do you mean that on Windows, long double has (by default) no more precision than double? I don't think so (it's confirmed by my experience). According to the paper of David Monniaux, only FreeBSD 4 sets double precision by default (but I know almost nothing about BSD). > > (1) A very simple solution: Use long double everywhere. > This avoids the bug, but this is not possible for software that requires double > precision exactly, e.g. XML tools that use XPath. Yes, of course. I don't say this can be used everywhere. > > (But be careful when transfering binary data in long double format between > > computers because this format is not standardized and so the concrete bit > > representations vary between different CPU architectures.) > Well, this is not specific to long double anyway: there exist 3 possible > endianess for the double format (x86, PowerPC, ARM). OK but David Monniaux mentions portability issues just in the case of long double, so the differences are probably more frequent in this case (maybe even within the x86 architecture). > Yes, but note that this is not the only problem with compilers. See e.g. > http://gcc.gnu.org/bugzilla/show_bug.cgi?id=36578 Thanks for info.
(In reply to comment #112) > It's true that double *precision* is available on x87. But not the *IEEE-754 > "double precision" type*. It is available when storing a result to memory. > Beside the precision of mantissa, this includes also the range of exponent. > On the x87, it is possible to set the precision of mantissa but not the range > of exponent. The IEEE754-1985 allows this. Section 4.3: "Normally, a result is rounded to the precision of its destination. However, some systems deliver results only to double or extended destinations. On such a system the user, which may be a high-level language compiler, shall be able to specify that a result be rounded instead to single precision, though it may be stored in the double or extended format with its wider exponent range. [...]" > That's why I believe it doesn't obey the IEEE. (I haven't ever seen the > IEEE-754 standard but I base on the work of David Monniaux.) See above. Also beware of subtilities in the wording used by David Monniaux. FYI, the IEEE754-1985 standard (with minor corrections) is available from the following page: http://www.validlab.com/754R/ (look at the end). AFAIK, the IEEE754-1985 standard was designed from the x87 implementation, so it would have been very surprising that x87 didn't conform to IEEE754-1985. > Do you mean that on Windows, long double has (by default) no more precision > than double? I don't think so (it's confirmed by my experience). I don't remember my original reference, but here's a new one: http://msdn.microsoft.com/en-us/library/aa289157(vs.71).aspx In fact, this depends on the architecture. I quote: "x86. Intermediate expressions are computed at the default 53-bit precision with an extended range provided by a 16-bit exponent. When these 53:16 values are "spilled" to memory (as can happen during a function call), the extended exponent range will be narrowed to 11-bits. That is, spilled values are cast to the standard double precision format with only an 11-bit exponent. A user may switch to extended 64-bit precision for intermediate rounding by altering the floating-point control word using _controlfp and by enabling FPU environment access (see The fpenv_access Pragma). However, when extended precision register-values are spilled to memory, the intermediate results will still be rounded to double precision. This particular semantic is subject to change." Note that the behavior has changed in some version of Windows (it was using the extended precision, then it switched to double precision for x86). Now, this may also depend on the compiler. > According to the paper of David Monniaux, only FreeBSD 4 sets double > precision by default (but I know almost nothing about BSD). I've noted that amongst the BSD's, NetBSD does this too (I don't remember if I've tried or got it from some document, and this might also depend on the NetBSD version and/or the processor).
(In reply to comment #113) > It is available when storing a result to memory. Yes, but this requires quite a complicated workaround (solution (4) in my comment #109). So you could say that the IEEE754 double precision type is available even on a processor without any FPU because this can be emulated using integers. Moreover, if we assess things pedantically, the workaround (4) still doesn't fully obey the IEEE single/double precision type(s), because there remains the problem of double rounding of denormals. > The IEEE754-1985 allows this. Section 4.3: "Normally, a result is rounded to > the precision of its destination. However, some systems deliver results only to > double or extended destinations. On such a system the user, which may be a > high-level language compiler, shall be able to specify that a result be rounded > instead to single precision, though it may be stored in the double or extended > format with its wider exponent range. [...]" > [...] > AFAIK, the IEEE754-1985 standard was designed from the x87 > implementation, so it would have been very surprising that x87 didn't conform > to IEEE754-1985. So it seems I was wrong but the IEEE754-1985 standard is also quite "wrong". > > Do you mean that on Windows, long double has (by default) no more precision > > than double? I don't think so (it's confirmed by my experience). > I don't remember my original reference, but here's a new one: > http://msdn.microsoft.com/en-us/library/aa289157(vs.71).aspx > In fact, this depends on the architecture. I quote: "x86. Intermediate > expressions are computed at the default 53-bit precision with an extended range > [...]" I quote, too: "Applies To Microsoft® Visual C++®"
That ® should be (R).
(In reply to comment #114) > Yes, but this requires quite a complicated workaround (solution (4) in my > comment #109). The problem is on the compiler side, which could store every result of a cast or an assignment to memory (this is inefficient, but that's what you get with the x87, and the ISO C language could be blamed too for *requiring* something like that instead of being more flexible). > So you could say that the IEEE754 double precision type is available even on > a processor without any FPU because this can be emulated using integers. Yes, but a conforming implementation would be the processor + a library, not just the processor with its instruction set. > Moreover, if we assess things pedantically, the workaround (4) still doesn't > fully obey the IEEE single/double precision type(s), because there remains the > problem of double rounding of denormals. As I said, in this particular case (underflow/overflow), double rounding is allowed by the IEEE standard. It may not be allowed by some languages (e.g. XPath, and Java in some mode) for good or bad reasons, but this is another problem. > I quote, too: > "Applies To > Microsoft® Visual C++®" Now I assume that it follows the MS-Windows API (though nothing is certain with Microsoft). And the other compilers under MS-Windows could (or should) do the same thing.
(In reply to comment #116) > > Yes, but this requires quite a complicated workaround (solution (4) in my > > comment #109). > > The problem is on the compiler side, which could store every result of a cast > or an assignment to memory (this is inefficient, but that's what you get with > the x87, and the ISO C language could be blamed too for *requiring* something > like that instead of being more flexible). > > > So you could say that the IEEE754 double precision type is available even on > > a processor without any FPU because this can be emulated using integers. > > Yes, but a conforming implementation would be the processor + a library, not > just the processor with its instruction set. > > > Moreover, if we assess things pedantically, the workaround (4) still doesn't > > fully obey the IEEE single/double precision type(s), because there remains the > > problem of double rounding of denormals. > > As I said, in this particular case (underflow/overflow), double rounding is > allowed by the IEEE standard. It may not be allowed by some languages (e.g. > XPath, and Java in some mode) for good or bad reasons, but this is another > problem. OK, thanks for explanation. I think now it's clear. > > I quote, too: > > "Applies To > > Microsoft® Visual C++®" > > Now I assume that it follows the MS-Windows API (though nothing is certain with > Microsoft). And the other compilers under MS-Windows could (or should) do the > same thing. By a lucky hit, I have found this in the GCC documentation: " -mpc32 -mpc64 -mpc80 Set 80387 floating-point precision to 32, 64 or 80 bits. When '-mpc32' is specified, the significands of results of floating-point operations are rounded to 24 bits (single precision); '-mpc64' rounds the the significands of results of floatingpoint operations to 53 bits (double precision) and '-mpc80' rounds the significands of results of floating-point operations to 64 bits (extended double precision), which is the default. When this option is used, floating-point operations in higher precisions are not available to the programmer without setting the FPU control word explicitly. [...]" So GCC sets extended precision by default. And it's easy to change it.
(In reply to comment #117) > By a lucky hit, I have found this in the GCC documentation: > " > -mpc32 > -mpc64 > -mpc80 OK, this is new in gcc 4.3. I haven't tried, but if gcc just changes the precision without changing the values of <float.h> macros to make them correct, this is just a workaround (better than nothing, though). Also, this is a problem for library code if it requires to have double precision instead of extended precision, as these options won't probably be taken into account at that point. (Unfortunately it's probably too late to have a clean ABI.)
Another example related with fp on x87? EXPECTED RESULT: 0 (with EPS accuracy) 0 (with EPS accuracy) 0 (with EPS accuracy) 0 (with EPS accuracy) REAL RESULT: 5.313991e+33 5.313991e+33 0.000000e+00 0.000000e+00 CODE #include <stdio.h> int main( void ) { /* register */ double d1 = 1e50; /* register */ double d2 = -2.7438011834107752e+51; /* register */ double s = 0.036445789368634796; /* register */ double d3 = -d1/s; /* register */ double d4 = s*d2; /* register */ double d5 = s*d3; printf( "%e\n", d1 + s*d2); printf( "%e\n", d1 + s*d3); printf( "%e\n", d1 + d4); printf( "%e\n", d1 + d5); return 0; }
(In reply to comment #119) > REAL RESULT: > 5.313991e+33 > 5.313991e+33 > 0.000000e+00 > 0.000000e+00 Only without optimizations. But since the ISO C standard allows expressions to be evaluated in a higher precision, there's no bug here (unless you show a contradiction with the value of FLT_EVAL_METHOD, but the FP_CONTRACT pragma should also be set to OFF -- though this currently has no effect on gcc).
Thank you Vincent. I fact after commenting I realized that this is a plain numerical error on the last digit of double in multiplication. I think that my comment was rather irrelevant and I am the more ashamed the more I cannot remove it from bugzilla :)
*** Bug 37390 has been marked as a duplicate of this bug. ***
Subject: Bug 323 Author: jsm28 Date: Tue Nov 4 13:24:30 2008 New Revision: 141578 URL: http://gcc.gnu.org/viewcvs?root=gcc&view=rev&rev=141578 Log: PR rtl-optimization/323 * c-common.c (convert_and_check, c_common_truthvalue_conversion): Handle EXCESS_PRECISION_EXPR. * c-common.def (EXCESS_PRECISION_EXPR): New. * c-cppbuiltin.c (builtin_define_float_constants): Define constants with enough digits for long double. * c-lex.c (interpret_float): Interpret constant with excess precision where appropriate. * c-opts.c (c_common_post_options): Set flag_excess_precision_cmdline. * c-parser.c (c_parser_conditional_expression): Handle excess precision in condition. * c-typeck.c (c_fully_fold): Handle EXCESS_PRECISION_EXPR. (c_fully_fold_internal): Disallow EXCESS_PRECISION_EXPR. (convert_arguments): Handle arguments with excess precision. (build_unary_op): Move excess precision outside operation. (build_conditional_expr): Likewise. (build_compound_expr): Likewise. (build_c_cast): Do cast on operand of EXCESS_PRECISION_EXPR. (build_modify_expr): Handle excess precision in RHS. (convert_for_assignment): Handle excess precision in converted value. (digest_init, output_init_element, process_init_element): Handle excess precision in initializer. (c_finish_return): Handle excess precision in return value. (build_binary_op): Handle excess precision in operands and add excess precision as needed for operation. * common.opt (-fexcess-precision=): New option. * config/i386/i386.h (X87_ENABLE_ARITH, X87_ENABLE_FLOAT): New. * config/i386/i386.md (float<SSEMODEI24:mode><X87MODEF:mode>2): For standard excess precision, output explicit conversion to and truncation from XFmode. (*float<SSEMODEI24:mode><X87MODEF:mode>2_1, *float<SSEMODEI24:mode><X87MODEF:mode>2_i387_with_temp, *float<SSEMODEI24:mode><X87MODEF:mode>2_i387, two unnamed define_splits, floatdi<X87MODEF:mode>2_i387_with_xmm, two unnamed define_splits, *floatunssi<mode>2_1, two unnamed define_splits, floatunssi<mode>2, add<mode>3, sub<mode>3, mul<mode>3, divdf3, divsf3, *fop_<mode>_comm_i387, *fop_<mode>_1_i387, *fop_<MODEF:mode>_2_i387, *fop_<MODEF:mode>_3_i387, *fop_df_4_i387, *fop_df_5_i387, *fop_df_6_i387, two unnamed define_splits, sqrt<mode>2): Disable where appropriate for standard excess precision. * convert.c (convert_to_real): Do not shorten arithmetic to type for which excess precision would be used. * doc/invoke.texi (-fexcess-precision=): Document option. (-mfpmath=): Correct index entry. * flags.h (enum excess_precision, flag_excess_precision_cmdline, flag_excess_precision): New. * langhooks.c (lhd_post_options): Set flag_excess_precision_cmdline. * opts.c (common_handle_option): Handle -fexcess-precision=. * toplev.c (flag_excess_precision_cmdline, flag_excess_precision, init_excess_precision): New. (lang_dependent_init_target): Call init_excess_precision. * tree.c (excess_precision_type): New. * tree.h (excess_precision_type): Declare. ada: * gcc-interface/misc.c (gnat_post_options): Set flag_excess_precision_cmdline. fortran: * options.c (gfc_post_options): Set flag_excess_precision_cmdline. java: * lang.c (java_post_options): Set flag_excess_precision_cmdline. testsuite: * gcc.target/i386/excess-precision-1.c, gcc.target/i386/excess-precision-2.c, gcc.target/i386/excess-precision-3.c, gcc.target/i386/excess-precision-4.c, gcc.target/i386/excess-precision-5.c, gcc.target/i386/excess-precision-6.c: New tests. Added: branches/c-4_5-branch/gcc/ada/ChangeLog.c45 branches/c-4_5-branch/gcc/fortran/ChangeLog.c45 branches/c-4_5-branch/gcc/java/ChangeLog.c45 branches/c-4_5-branch/gcc/testsuite/gcc.target/i386/excess-precision-1.c branches/c-4_5-branch/gcc/testsuite/gcc.target/i386/excess-precision-2.c branches/c-4_5-branch/gcc/testsuite/gcc.target/i386/excess-precision-3.c branches/c-4_5-branch/gcc/testsuite/gcc.target/i386/excess-precision-4.c branches/c-4_5-branch/gcc/testsuite/gcc.target/i386/excess-precision-5.c branches/c-4_5-branch/gcc/testsuite/gcc.target/i386/excess-precision-6.c Modified: branches/c-4_5-branch/gcc/ChangeLog.c45 branches/c-4_5-branch/gcc/ada/gcc-interface/misc.c branches/c-4_5-branch/gcc/c-common.c branches/c-4_5-branch/gcc/c-common.def branches/c-4_5-branch/gcc/c-cppbuiltin.c branches/c-4_5-branch/gcc/c-lex.c branches/c-4_5-branch/gcc/c-opts.c branches/c-4_5-branch/gcc/c-parser.c branches/c-4_5-branch/gcc/c-typeck.c branches/c-4_5-branch/gcc/common.opt branches/c-4_5-branch/gcc/config/i386/i386.h branches/c-4_5-branch/gcc/config/i386/i386.md branches/c-4_5-branch/gcc/convert.c branches/c-4_5-branch/gcc/doc/invoke.texi branches/c-4_5-branch/gcc/flags.h branches/c-4_5-branch/gcc/fortran/options.c branches/c-4_5-branch/gcc/java/lang.c branches/c-4_5-branch/gcc/langhooks.c branches/c-4_5-branch/gcc/opts.c branches/c-4_5-branch/gcc/testsuite/ChangeLog.c45 branches/c-4_5-branch/gcc/toplev.c branches/c-4_5-branch/gcc/tree.c branches/c-4_5-branch/gcc/tree.h
Vincent Lefèvre is right: the issue is quite subtle. (I should mention that Vincent is an expert in computer arithmetics, which I'm not.) As he rightly points, conformance to IEEE-754 should be evaluated for a whole software/hardware system - it is possible to have a IEEE-754 system entirely implemented in software. It seems like the C99 standard prohibits double rounding, and prohibits having values depend on the vagaries of register spilling. Except that this prohibition is explicit only in non-normative sections. "Language lawyers" have sent me justifications that this prohibition is implied by various normative prescriptions of the standard. I think that in any case we should not spend too much energy trying to assign blame (who conforms to the standard) but rather try to find solutions. The short answer is that no compiler, be it gcc, will be modified so that complex sequences of operations are used for floating-point operations in lieu of directly using x87 instructions! At least for two reasons: * x87 is now fading away (its use is deprecated on x86-64, it's not used by default on Intel Macintosh...) * Most people don't want to pay the performance hit. In addition, I think there are more urgent things to fix in gcc's floating-point system, such as support for #pragma STDC FENV ACCESS Now for some additional facts: * It is possible to force the x87 to use reduced precision for the mantissa (with inline asm or even now with gcc options). * This setting does not affect the range of exponents. so you can still have surprises if a very tiny nonzero value is kept in a register then is rounded to zero when spilled to memory. * In some rare cases, you can have "double rounding on underflow": you get a different result by computing on SSE in double precision mode on the one hand, and by computing on x87 in "double precision" then writing to a double variable in memory.
(In reply to comment #124) > It seems like the C99 standard prohibits double rounding, only if Annex F is claimed to be supported (note: Annex F is not just IEEE 754, it also contains specific bindings). IEEE 754 doesn't prohibit double rounding either (this depends on the bindings), but with C99 + Annex F, double rounding is prohibited. Now, bug 323 is not about double rounding specifically. There are two potential problems: 1. A "double" variable (or result of a cast) contains a "long double" value (not exactly representable in a "double"). This is prohibited by C99 (5.1.2.3#12, 6.3.1.5#2 and 6.3.1.8#2[52]). This problem seems to be fixed by Joseph Myers' patch mentioned in comment #123 (but I haven't tried). 2. Computations on "double" expressions are carried out in extended precision. This is allowed by C99 (except for casts and assignments), e.g. when FLT_EVAL_METHOD=2. But if the implementation (i.e. here compiler + library + ...) claims to support Annex F, then this is prohibited. This point is rather tricky because the compiler (GCC) and library (e.g. GNU libc) settings must be consistent, so their developers need to talk with each other. FYI, I reported the following bug concerning glibc: http://sourceware.org/bugzilla/show_bug.cgi?id=6981 because it sets __STDC_IEC_559__ to 1 unconditionally. > The short answer is that no compiler, be it gcc, will be modified so that > complex sequences of operations are used for floating-point operations in lieu > of directly using x87 instructions! At least for two reasons: > * x87 is now fading away (its use is deprecated on x86-64, it's not used by > default on Intel Macintosh...) > * Most people don't want to pay the performance hit. That's why in Joseph's patch, it's just an option (disabled by default, but enabled by -std=c99 because one should assume that if a user wants C99, then he really wants it, and if he is able to add an option, then he is also able to add another one if he wants to disable this fix in case he knows it is useless for his application -- this is also true for -ffast-math). GCC already supports SSE, but this patch is for processors that don't. Also the performance hit depends very much on the application. Performance hit is reduced in applications that do not use intensive FP or mostly interactive applications. > In addition, I think there are more urgent things to fix in gcc's > floating-point system, such as support for #pragma STDC FENV ACCESS FYI, this is bug 34678. And I submitted bug 37845 concerning the FP_CONTRACT pragma. > * It is possible to force the x87 to use reduced precision for the mantissa > (with inline asm or even now with gcc options). Unfortunately, this means that "long double" wouldn't behave as expected, and the behavior is not controllable enough (e.g. due to libraries, plugins...). Such a change should have been system-wide. Now, this is needed in software where double rounding is prohibited (e.g. XSLT processor).
Subject: Bug 323 Author: jsm28 Date: Mon Mar 30 01:50:44 2009 New Revision: 145272 URL: http://gcc.gnu.org/viewcvs?root=gcc&view=rev&rev=145272 Log: PR rtl-optimization/323 * c-common.c (c_fully_fold, convert_and_check, c_common_truthvalue_conversion): Handle EXCESS_PRECISION_EXPR. (c_fully_fold_internal): Disallow EXCESS_PRECISION_EXPR. * c-common.def (EXCESS_PRECISION_EXPR): New. * c-cppbuiltin.c (builtin_define_float_constants): Define constants with enough digits for long double. * c-lex.c (interpret_float): Interpret constant with excess precision where appropriate. * c-opts.c (c_common_post_options): Set flag_excess_precision_cmdline. Give an error for -fexcess-precision=standard for C++ for processors where the option is significant. * c-parser.c (c_parser_conditional_expression): Handle excess precision in condition. * c-typeck.c (convert_arguments): Handle arguments with excess precision. (build_unary_op): Move excess precision outside operation. (build_conditional_expr): Likewise. (build_compound_expr): Likewise. (build_c_cast): Do cast on operand of EXCESS_PRECISION_EXPR. (build_modify_expr): Handle excess precision in RHS. (convert_for_assignment): Handle excess precision in converted value. (digest_init, output_init_element, process_init_element): Handle excess precision in initializer. (c_finish_return): Handle excess precision in return value. (build_binary_op): Handle excess precision in operands and add excess precision as needed for operation. * common.opt (-fexcess-precision=): New option. * config/i386/i386.h (X87_ENABLE_ARITH, X87_ENABLE_FLOAT): New. * config/i386/i386.md (float<SSEMODEI24:mode><X87MODEF:mode>2): For standard excess precision, output explicit conversion to and truncation from XFmode. (*float<SSEMODEI24:mode><X87MODEF:mode>2_1, *float<SSEMODEI24:mode><X87MODEF:mode>2_i387_with_temp, *float<SSEMODEI24:mode><X87MODEF:mode>2_i387, two unnamed define_splits, floatdi<X87MODEF:mode>2_i387_with_xmm, two unnamed define_splits, *floatunssi<mode>2_1, two unnamed define_splits, floatunssi<mode>2, add<mode>3, sub<mode>3, mul<mode>3, divdf3, divsf3, *fop_<mode>_comm_i387, *fop_<mode>_1_i387, *fop_<MODEF:mode>_2_i387, *fop_<MODEF:mode>_3_i387, *fop_df_4_i387, *fop_df_5_i387, *fop_df_6_i387, two unnamed define_splits, sqrt<mode>2): Disable where appropriate for standard excess precision. * convert.c (convert_to_real): Do not shorten arithmetic to type for which excess precision would be used. * defaults.h (TARGET_FLT_EVAL_METHOD_NON_DEFAULT): Define. * doc/invoke.texi (-fexcess-precision=): Document option. (-mfpmath=): Correct index entry. * flags.h (enum excess_precision, flag_excess_precision_cmdline, flag_excess_precision): New. * langhooks.c (lhd_post_options): Set flag_excess_precision_cmdline. * opts.c (common_handle_option): Handle -fexcess-precision=. * toplev.c (flag_excess_precision_cmdline, flag_excess_precision, init_excess_precision): New. (lang_dependent_init_target): Call init_excess_precision. * tree.c (excess_precision_type): New. * tree.h (excess_precision_type): Declare. ada: * gcc-interface/misc.c (gnat_post_options): Set flag_excess_precision_cmdline. Give an error for -fexcess-precision=standard for processors where the option is significant. fortran: * options.c (gfc_post_options): Set flag_excess_precision_cmdline. Give an error for -fexcess-precision=standard for processors where the option is significant. java: * lang.c (java_post_options): Set flag_excess_precision_cmdline. Give an error for -fexcess-precision=standard for processors where the option is significant. testsuite: * gcc.target/i386/excess-precision-1.c, gcc.target/i386/excess-precision-2.c, gcc.target/i386/excess-precision-3.c, gcc.target/i386/excess-precision-4.c, gcc.target/i386/excess-precision-5.c, gcc.target/i386/excess-precision-6.c: New tests. Added: trunk/gcc/testsuite/gcc.target/i386/excess-precision-1.c trunk/gcc/testsuite/gcc.target/i386/excess-precision-2.c trunk/gcc/testsuite/gcc.target/i386/excess-precision-3.c trunk/gcc/testsuite/gcc.target/i386/excess-precision-4.c trunk/gcc/testsuite/gcc.target/i386/excess-precision-5.c trunk/gcc/testsuite/gcc.target/i386/excess-precision-6.c Modified: trunk/gcc/ChangeLog trunk/gcc/ada/ChangeLog trunk/gcc/ada/gcc-interface/misc.c trunk/gcc/c-common.c trunk/gcc/c-common.def trunk/gcc/c-cppbuiltin.c trunk/gcc/c-lex.c trunk/gcc/c-opts.c trunk/gcc/c-parser.c trunk/gcc/c-typeck.c trunk/gcc/common.opt trunk/gcc/config/i386/i386.h trunk/gcc/config/i386/i386.md trunk/gcc/convert.c trunk/gcc/defaults.h trunk/gcc/doc/invoke.texi trunk/gcc/flags.h trunk/gcc/fortran/ChangeLog trunk/gcc/fortran/options.c trunk/gcc/java/ChangeLog trunk/gcc/java/lang.c trunk/gcc/langhooks.c trunk/gcc/opts.c trunk/gcc/testsuite/ChangeLog trunk/gcc/toplev.c trunk/gcc/tree.c trunk/gcc/tree.h
Fixed for C (and ObjC) for 4.5 with the new -fexcess-precision=standard support. The issue remains for other languages (and maybe for some m68k processors); I quote from my original message <http://gcc.gnu.org/ml/gcc-patches/2008-11/msg00105.html>: It would be possible to implement the option for non-C languages, to provide whatever predictable semantics are appropriate for those languages (whether or not described in their standards). Note that bug 323 was originally reported with a C++ testcase. If implemented for all languages, the option might supersede -ffloat-store. Right now, -ffloat-store checks are scattered about the optimizers and it seems unlikely that -ffloat-store really implements any form of predictable semantics now; such semantic effect as it was intended to have could be better represented as an alias for a -fexcess-precision=standard option supported for all languages. It would probably be most appropriate not to close bug 323 without having some form of predictable semantics available for each language. and: I have not changed the m68k back end in this patch. Thus the option may not be fully effective for the affected m68k processors (classic m68k with 68881, before 68040, only, not ColdFire, not 68040 or later). If anyone wishes to make it fully effective for such processors they should copy the testcases to gcc.target/m68k/ and go through m68k insn patterns appropriately adjusting them.
*** Bug 40186 has been marked as a duplicate of this bug. ***
I am a bit wondering if this bug is also for the case (a < b) && (b < a) == true. Is it? Because if so, this becomes way more serious, as for example std::set<double> is broken then (and depending on the STL implementation, it will throw assertions then). If not, I guess my bug #40186 is not a duplicate of this bug.