13 ms·

GCC optimized code gives strange floating point results

145 points by tbrowbdidnso 10y ago

castratikron 10y ago
More proof that floats are evil. Moral of the story is to never use equality to compare floats.
- trendia 10y ago
  Pretty much. If you want to compare two floats a and b within some precision eps, use: abs(a - b) <= eps*0.5*(a + b) or even (abs(a - b) <= eps*abs(a)) && (abs(a - b) <= eps*abs(b)) which is implemented in Boost [0]. [0] http://www.boost.org/doc/libs/1_34_0/libs/test/doc/components/test_tools/floating_point_comparison.html http://www.boost.org/doc/libs/1_34_0/libs/test/doc/component...
  - ratboy666 10y ago
    I thought it was (a / b) != (a / b)
- geezerjay 10y ago
  > More proof that floats are evil. Just because someone doesn't know the basics of floating point numbers it doesn't mean that the floating point format is somehow evil. Testing floats for equality is one of those no-nos which is taught right in the first class of introduction to FP numbers. Another very basic issue covered right in the intro clsses are problems caused by using decimal values to define numbers encoded in binary format, which is often explained through the famous "0.1 doesn't exist in binary" trivia bit.
  - thechao 10y ago
    Part of my day job is to test floating point hardware and write compilers and assemblers. Testing equality---with some care---is completely fine. I can assure you this is a bug in GCC. The GCC bug is compounded by the nonstandard behavior of x87. The x87 unit is not "fixable" but GCC's behavior is.
    - nurettin 10y ago
      That's interesting, you make it sound like the production of new FPUs and new compilers is common. Which industry do you work in?
- nom 10y ago
  Floats are not evil. They are just not as intuitive as integers because it is impossible to represent every possible rational number in the range with a fixed number of bits. Once you get that, everything makes sense.
  - monochromatic 10y ago
    That does not make sense of getting a failure from assert x/y == x/y
    - mikeash 10y ago
      Properly implemented floats will never fail that (barring NaN). That gcc sometimes does is a compiler bug, not a float problem.
    - cellularmitosis 10y ago
      I agree. Explaining that assertion requires explaining that the x87 FPU is providing a leaky abstraction. If you ask for double math, it gives you a result which is a long double while it resides in a register, and then becomes a double when moved to RAM. There are some cases where direct comparison of floating point values is what you want, and this issue can still bite you even then.
      - int_19h 10y ago
        For x87, it's not a leaky abstraction - it fully supports 80-bit floats as a type, and documents that its registers are 80-bit, and so are operations on them. There's nothing abstract about it. The bug is that GCC is utilizing the (documented, by-design) x87 behavior in a way that causes it to produce the wrong result, for the sake of performance.
        abstractbeliefs 10y ago
        >compiler developers do unintuitive but technically legal thing in the name of performance I still think the relationship between compiler developers and developer-users is toxic as fuck, but as long as we keep comparing compilers by increasingly anal microbenchmarks instead of increasingly anal "correct in practice" tests, I don't see why we continue to find this a surprise. When we demand better performance, compiler developers are going to exploit undefined behavior as well as other language "loopholes" that don't make intuitive sense. Perhaps a new focus on correctness in practice would change this.
- pkaye 10y ago
  What would you suggest as an alternative to float?
  - CalChris 10y ago
    UNUMs. https://en.wikipedia.org/wiki/Unum_(number_format) https://en.wikipedia.org/wiki/Unum_(number_format)
  - dgudkov 10y ago
    128-bit decimal.
- yongjik 10y ago
  There are numbers that are safe to compare against. Storing 1.0 (or 0.0 or 0.25 or any small integer), and then later comparing for equality to check "Did I assign this value before?" is completely fine. Of course, it's not fine to ask "Does the result of this computation equal this value?", but everyone knows that.
  - jobigoud 10y ago
    > Of course, it's not fine to ask "Does the result of this computation equal this value?", but everyone knows that. How about "Does the result of this computation equal the result of this exact same computation?"? That's the issue here.
- fdej 10y ago
  That's not the issue here at all. The IEEE 754 standard does in principle make floating-point operations completely deterministic, i.e. if you do a specific sequence of floating-point operations, then the result may indeed be slightly off from the result you would expect with infinite precision arithmetic, but it should be off in just the same way every time. For example (x + 0.1) + 0.1 and x + 0.2 can give different results due to rounding, but the respective result should be constant for a given x regardless of the machine, compiler, or phase of the moon. Knowing this, there are situations where exact comparisons are appropriate. The "bug" here is that a single operation x + 0.1 and x + 0.1 can give two different results depending on the mood of the compiler. It's not actually a bug, because the IEEE 754 standard permits this kind of sloppiness by default. You can get stricter semantics, but you have to ask for it explicitly by setting the right compiler flags. This is why I say that the IEEE 754 standard makes floating-point arithmetic deterministic "in principle" - in practice, you have to work for it. (Changes to the rounding mode are another big issue here.) An extra error of 1 ulp somewhere might not seem like much for your average numerical code, but it's a big deal in some situations: * Assumptions about the precise behavior of floating-point cancellations are often used in numerical libraries, and deviations can cause disastrous loss of accuracy. For example, Kahan summation is a widely used algorithm for computing sums with increased accuracy. If compilers optimize the code too aggressively, the algorithm doesn't work and you get worse accuracy. And ironically, x87 long doubles can make general numerical calculations less accurate (not more accurate as intended) due to the double rounding issue. * If you use floats in game logic, you might want it to behave exactly the same way on different platforms, e.g. for synchronization or replays. * In scientific computing, you might want results to be precisely reproducible by others. * Debugging!
  - tbirdz 10y ago
    In practice, I don't think IEEE 754 floating point is deterministic, especially not across multiple ISAs, multiple vendors of the same ISA (Intel, AMD), different generations of the same vendor's cpus, different compilers and different math libraries, even just C's libm. So just a warning to everyone out there: unless you are very, very careful your floating point code will not be deterministic. If you thought you can just use floats and expect it to "just work" the same, then you're out of luck. Here's a couple of articles which go into more detail about problems with floating point determinism: http://gafferongames.com/networking-for-game-programmers/floating-point-determinism/ http://gafferongames.com/networking-for-game-programmers/flo... https://randomascii.wordpress.com/2013/07/16/floating-point-determinism/ https://randomascii.wordpress.com/2013/07/16/floating-point-...
- andrepd 10y ago
  Floats are evil? o.O What do you propose instead if you want to deal with non-integer values?
  - rnhmjoj 10y ago
    Continued fractions are perfect if you have infinite memory.
    - kccqzy 10y ago
      A real number can contain an infinite amount of information. If you have infinite memory, you would probably start by asking whether you have countably infinite or uncountably infinite memory.
  - tbirdz 10y ago
    There's always fixed point integers[0]. They're less flexible than floats, but at least the amount of precision available is much more clear and uniform, and all the rounding is pretty explicit. https://en.wikipedia.org/wiki/Fixed-point_arithmetic https://en.wikipedia.org/wiki/Fixed-point_arithmetic https://en.wikipedia.org/wiki/Q_%28number_format%29 https://en.wikipedia.org/wiki/Q_%28number_format%29
    - rwallace 10y ago
      Back in the late eighties, early nineties, I was of the opinion that fixed point is really the way to go, the only thing holding it back was lack of adequate compiler support. Until I realized what would actually happen if the world started using fixed point: programmers would always set the dynamic range too small, so we'd end up with a world full of programs that were constantly failing because they needed to be run on data larger than the test cases. The headaches of floating point, annoying though they can be at times, are much less onerous.
      - MaulingMonkey 10y ago
        > Until I realized what would actually happen if the world started using fixed point: programmers would always set the dynamic range too small, so we'd end up with a world full of programs that were constantly failing because they needed to be run on data larger than the test cases. Floating point doesn't actually fix this issue. A fun bug I dealt with was a black screen caused by a few pixels having NaN values - which were then picked up by a blur pass which would poison basically the entire back buffer. The original NaNs were caused by divides by near-zero values causing out-of-range results. You now add precision issues: Multiple coordinate systems are a common hack-around to maintain local precision when using single-precision based physics engines. Stuttering effects when time is represented as "floating point seconds since uptime" when integer wraparound would avoid any problem in the first place. And this is before even adding e.g. half precision or worse floats into the mix. There's a reason floating point isn't used for a lot of serious banking calculations.
yokohummer7 10y ago
Probably a bug report with the most duplicates? > 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!!! TIL long double is a thing.
- vortico 10y ago
  So would turning off denormals fix the problem? That concept has given me 100 problems in the past and had helped me 0 times.
  - wtallis 10y ago
    Turning off denormals won't fix the issue that a long double has a wider exponent field than a double, and can represent smaller magnitudes without relying on denormals.
  - kr7 10y ago
    The solution is to compile with SSE2 on x86. (flags: -mfpmath=sse -msse -msse2) On x86-64, the compiler should default to SSE2. SSE2 is ~16 years old so compatibility shouldn't be an issue.
Gibbon1 10y ago
https://randomascii.wordpress.com/2013/07/16/floating-point-determinism/ https://randomascii.wordpress.com/2013/07/16/floating-point-...
Animats 10y ago
*DIAGNOSTIC* THE TEST FOR EQUALITY BETWEEN NON-INTEGERS MAY NOT BE MEANINGFUL UNIVAC 1107 FORTRAN IV error message for this, 1967.
Negitivefrags 10y ago
This is the bug that makes std::set< float > crash randomly under GCC with default compile options. The fact that they consider this reasonable is insane.
- emmelaich 10y ago
  Is having a set of floats reasonable?
  - Negitivefrags 10y ago
    Yes. There is nothing about floats that should make them unusable in a set. In particular though, it's fairly common to use std::sets (or std::maps which use the same underlying data structure) to maintain something in a sorted order. I ran into this myself when implementing A*. I was using an std::map with the float as key for the priority queue. Under GCC with default compile arguments it would randomly hang or crash due to GCC breaking the set invariant.
    - Sharlin 10y ago
      > Yes. There is nothing about floats that should make them unusable in a set. Strictly speaking, there is. The < operator for floating-point types does not impose a strict weak order because NaNs are a thing. In any case, comparing floating-point values for strict equality (instead of using an appropriate epsilon) is more often a bad idea than not.
      - temac 10y ago
        Well the thing is, it should not fail because of obscure reasons completely unrelated to NaNs...
      - Someone 10y ago
        But std::set doesn't use operator< for comparing values, but !(a<b) && !(b<a). http://www.cplusplus.com/reference/set/set http://www.cplusplus.com/reference/set/set: "Compare A binary predicate that takes two arguments of the same type as the elements and returns a bool. The expression comp(a,b), where comp is an object of this type and a and b are key values, shall return true if a is considered to go before b in the strict weak ordering the function defines. The set object uses this expression to determine both the order the elements follow in the container and whether two element keys are equivalent (by comparing them reflexively: they are equivalent if !comp(a,b) && !comp(b,a)). No two elements in a set container can be equivalent. This can be a function pointer or a function object (see constructor for an example). This defaults to less<T>, which returns the same as applying the less-than operator (a<b)." I think that guarantees that std::set works for IEEE floating point values, and that NaNs are at the end of the set.
gruntz 10y ago
https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.h...
twic 10y ago
In case anyone is wondering why this showed up right now, it's because tbrowbdidnso was reading the 'examples of unexpected mathematical images' article [1] and came across the one about FPU-dependent aliasing artifacts [2]. [1] https://news.ycombinator.com/item?id=13736568 https://news.ycombinator.com/item?id=13736568 [2] http://www.cgl.uwaterloo.ca/csk/projects/alias/ http://www.cgl.uwaterloo.ca/csk/projects/alias/
- tbrowbdidnso 10y ago
  Eh, was wondering how I wandered to the page myself.
cperciva 10y ago
I met William Kahan about 15 years ago and asked him about this. His response: "If your program says double then yes, of course it's a bug if the compiler does long-double precision arithmetic instead. Isn't that obvious? Why are you even asking this?" If the author of the IEEE 754 standard says that the way you do floating-point arithmetic is obviously a bug, you're probably doing it wrong.
- acqq 10y ago
  And I believe he's right, if he was asked specifically what that program should do. If all the statements in the program are really just: #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); } Then the comparison has to be done on doubles, and "but it's extended precision internally" mentioned in some answers is the wrong answer. All the variables in the program are stated to be doubles. If I'd write the x86 and x87 asm instructions I could easily write them so that they work as expected, even without SSE2. However if you compile with some "special" optimizations that are documented to relax some specific code generation rules, then it can be OK that these specific optimizations behave differently, but then: you asked for it by choosing that documented optimization. In the answers there, Richard Henderson wrote 2001-01-16 06:30:04 UTC: "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." He's also right, if the -ffloat-store is documented to be turned off by -O and if it's documented that that flag also affects the precomputed calculations that are done in the compile time. If the second detail is not documented, the results could be unexpected, but again depends of what you expect from the optimizer. The properly written compiler and optimizer should perform the calculations during the code generation correct to the definition and not cut the corners. If the generated code was generated to actually perform the != during the run time, then the compile-time calculations are irrelevant, but the sane default should be something like -ffloat-store when not using SSE2 even with -O, if -O has the semantic "the optimizations that should not produce really unexpected effects." It sounds complicated, but it has its reasons. P.S. answer to nsajko's post below: I mention -ffloat-store not ffast-math. Read the description of ffloat-store please: https://gcc.gnu.org/onlinedocs/gcc-3.2/gcc/Optimize-Options.html https://gcc.gnu.org/onlinedocs/gcc-3.2/gcc/Optimize-Options.... and also "Note on x86 and m68080 floating-point math" here https://gcc.gnu.org/wiki/FloatingPointMath https://gcc.gnu.org/wiki/FloatingPointMath
solarexplorer 10y ago
Always a good read: David Goldberg: What Every Computer Scientist Should Know About Floating-Point Arithmetic (1991) https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.h...
77pt77 10y ago
Comment 109 has a good exposition I'll just copy here: Jan Lachnitt 2008-05-20 16:59:31 UTC 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 = 2s; } . 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.
shruubi 10y ago
> Regardless of where one wishes to put the blame, this problem will _not_ be fixed. Period. I realise that floating-point numbers are, to use the technical term, "whack", and I also realise that triaging the issue seems to point to a hardware error, not a software error, but I can't get passed this "yeah, there is a problem, so what?" attitude. Am I missing something critical here that makes this unreasonable to work around or fix?
- x0x0 10y ago
  It's not a hardware error at all. If you want to use ieee 784 on a processor, you have to learn how it works; there are some subtleties whose root cause is attempting to approximate infinite precision mathematics in a finite precision processor.
- pizza234 10y ago
  According [indirectly] to one of the comments, the issue is that there is no solution which solves all the problems in the context. Reference: https://gcc.gnu.org/ml/gcc/2003-08/msg01257.html https://gcc.gnu.org/ml/gcc/2003-08/msg01257.html Specifically: There are two things missing. The abililty to turn on the workaround in c (i.e. emit FP rounding instruction after every operation), and a bug fix for the register spills. Even with those things, I think we are still in trouble. In the first case, having explicit rounding instructions eliminates the excess precision problem, but it introduces a double-rounding problem. So we have lost again there. This is probably not fixable without changing the hardware. In the second case, fixing the problem with reload spills eliminates one source of unpredictable rounding, however, we still have the problem that local variables get rounded if they are allocated to the stack and not rounded if they are allocated to an FP reg-stack register. Thus we still have the problem that we get different results at different optimization levels. So we still lose there again also. This might be fixable by promoting all stack locals to long double to avoid unexpected rounding, but that will probably cause other problems in turn. It will break all programs that expect assignments doubles will round to double for instance. If we don't promote stack locals, then we need explicit rounding for them, and then we have double-rounding again. I really see no way to fix this problem other than by fixing the hardware. The hardware has to have explicit float and double operations.
Kenji 10y ago
This is well-known by everyone who works with float determinism. You have to set the FPU flag with assembler to force the desired width (among other things). See: fldcw, LDMXCSR, etc. Take into account rounding modes and denormals too, and make sure to disable float optimisations that break determinism in your compiler.
nayuki 10y ago
This idea is related to the strictfp modifier in Java. https://en.wikipedia.org/wiki/Strictfp https://en.wikipedia.org/wiki/Strictfp
ChuckMcM 10y ago
Hah! That is such a great read, all of the "this isn't a bug" and the tombstones of all the "duplicate" bug reports. Reminds me how challenging it is to do decimal floating point in binary. Mike Cowlishaw could go on for literally hours how just the concept is so broken and he wrote and shipped a decimal based floating point library which has none of these issues[1]. I still remember when I was told to never compare my floating point result to zero because even when it should be it probably won't be, always compare it to be less than a small sigma. [1] http://speleotrove.com/decimal/ http://speleotrove.com/decimal/
- mauvehaus 10y ago
  > always compare it to be less than a small sigma. Unfortunately, this approach can result in problems with the usual expectation that equality is transitive. Consider 0 < x < y < sigma. 0 == x, 0 == y, but x != y. This problem applies generally any time you're doing equality comparisons with sigma.
  - SolarNet 10y ago
    You aren't following the instructions. You have to rewrite all comparisons to use a small sigma. `0.f == x` becomes `SIGMA > x` `0.f == y` becomes `SIGMA > y` `x == y` becomes `SIGMA > abs(x - y)` All the expression work if written correctly.
    - mauvehaus 10y ago
      I probably should have given a better example instead of the hasty one. -sigma < x < 0 < y < sigma, abs(x - y) > sigma. I hope we can agree that this causes problems with transitivity :-) There are also problems with granularity. Suppose x and y are both much, much larger than sigma, and the smallest representable difference between x and y is e.g. 1.0. Now there's no possible way that abs(x - y) < sigma. I suspect that the resolution to this is to compare equality within some number of ULPs, but: 1) You still have the transitivity problem and 2) This is not my field of expertise, and there may be a more appropriate solution.
trelliscoded 10y ago
I remember that in university we had a test on a toy 16-bit floating point ALU. We had to walk through each of the steps in the ALU for the usual operations, showing the bit patterns at each step and keeping track of precision loss, etc. Because we had to track state and some aspects of the test were randomized, there were a few cases where identical operations on identical inputs resulted in different outputs due to the residual state. During the next lecture, someone sticks their hand up and says, "this FPU thing is so nondeterministic, who in their right mind would use floating point for 99% of the kinds of work that we use computers for?" The professor says, "you should think real hard about that the next time you decide to use a floating point type in your program."
- brucedawson 10y ago
  > identical operations on identical inputs resulted > in different outputs due to the residual state. What? What residual state? Floating-point math is tricky, but on a particular machine/compiler/compile-options set you should expect that identical operations on identical inputs should give identical results. What is this "residual state" you mention? If you have an FPU that lets you change precision or rounding modes then you can get different results from the same operations - is that what you're talking about? I get frustrated when people think that floating-point just randomly changes bits for no reason. It doesn't. It randomly changes bits for specific reasons, which can be understood.
  - effie 10y ago
    trelliscoded probably meant bits in scratchpad, which are not accessible to the programmer. I read somewhere that those hidden bits may not be saved and regenerated properly as part of cpu resuming execution of the program (task switching). I do not know if that can happen, but even if not, achieving consistency on subsequent runs is quite difficult - apparently one has to make sure the program and all libraries it uses were compiled with special compiler flags that prevent data alignment influencing results of FPU operations. https://software.intel.com/articles/consistency-of-floating-point-results-using-the-intel-compiler/ https://software.intel.com/articles/consistency-of-floating-...
greglindahl 10y ago
For the PathScale compiler's initial release on x86, in the unlikely event that the x87 was used at all, I set -mpc64 as the default. We never got a customer complaint about it, despite having many numerically-savvy customers.
beached_whale 10y ago
I have found that I often end up forsaking the performance and using a decimal type. It is definitely slower, but the correctness is often more desirable.
- Gibbon1 10y ago
  A friend mentioned being in a meeting where Mike Cowlishaw argued that Javascript's floating point type should be decimal floats not binary floats. Everyone including my friend thought he was out to lunch. My friend says now, 'a great opportunity was pissed away' Because almost no one actually wants wants binary floating point.
  - brucedawson 10y ago
    Why not? Decimal floats have definite advantages if you are dealing with decimal data, but that is not the only type of data that floats are used for, and it probably isn't the majority. Binary floating point is better than decimal floating point if you don't care about decimal-interop (more precision per bit, more consistent precision), and decimal floating point has all the same issues (except for confusing conversions to decimals). So, if you are dealing with raw sensor data, neural nets, character models in a game, numerical simulations, and many other types of data, the extra precision of binary floating-point is far more important than any confusion that arises during those (rare) conversions to decimal.
    - beached_whale 10y ago
      It's all about costs and benefits though. I do think that a decimal type should be standard along side a float so that the most appropriate one can be chosen. But in many contexts, one needs to go outside the immediate ecosystem to third party libraries to get an alternative to floating point. Currency is a good example of a type that should not use floating point arithmetic but is far too often.
      - Gibbon1 10y ago
        I think the difference is binary floating point tends to be on massive automatically generated sets of data. And binary floats are mostly perfect for that. However the bias towards scientific computing is a problem. because there are large numbers of use cases where the inputs and outputs are decimal floats because you're dealing with humans and their preference for base 10.
phkahler 10y ago
After reading much of the discussion, my conclusion is that it's past time to deprecate the x87 FPU. Ever since x86_64, it is expected that SSE2 will be available and compilers are expected to use it over the old FPU. Face it, this is an old Intel hardware "feature" that causes problems and has been obsolete for a long time.
hzhou321 10y ago
People need understand that floating point comparison at floating precision is undefined, so it is not a bug. The proper floating number representation should track its precision. Since that is more complicated and potentially performance affecting, so until then the programmer should manage the precision tracking on their own (similar to memory allocations), that is, one should do: if ( f_a < f_b-f_eps) if ( f_a > f_b+f_eps) if (fabs(f_a-f_b)<f_eps) If you don't do this, then you are on your own. Sometimes it is fine, but sometimes it may give you grief, just like the dangers of other undefined behaviors. PS: I also have a problem that people think "undefined behavior" is undefined because certain standards say so. Undefined behavior is due to lack of agreement of common sense or defies common sense. For example, divide by zero, there is no common sense doing so, and define its behavior does not change anything. Signed integer overflow, there is no common sense agreement for that behavior either. The common sense for undefined behaviors is to avoid it. (On the other hand, there seems to be some common sense agreement on unsigned integer overflow, and people are writing code based on that behavior. But that agreement is weak, so I would avoid it if I can regardless what the standards say.) PS2: The reason that floating point comparison is undefined behavior because, unlike fixed precision or integers, floating point number has floating precision. Floating precision does not have common sense correspondence. In any real application (where common sense is based), a finite, and above all, consistent precision, is always assumed. And the solution, by manually specifying a f_eps, we put that comparison into a common-sense understanding that you know the result of comparisons is only meaningful above f_eps. With that common sense, having an internal calculation in higher precision is always OK -- isn't that common sense? PS3: The question is why not turn undefined behaviors into errors? Undefined behavior does not mean there are no behaviors, and the behaviors can be used, or more often, ignored to achieve certain efficiency. When programmers write code with undefined behaviors, they are extracting their code out of common sense domain, but into a special sense domain. Special sense domain is practically fine. Everyone has their own efficient quirks, just do not expect your quirks to be understood or supported by out large where common sense rules (unless you carefully define it and document it). A language that allows for undefined behaviors allows for individual efficiency, which sometimes it is rather good. For example, it is ok to write programs with direct floating point number comparisons as long as you understand that your program does not care when the comparison falls into the ambiguous range. Only you who understand your special application can assess that. If you are the only one who is responsible for the result of that code, then it is perfectly fine!
nortiero 10y ago
I think that this is a case of mismatch between two standards, not a bug. The C standard allows higher precision values to be used in place of lesser ones (extended precision vs double) and REQUIRES conversion to the correct precision during assignment (or casting). Also, the C abstract machine has the option to store those values when/if it deems opportune, it won't change observables as defined there. On the other side, IEEE 754 allows extended precision to be used in place of double and of course requires that any chosen precision be kept or else. But C Standard mandates IEEE 754 , too! It seems to me that modern C deliberately chose to ignore such kind of mismatch in the name of (substantial) performance gains. K&R, good or bad, was way simpler!
- leni536 10y ago
  > But C Standard mandates IEEE 754 , too! Are you sure? I can't find it in the standard, specifically not in [1] page 28., section 5.2.4.2.2. [1] (the final draft of C1X, dated 12 April 2011) http://www.open-std.org/jtc1/sc22/wg14/www/docs/n1570.pdf http://www.open-std.org/jtc1/sc22/wg14/www/docs/n1570.pdf
  - nortiero 10y ago
    You are absolutely correct, "mandates" is too strong a word. Annex F, which is normative, have a way out by not setting __STDC_IEC_559__, GCC6.3.0 does set a slightly different GCC_IEC_559. Yet I think that my argument still holds: floating point calculations can be executed at wider range (§5.1.2.3.12); assignments and casts are under obligation to wipe out that extra precision (§ same paragraph); I am no language lawyer, but .. given the issue of program observable effects (§5.1.2.3.12) and the "as is" implied rule that governs optimizations, how possibly could equality be stable over time?
    - leni536 10y ago
      Yeah, you are right. It seems to be that the problem is that we cram two standards into too few types. Additionally to float, double, and long double there could be optional iec559_single_t, iec559_double_t and iec559_extended_t (similarly to stdint.h integer types). The operations on variables of these types could be strictly iec559 conformant.
caf 10y ago
The -fexcess-precision=standard option, which is enabled by -std=c99 (requesting ICO C99 compliance), forces excess precision to be rounded away on all assignments and casts to narrower precision types. See https://gcc.gnu.org/onlinedocs/gcc/Optimize-Options.html#index-fexcess-precision-960 https://gcc.gnu.org/onlinedocs/gcc/Optimize-Options.html#ind...
Loic 10y ago
Because I am doing a lot of numerical code in Fortran using the Intel compiler on Windows for production and GCC/GNU Fortran for prototyping on Linux I was a bit scared by the title. This is a problem with the equality comparison with doubles, something I try to avoid as much as possible, as this was the first thing I learnt to do 25 years ago... never compare double for equality. Anyway, the bug report is linking to a very nice research paper[0]: "The pitfalls of verifying floating-point computations", which even so very technical is a good read on the subject. [0]: https://hal.archives-ouvertes.fr/hal-00128124 https://hal.archives-ouvertes.fr/hal-00128124
vanderZwan 10y ago
I'm very surprised nobody has brought up John Gustafson's Unum work. His latest iteration is the posit, which he presented at Stanford earlier this month[0][1] EDIT: apologies to CalChris[2]. What I should have said is discussed Gustafson's Unum work. [0] https://www.youtube.com/watch?v=aP0Y1uAA-2Y https://www.youtube.com/watch?v=aP0Y1uAA-2Y (relevant question at 1h22m37s) [1] http://web.stanford.edu/class/ee380/Abstracts/170201.html http://web.stanford.edu/class/ee380/Abstracts/170201.html [2] https://news.ycombinator.com/item?id=13739617 https://news.ycombinator.com/item?id=13739617
jlgustafson 10y ago
I have been following this bug in its various form for 30 years, and it is covered in pages 55-62 of The End of Error: Unum Computing. It is the "hidden scratchpad" bug where designers of a computing environment try to be helpful by providing more precision than what is explicitly requested by the programmer. If you read the blog postings of William Kahan about this issue, you will find he speaks fervently on both sides, sometimes defending it and other times vehemently deriding it. Thus the confusion. Kahan favors any technique that improves accuracy, even if it creates behavior that can only be explained through careful numerical analysis (which very few programmers have the time or skill set to perform). Unums and posits, incidentally, are designed to eliminate all sources of irreproducibility and inexplicable behavior like this gcc bug. Hardware support for posits is on the way, maybe before the year is out. There is a way out of this mess, but it will take a few years to become mainstream.
- wbl 10y ago
  No, Kahan advocates for what the C standard does which GCC doesn't implement correctly.
gasgiant 10y ago
As a system administrator, my favorite part of the Bug 323 thread is the mail bounces from postmaster.
laurent123456 10y ago
Should add the date to the title (2000)

145 points by tbrowbdidnso 10y ago

castratikron 10y ago

More proof that floats are evil. Moral of the story is to never use equality to compare floats.

trendia 10y ago
Pretty much. If you want to compare two floats a and b within some precision eps, use: abs(a - b) <= eps*0.5*(a + b) or even (abs(a - b) <= eps*abs(a)) && (abs(a - b) <= eps*abs(b)) which is implemented in Boost [0]. [0] http://www.boost.org/doc/libs/1_34_0/libs/test/doc/components/test_tools/floating_point_comparison.html http://www.boost.org/doc/libs/1_34_0/libs/test/doc/component...
- ratboy666 10y ago
  I thought it was (a / b) != (a / b)
geezerjay 10y ago
> More proof that floats are evil. Just because someone doesn't know the basics of floating point numbers it doesn't mean that the floating point format is somehow evil. Testing floats for equality is one of those no-nos which is taught right in the first class of introduction to FP numbers. Another very basic issue covered right in the intro clsses are problems caused by using decimal values to define numbers encoded in binary format, which is often explained through the famous "0.1 doesn't exist in binary" trivia bit.
- thechao 10y ago
  Part of my day job is to test floating point hardware and write compilers and assemblers. Testing equality---with some care---is completely fine. I can assure you this is a bug in GCC. The GCC bug is compounded by the nonstandard behavior of x87. The x87 unit is not "fixable" but GCC's behavior is.
  - nurettin 10y ago
    That's interesting, you make it sound like the production of new FPUs and new compilers is common. Which industry do you work in?
nom 10y ago
Floats are not evil. They are just not as intuitive as integers because it is impossible to represent every possible rational number in the range with a fixed number of bits. Once you get that, everything makes sense.
- monochromatic 10y ago
  That does not make sense of getting a failure from assert x/y == x/y
  - mikeash 10y ago
    Properly implemented floats will never fail that (barring NaN). That gcc sometimes does is a compiler bug, not a float problem.
  - cellularmitosis 10y ago
    I agree. Explaining that assertion requires explaining that the x87 FPU is providing a leaky abstraction. If you ask for double math, it gives you a result which is a long double while it resides in a register, and then becomes a double when moved to RAM. There are some cases where direct comparison of floating point values is what you want, and this issue can still bite you even then.
    - int_19h 10y ago
      For x87, it's not a leaky abstraction - it fully supports 80-bit floats as a type, and documents that its registers are 80-bit, and so are operations on them. There's nothing abstract about it. The bug is that GCC is utilizing the (documented, by-design) x87 behavior in a way that causes it to produce the wrong result, for the sake of performance.
      - abstractbeliefs 10y ago
        >compiler developers do unintuitive but technically legal thing in the name of performance I still think the relationship between compiler developers and developer-users is toxic as fuck, but as long as we keep comparing compilers by increasingly anal microbenchmarks instead of increasingly anal "correct in practice" tests, I don't see why we continue to find this a surprise. When we demand better performance, compiler developers are going to exploit undefined behavior as well as other language "loopholes" that don't make intuitive sense. Perhaps a new focus on correctness in practice would change this.
pkaye 10y ago
What would you suggest as an alternative to float?
- CalChris 10y ago
  UNUMs. https://en.wikipedia.org/wiki/Unum_(number_format) https://en.wikipedia.org/wiki/Unum_(number_format)
- dgudkov 10y ago
  128-bit decimal.
yongjik 10y ago
There are numbers that are safe to compare against. Storing 1.0 (or 0.0 or 0.25 or any small integer), and then later comparing for equality to check "Did I assign this value before?" is completely fine. Of course, it's not fine to ask "Does the result of this computation equal this value?", but everyone knows that.
- jobigoud 10y ago
  > Of course, it's not fine to ask "Does the result of this computation equal this value?", but everyone knows that. How about "Does the result of this computation equal the result of this exact same computation?"? That's the issue here.
fdej 10y ago
That's not the issue here at all. The IEEE 754 standard does in principle make floating-point operations completely deterministic, i.e. if you do a specific sequence of floating-point operations, then the result may indeed be slightly off from the result you would expect with infinite precision arithmetic, but it should be off in just the same way every time. For example (x + 0.1) + 0.1 and x + 0.2 can give different results due to rounding, but the respective result should be constant for a given x regardless of the machine, compiler, or phase of the moon. Knowing this, there are situations where exact comparisons are appropriate. The "bug" here is that a single operation x + 0.1 and x + 0.1 can give two different results depending on the mood of the compiler. It's not actually a bug, because the IEEE 754 standard permits this kind of sloppiness by default. You can get stricter semantics, but you have to ask for it explicitly by setting the right compiler flags. This is why I say that the IEEE 754 standard makes floating-point arithmetic deterministic "in principle" - in practice, you have to work for it. (Changes to the rounding mode are another big issue here.) An extra error of 1 ulp somewhere might not seem like much for your average numerical code, but it's a big deal in some situations: * Assumptions about the precise behavior of floating-point cancellations are often used in numerical libraries, and deviations can cause disastrous loss of accuracy. For example, Kahan summation is a widely used algorithm for computing sums with increased accuracy. If compilers optimize the code too aggressively, the algorithm doesn't work and you get worse accuracy. And ironically, x87 long doubles can make general numerical calculations less accurate (not more accurate as intended) due to the double rounding issue. * If you use floats in game logic, you might want it to behave exactly the same way on different platforms, e.g. for synchronization or replays. * In scientific computing, you might want results to be precisely reproducible by others. * Debugging!
- tbirdz 10y ago
  In practice, I don't think IEEE 754 floating point is deterministic, especially not across multiple ISAs, multiple vendors of the same ISA (Intel, AMD), different generations of the same vendor's cpus, different compilers and different math libraries, even just C's libm. So just a warning to everyone out there: unless you are very, very careful your floating point code will not be deterministic. If you thought you can just use floats and expect it to "just work" the same, then you're out of luck. Here's a couple of articles which go into more detail about problems with floating point determinism: http://gafferongames.com/networking-for-game-programmers/floating-point-determinism/ http://gafferongames.com/networking-for-game-programmers/flo... https://randomascii.wordpress.com/2013/07/16/floating-point-determinism/ https://randomascii.wordpress.com/2013/07/16/floating-point-...
andrepd 10y ago
Floats are evil? o.O What do you propose instead if you want to deal with non-integer values?
- rnhmjoj 10y ago
  Continued fractions are perfect if you have infinite memory.
  - kccqzy 10y ago
    A real number can contain an infinite amount of information. If you have infinite memory, you would probably start by asking whether you have countably infinite or uncountably infinite memory.
- tbirdz 10y ago
  There's always fixed point integers[0]. They're less flexible than floats, but at least the amount of precision available is much more clear and uniform, and all the rounding is pretty explicit. https://en.wikipedia.org/wiki/Fixed-point_arithmetic https://en.wikipedia.org/wiki/Fixed-point_arithmetic https://en.wikipedia.org/wiki/Q_%28number_format%29 https://en.wikipedia.org/wiki/Q_%28number_format%29
  - rwallace 10y ago
    Back in the late eighties, early nineties, I was of the opinion that fixed point is really the way to go, the only thing holding it back was lack of adequate compiler support. Until I realized what would actually happen if the world started using fixed point: programmers would always set the dynamic range too small, so we'd end up with a world full of programs that were constantly failing because they needed to be run on data larger than the test cases. The headaches of floating point, annoying though they can be at times, are much less onerous.
    - MaulingMonkey 10y ago
      > Until I realized what would actually happen if the world started using fixed point: programmers would always set the dynamic range too small, so we'd end up with a world full of programs that were constantly failing because they needed to be run on data larger than the test cases. Floating point doesn't actually fix this issue. A fun bug I dealt with was a black screen caused by a few pixels having NaN values - which were then picked up by a blur pass which would poison basically the entire back buffer. The original NaNs were caused by divides by near-zero values causing out-of-range results. You now add precision issues: Multiple coordinate systems are a common hack-around to maintain local precision when using single-precision based physics engines. Stuttering effects when time is represented as "floating point seconds since uptime" when integer wraparound would avoid any problem in the first place. And this is before even adding e.g. half precision or worse floats into the mix. There's a reason floating point isn't used for a lot of serious banking calculations.

yokohummer7 10y ago

Probably a bug report with the most duplicates? > 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!!! TIL long double is a thing.

vortico 10y ago
So would turning off denormals fix the problem? That concept has given me 100 problems in the past and had helped me 0 times.
- wtallis 10y ago
  Turning off denormals won't fix the issue that a long double has a wider exponent field than a double, and can represent smaller magnitudes without relying on denormals.
- kr7 10y ago
  The solution is to compile with SSE2 on x86. (flags: -mfpmath=sse -msse -msse2) On x86-64, the compiler should default to SSE2. SSE2 is ~16 years old so compatibility shouldn't be an issue.

Gibbon1 10y ago

https://randomascii.wordpress.com/2013/07/16/floating-point-determinism/ https://randomascii.wordpress.com/2013/07/16/floating-point-...

Animats 10y ago

*DIAGNOSTIC* THE TEST FOR EQUALITY BETWEEN NON-INTEGERS MAY NOT BE MEANINGFUL UNIVAC 1107 FORTRAN IV error message for this, 1967.

Negitivefrags 10y ago

This is the bug that makes std::set< float > crash randomly under GCC with default compile options. The fact that they consider this reasonable is insane.

emmelaich 10y ago
Is having a set of floats reasonable?
- Negitivefrags 10y ago
  Yes. There is nothing about floats that should make them unusable in a set. In particular though, it's fairly common to use std::sets (or std::maps which use the same underlying data structure) to maintain something in a sorted order. I ran into this myself when implementing A*. I was using an std::map with the float as key for the priority queue. Under GCC with default compile arguments it would randomly hang or crash due to GCC breaking the set invariant.
  - Sharlin 10y ago
    > Yes. There is nothing about floats that should make them unusable in a set. Strictly speaking, there is. The < operator for floating-point types does not impose a strict weak order because NaNs are a thing. In any case, comparing floating-point values for strict equality (instead of using an appropriate epsilon) is more often a bad idea than not.
    - temac 10y ago
      Well the thing is, it should not fail because of obscure reasons completely unrelated to NaNs...
    - Someone 10y ago
      But std::set doesn't use operator< for comparing values, but !(a<b) && !(b<a). http://www.cplusplus.com/reference/set/set http://www.cplusplus.com/reference/set/set: "Compare A binary predicate that takes two arguments of the same type as the elements and returns a bool. The expression comp(a,b), where comp is an object of this type and a and b are key values, shall return true if a is considered to go before b in the strict weak ordering the function defines. The set object uses this expression to determine both the order the elements follow in the container and whether two element keys are equivalent (by comparing them reflexively: they are equivalent if !comp(a,b) && !comp(b,a)). No two elements in a set container can be equivalent. This can be a function pointer or a function object (see constructor for an example). This defaults to less<T>, which returns the same as applying the less-than operator (a<b)." I think that guarantees that std::set works for IEEE floating point values, and that NaNs are at the end of the set.

gruntz 10y ago

https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.h...

twic 10y ago

In case anyone is wondering why this showed up right now, it's because tbrowbdidnso was reading the 'examples of unexpected mathematical images' article [1] and came across the one about FPU-dependent aliasing artifacts [2]. [1] https://news.ycombinator.com/item?id=13736568 https://news.ycombinator.com/item?id=13736568 [2] http://www.cgl.uwaterloo.ca/csk/projects/alias/ http://www.cgl.uwaterloo.ca/csk/projects/alias/

tbrowbdidnso 10y ago
Eh, was wondering how I wandered to the page myself.

cperciva 10y ago

I met William Kahan about 15 years ago and asked him about this. His response: "If your program says double then yes, of course it's a bug if the compiler does long-double precision arithmetic instead. Isn't that obvious? Why are you even asking this?" If the author of the IEEE 754 standard says that the way you do floating-point arithmetic is obviously a bug, you're probably doing it wrong.

acqq 10y ago
And I believe he's right, if he was asked specifically what that program should do. If all the statements in the program are really just: #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); } Then the comparison has to be done on doubles, and "but it's extended precision internally" mentioned in some answers is the wrong answer. All the variables in the program are stated to be doubles. If I'd write the x86 and x87 asm instructions I could easily write them so that they work as expected, even without SSE2. However if you compile with some "special" optimizations that are documented to relax some specific code generation rules, then it can be OK that these specific optimizations behave differently, but then: you asked for it by choosing that documented optimization. In the answers there, Richard Henderson wrote 2001-01-16 06:30:04 UTC: "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." He's also right, if the -ffloat-store is documented to be turned off by -O and if it's documented that that flag also affects the precomputed calculations that are done in the compile time. If the second detail is not documented, the results could be unexpected, but again depends of what you expect from the optimizer. The properly written compiler and optimizer should perform the calculations during the code generation correct to the definition and not cut the corners. If the generated code was generated to actually perform the != during the run time, then the compile-time calculations are irrelevant, but the sane default should be something like -ffloat-store when not using SSE2 even with -O, if -O has the semantic "the optimizations that should not produce really unexpected effects." It sounds complicated, but it has its reasons. P.S. answer to nsajko's post below: I mention -ffloat-store not ffast-math. Read the description of ffloat-store please: https://gcc.gnu.org/onlinedocs/gcc-3.2/gcc/Optimize-Options.html https://gcc.gnu.org/onlinedocs/gcc-3.2/gcc/Optimize-Options.... and also "Note on x86 and m68080 floating-point math" here https://gcc.gnu.org/wiki/FloatingPointMath https://gcc.gnu.org/wiki/FloatingPointMath

solarexplorer 10y ago

Always a good read: David Goldberg: What Every Computer Scientist Should Know About Floating-Point Arithmetic (1991) https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.h...

77pt77 10y ago

Comment 109 has a good exposition I'll just copy here: Jan Lachnitt 2008-05-20 16:59:31 UTC 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 = 2s; } . 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.

shruubi 10y ago

> Regardless of where one wishes to put the blame, this problem will _not_ be fixed. Period. I realise that floating-point numbers are, to use the technical term, "whack", and I also realise that triaging the issue seems to point to a hardware error, not a software error, but I can't get passed this "yeah, there is a problem, so what?" attitude. Am I missing something critical here that makes this unreasonable to work around or fix?

x0x0 10y ago
It's not a hardware error at all. If you want to use ieee 784 on a processor, you have to learn how it works; there are some subtleties whose root cause is attempting to approximate infinite precision mathematics in a finite precision processor.
pizza234 10y ago
According [indirectly] to one of the comments, the issue is that there is no solution which solves all the problems in the context. Reference: https://gcc.gnu.org/ml/gcc/2003-08/msg01257.html https://gcc.gnu.org/ml/gcc/2003-08/msg01257.html Specifically: There are two things missing. The abililty to turn on the workaround in c (i.e. emit FP rounding instruction after every operation), and a bug fix for the register spills. Even with those things, I think we are still in trouble. In the first case, having explicit rounding instructions eliminates the excess precision problem, but it introduces a double-rounding problem. So we have lost again there. This is probably not fixable without changing the hardware. In the second case, fixing the problem with reload spills eliminates one source of unpredictable rounding, however, we still have the problem that local variables get rounded if they are allocated to the stack and not rounded if they are allocated to an FP reg-stack register. Thus we still have the problem that we get different results at different optimization levels. So we still lose there again also. This might be fixable by promoting all stack locals to long double to avoid unexpected rounding, but that will probably cause other problems in turn. It will break all programs that expect assignments doubles will round to double for instance. If we don't promote stack locals, then we need explicit rounding for them, and then we have double-rounding again. I really see no way to fix this problem other than by fixing the hardware. The hardware has to have explicit float and double operations.

Kenji 10y ago

This is well-known by everyone who works with float determinism. You have to set the FPU flag with assembler to force the desired width (among other things). See: fldcw, LDMXCSR, etc. Take into account rounding modes and denormals too, and make sure to disable float optimisations that break determinism in your compiler.

nayuki 10y ago

This idea is related to the strictfp modifier in Java. https://en.wikipedia.org/wiki/Strictfp https://en.wikipedia.org/wiki/Strictfp

ChuckMcM 10y ago

Hah! That is such a great read, all of the "this isn't a bug" and the tombstones of all the "duplicate" bug reports. Reminds me how challenging it is to do decimal floating point in binary. Mike Cowlishaw could go on for literally hours how just the concept is so broken and he wrote and shipped a decimal based floating point library which has none of these issues[1]. I still remember when I was told to never compare my floating point result to zero because even when it should be it probably won't be, always compare it to be less than a small sigma. [1] http://speleotrove.com/decimal/ http://speleotrove.com/decimal/

mauvehaus 10y ago
> always compare it to be less than a small sigma. Unfortunately, this approach can result in problems with the usual expectation that equality is transitive. Consider 0 < x < y < sigma. 0 == x, 0 == y, but x != y. This problem applies generally any time you're doing equality comparisons with sigma.
- SolarNet 10y ago
  You aren't following the instructions. You have to rewrite all comparisons to use a small sigma. `0.f == x` becomes `SIGMA > x` `0.f == y` becomes `SIGMA > y` `x == y` becomes `SIGMA > abs(x - y)` All the expression work if written correctly.
  - mauvehaus 10y ago
    I probably should have given a better example instead of the hasty one. -sigma < x < 0 < y < sigma, abs(x - y) > sigma. I hope we can agree that this causes problems with transitivity :-) There are also problems with granularity. Suppose x and y are both much, much larger than sigma, and the smallest representable difference between x and y is e.g. 1.0. Now there's no possible way that abs(x - y) < sigma. I suspect that the resolution to this is to compare equality within some number of ULPs, but: 1) You still have the transitivity problem and 2) This is not my field of expertise, and there may be a more appropriate solution.

trelliscoded 10y ago

I remember that in university we had a test on a toy 16-bit floating point ALU. We had to walk through each of the steps in the ALU for the usual operations, showing the bit patterns at each step and keeping track of precision loss, etc. Because we had to track state and some aspects of the test were randomized, there were a few cases where identical operations on identical inputs resulted in different outputs due to the residual state. During the next lecture, someone sticks their hand up and says, "this FPU thing is so nondeterministic, who in their right mind would use floating point for 99% of the kinds of work that we use computers for?" The professor says, "you should think real hard about that the next time you decide to use a floating point type in your program."

brucedawson 10y ago
> identical operations on identical inputs resulted > in different outputs due to the residual state. What? What residual state? Floating-point math is tricky, but on a particular machine/compiler/compile-options set you should expect that identical operations on identical inputs should give identical results. What is this "residual state" you mention? If you have an FPU that lets you change precision or rounding modes then you can get different results from the same operations - is that what you're talking about? I get frustrated when people think that floating-point just randomly changes bits for no reason. It doesn't. It randomly changes bits for specific reasons, which can be understood.
- effie 10y ago
  trelliscoded probably meant bits in scratchpad, which are not accessible to the programmer. I read somewhere that those hidden bits may not be saved and regenerated properly as part of cpu resuming execution of the program (task switching). I do not know if that can happen, but even if not, achieving consistency on subsequent runs is quite difficult - apparently one has to make sure the program and all libraries it uses were compiled with special compiler flags that prevent data alignment influencing results of FPU operations. https://software.intel.com/articles/consistency-of-floating-point-results-using-the-intel-compiler/ https://software.intel.com/articles/consistency-of-floating-...

greglindahl 10y ago

For the PathScale compiler's initial release on x86, in the unlikely event that the x87 was used at all, I set -mpc64 as the default. We never got a customer complaint about it, despite having many numerically-savvy customers.

beached_whale 10y ago

I have found that I often end up forsaking the performance and using a decimal type. It is definitely slower, but the correctness is often more desirable.

Gibbon1 10y ago
A friend mentioned being in a meeting where Mike Cowlishaw argued that Javascript's floating point type should be decimal floats not binary floats. Everyone including my friend thought he was out to lunch. My friend says now, 'a great opportunity was pissed away' Because almost no one actually wants wants binary floating point.
- brucedawson 10y ago
  Why not? Decimal floats have definite advantages if you are dealing with decimal data, but that is not the only type of data that floats are used for, and it probably isn't the majority. Binary floating point is better than decimal floating point if you don't care about decimal-interop (more precision per bit, more consistent precision), and decimal floating point has all the same issues (except for confusing conversions to decimals). So, if you are dealing with raw sensor data, neural nets, character models in a game, numerical simulations, and many other types of data, the extra precision of binary floating-point is far more important than any confusion that arises during those (rare) conversions to decimal.
  - beached_whale 10y ago
    It's all about costs and benefits though. I do think that a decimal type should be standard along side a float so that the most appropriate one can be chosen. But in many contexts, one needs to go outside the immediate ecosystem to third party libraries to get an alternative to floating point. Currency is a good example of a type that should not use floating point arithmetic but is far too often.
    - Gibbon1 10y ago
      I think the difference is binary floating point tends to be on massive automatically generated sets of data. And binary floats are mostly perfect for that. However the bias towards scientific computing is a problem. because there are large numbers of use cases where the inputs and outputs are decimal floats because you're dealing with humans and their preference for base 10.

phkahler 10y ago

After reading much of the discussion, my conclusion is that it's past time to deprecate the x87 FPU. Ever since x86_64, it is expected that SSE2 will be available and compilers are expected to use it over the old FPU. Face it, this is an old Intel hardware "feature" that causes problems and has been obsolete for a long time.

hzhou321 10y ago

People need understand that floating point comparison at floating precision is undefined, so it is not a bug. The proper floating number representation should track its precision. Since that is more complicated and potentially performance affecting, so until then the programmer should manage the precision tracking on their own (similar to memory allocations), that is, one should do: if ( f_a < f_b-f_eps) if ( f_a > f_b+f_eps) if (fabs(f_a-f_b)<f_eps) If you don't do this, then you are on your own. Sometimes it is fine, but sometimes it may give you grief, just like the dangers of other undefined behaviors. PS: I also have a problem that people think "undefined behavior" is undefined because certain standards say so. Undefined behavior is due to lack of agreement of common sense or defies common sense. For example, divide by zero, there is no common sense doing so, and define its behavior does not change anything. Signed integer overflow, there is no common sense agreement for that behavior either. The common sense for undefined behaviors is to avoid it. (On the other hand, there seems to be some common sense agreement on unsigned integer overflow, and people are writing code based on that behavior. But that agreement is weak, so I would avoid it if I can regardless what the standards say.) PS2: The reason that floating point comparison is undefined behavior because, unlike fixed precision or integers, floating point number has floating precision. Floating precision does not have common sense correspondence. In any real application (where common sense is based), a finite, and above all, consistent precision, is always assumed. And the solution, by manually specifying a f_eps, we put that comparison into a common-sense understanding that you know the result of comparisons is only meaningful above f_eps. With that common sense, having an internal calculation in higher precision is always OK -- isn't that common sense? PS3: The question is why not turn undefined behaviors into errors? Undefined behavior does not mean there are no behaviors, and the behaviors can be used, or more often, ignored to achieve certain efficiency. When programmers write code with undefined behaviors, they are extracting their code out of common sense domain, but into a special sense domain. Special sense domain is practically fine. Everyone has their own efficient quirks, just do not expect your quirks to be understood or supported by out large where common sense rules (unless you carefully define it and document it). A language that allows for undefined behaviors allows for individual efficiency, which sometimes it is rather good. For example, it is ok to write programs with direct floating point number comparisons as long as you understand that your program does not care when the comparison falls into the ambiguous range. Only you who understand your special application can assess that. If you are the only one who is responsible for the result of that code, then it is perfectly fine!

nortiero 10y ago

I think that this is a case of mismatch between two standards, not a bug. The C standard allows higher precision values to be used in place of lesser ones (extended precision vs double) and REQUIRES conversion to the correct precision during assignment (or casting). Also, the C abstract machine has the option to store those values when/if it deems opportune, it won't change observables as defined there. On the other side, IEEE 754 allows extended precision to be used in place of double and of course requires that any chosen precision be kept or else. But C Standard mandates IEEE 754 , too! It seems to me that modern C deliberately chose to ignore such kind of mismatch in the name of (substantial) performance gains. K&R, good or bad, was way simpler!

leni536 10y ago
> But C Standard mandates IEEE 754 , too! Are you sure? I can't find it in the standard, specifically not in [1] page 28., section 5.2.4.2.2. [1] (the final draft of C1X, dated 12 April 2011) http://www.open-std.org/jtc1/sc22/wg14/www/docs/n1570.pdf http://www.open-std.org/jtc1/sc22/wg14/www/docs/n1570.pdf
- nortiero 10y ago
  You are absolutely correct, "mandates" is too strong a word. Annex F, which is normative, have a way out by not setting __STDC_IEC_559__, GCC6.3.0 does set a slightly different GCC_IEC_559. Yet I think that my argument still holds: floating point calculations can be executed at wider range (§5.1.2.3.12); assignments and casts are under obligation to wipe out that extra precision (§ same paragraph); I am no language lawyer, but .. given the issue of program observable effects (§5.1.2.3.12) and the "as is" implied rule that governs optimizations, how possibly could equality be stable over time?
  - leni536 10y ago
    Yeah, you are right. It seems to be that the problem is that we cram two standards into too few types. Additionally to float, double, and long double there could be optional iec559_single_t, iec559_double_t and iec559_extended_t (similarly to stdint.h integer types). The operations on variables of these types could be strictly iec559 conformant.

caf 10y ago

The -fexcess-precision=standard option, which is enabled by -std=c99 (requesting ICO C99 compliance), forces excess precision to be rounded away on all assignments and casts to narrower precision types. See https://gcc.gnu.org/onlinedocs/gcc/Optimize-Options.html#index-fexcess-precision-960 https://gcc.gnu.org/onlinedocs/gcc/Optimize-Options.html#ind...

Loic 10y ago

Because I am doing a lot of numerical code in Fortran using the Intel compiler on Windows for production and GCC/GNU Fortran for prototyping on Linux I was a bit scared by the title. This is a problem with the equality comparison with doubles, something I try to avoid as much as possible, as this was the first thing I learnt to do 25 years ago... never compare double for equality. Anyway, the bug report is linking to a very nice research paper[0]: "The pitfalls of verifying floating-point computations", which even so very technical is a good read on the subject. [0]: https://hal.archives-ouvertes.fr/hal-00128124 https://hal.archives-ouvertes.fr/hal-00128124

vanderZwan 10y ago

I'm very surprised nobody has brought up John Gustafson's Unum work. His latest iteration is the posit, which he presented at Stanford earlier this month[0][1] EDIT: apologies to CalChris[2]. What I should have said is discussed Gustafson's Unum work. [0] https://www.youtube.com/watch?v=aP0Y1uAA-2Y https://www.youtube.com/watch?v=aP0Y1uAA-2Y (relevant question at 1h22m37s) [1] http://web.stanford.edu/class/ee380/Abstracts/170201.html http://web.stanford.edu/class/ee380/Abstracts/170201.html [2] https://news.ycombinator.com/item?id=13739617 https://news.ycombinator.com/item?id=13739617

jlgustafson 10y ago

I have been following this bug in its various form for 30 years, and it is covered in pages 55-62 of The End of Error: Unum Computing. It is the "hidden scratchpad" bug where designers of a computing environment try to be helpful by providing more precision than what is explicitly requested by the programmer. If you read the blog postings of William Kahan about this issue, you will find he speaks fervently on both sides, sometimes defending it and other times vehemently deriding it. Thus the confusion. Kahan favors any technique that improves accuracy, even if it creates behavior that can only be explained through careful numerical analysis (which very few programmers have the time or skill set to perform). Unums and posits, incidentally, are designed to eliminate all sources of irreproducibility and inexplicable behavior like this gcc bug. Hardware support for posits is on the way, maybe before the year is out. There is a way out of this mess, but it will take a few years to become mainstream.

wbl 10y ago
No, Kahan advocates for what the C standard does which GCC doesn't implement correctly.

gasgiant 10y ago

As a system administrator, my favorite part of the Bug 323 thread is the mail bounces from postmaster.

laurent123456 10y ago

Should add the date to the title (2000)