10 ms·
GNU GCC does not round floating-point divisions to the nearest value
- wiml 6y agoThis seems like a property of the Intel 8087 math coprocessor and its distant descendants found in x86 chips, not anything to do with gcc?
- bhouston 6y agoBut then why does python and JavaScript on the same machines get it right?
- nitrogen 6y agoUnsure: compiled to use SSE instead of x87 for fpmath?
- stefan_ 6y agoBecause it's now getting exceedingly rare to find a program that was not compiled with at least SSE instructions enabled. So the x87 FPU is simply not used then.
- nwallin 6y agoBecause they're compiled to use SSE instructions instead of x87 FPU instructions. If you are still writing 32 bit x86 C code, it's probably some embedded or legacy CPU, and it's wrong to assume SSE instructions are available. So GCC correctly defaults to using FPU instructions instead of SSE instructions on 32 bit x86. Javascript and python will basically never be run in embedded environments, and their legacy environments probably aren't that legacy that SSE instructions aren't available. So it's reasonable to default to SSE.
- pwdisswordfish2 6y agoPython IIRC stores intermediate results in memory, it doesn't perform long chains of calculations entirely in FPU registers. This ensures the intermediate results are predictably rounded. JavaScript JITs blindly use SSE without checking if it is available. Bytecode interpreters behave the same as above.
- dpwm 6y agoMy understanding is that in this case it's GCC's optimizer evaluating the expression at compile time in a way that gives a different result to if it's done at runtime. EDIT: GCC gets it wrong with -O0 (when it's evaluated at runtime) and right with -O2 (when it evaluates it at compile time) Clang appears to get it right when evaluated at compile time or runtime. Clang uses the divsd instruction; GCC uses the fdiv instruction – so this really is SSE vs x87 FPU.
- nwallin 6y agoclang gets it "right" because it uses the sse divsd instruction to perform the division, gcc gets it "wrong" because it uses the x87 fdiv instruction. You can configure gcc to use sse math with the -mfpmath=sse instruction. The author states that gcc still gets it "wrong" but he simply is not correct. gcc will use divsd and gets 0.501782303180000055. gcc defaults to use the x87 fpu instead of sse because not all 32 bit x86 cpus have SSE instructions. It's a safer default. When compiled in 64 bit mode, it uses up to SSE2 instructions, because all x86_64 CPUs have SSE2. https://godbolt.org/z/3PfKQU https://godbolt.org/z/3PfKQU
- ndesaulniers 6y agoI think `-mno-see -mno-sse2` are of interest here, too.
- montalbano 6y agoLearning C recently, I ran into confusion around FLT_EVAL_METHOD as my dirty (not recommended) trick for getting machine epsilon didn't work (7.0/3.0 - 4.0/3.0 - 1.0). With strict IEEE 754 floating points it works fine (Python, Julia, Go all work as expected). But gcc, with c11 flags at least, performed the computation using the 80 bit upsampling mentioned in this post. Caused some head-scratching initially. Edit: Also, curiously, things worked as expected for me using -std=gnu11. This was a fairly old version of GCC though, things might have changed.
- mikedilger 6y agorust (default settings, amd64) gets 0.501782303180000055
- MayeulC 6y agoDoes Rust still uses LLVM? If so, I'd tend to expect a result similar to Clang.
- deleted 6y ago[deleted]
- torstenvl 6y agoClang with default settings seems to get it right: https://imgur.com/a/d3LBlLK https://imgur.com/a/d3LBlLK
- zucker42 6y agoDoes it work with GCC on your computer? This might be 32-bit specific.
- bonzini 6y agoIt may gets it right for this particular test case, but it will get others wrong. It's impossible to say "it gets it right" from a single testcase or three, you have to look at the compiler innards. See https://gcc.gnu.org/bugzilla/show_bug.cgi?id=323 https://gcc.gnu.org/bugzilla/show_bug.cgi?id=323 for the gory details; clang works the same as GCC's C++ frontend and the same that GCC <4.5 used to behave for C.
- 1anna69 6y agoHello guys my hot photos are here .. http://vivy.us/isxxx http://vivy.us/isxxx username annaxx))
- saagarjha 6y agoI believe the issue is that some of this ends up evaluated at compile-time, and the compiler really has no idea how to do rounding correctly at that point. When implementing floating point tests for an x86 emulator, we had to create a macro that would ensure that the calculation was done at runtime rather than at compile time for this exact reason: https://github.com/ish-app/ish/blob/18176b6931d69bdabe48f137e48d8e8118c22122/emu/float80-test.c#L11 https://github.com/ish-app/ish/blob/18176b6931d69bdabe48f137.... If anyone knows how to control floating point rounding in C reliably, I'd be interesting in hearing it :)
- deleted 6y ago[deleted]
- slavik81 6y agoIt's been a while since he's done that sort of thing, but Bruce Dawson could probably tell you some stories. I don't remember seeing that specific problem discussed, but his blog is a treasure trove of floating point information: https://randomascii.wordpress.com/category/floating-point/ https://randomascii.wordpress.com/category/floating-point/
- josefx 6y ago> the compiler really has no idea how to do rounding correctly at that point. Flags, pragmas, function attributes, there are so many ways this could be communicated to the compiler. If this behavior is not a bug, then having a way to change it should be a feature.
- avar 6y ago> If anyone knows how to control floating point rounding in C reliably, I'd be interesting in hearing it I skimmed over that code and it seems that hack is required because a) you're not happy with what the compiler does with a double b) insisting on using a "double" type. It's easy to get your own semantics for numbers in C, it's not easy to do so while insisting on using C's own types whose behavior are contrary to what you want. If you want your own floating point semantics why not make a struct with a significand & base and pass that around? Something like that is how established numbers-in-C-without-native-C-semantics libraries do it, e.g. GMP.
- amluto 6y agoPerhaps the title should mention that this is 32-bit x86, which is rather obsolete. This should work okay on x86_64.
- tbodt 6y agoAre you sure? The FPU instruction set is the same.
- dpwm 6y agoI just checked on the linked Compiler Explorer session [0], and if you take out the -m32 then GCC does indeed use divsd instead of fdiv. [0] https://godbolt.org/z/py3Dw0 https://godbolt.org/z/py3Dw0
- nwallin 6y agoYes. Toggle -m32 (compiles to 32 bit) or deleting it. (defaults to 64 bit, because it's not 2005 anymore) It also works with -m32 -mfpmath=sse, even though the author of the article incorrectly states otherwise. https://godbolt.org/z/py3Dw0 https://godbolt.org/z/py3Dw0
- olliej 6y agoThe platform ABI can differ on the same cpu. I originally had the wrong answer here, but I'm including below for historical purposes and embarrassment. Re-reading on my laptop meant I didn't misread/read-the-expected? Basically by default when compiling for IA32 gcc assumes that there is no SSE unit, and so must use the x87 fpu. Alas the x87 unit doesn't have either 32 or 64 bit ieee floating point types, so the calculation at runtime must be rounded to float/double at the end of the computation. That results in a double round that is incorrect. Minor additional comment - x87's reduced precision modes (that reduce mantissa precision to ieee754-32/64) don't reduce the exponent size so it's conceivable you could get multiple-rounding derived errors even then. But no one would be toggling x87 state like that anyway, so theoretically not something that would matter in most cases. Old incorrect answer: On ia32 long double is the x87 80bit ieee754 numeric type. On x86_64 on windows, and maybe Linux, long double is just a 64bit ieee float. So my guess is that GCC computes constant evaluation of floating point expressions in “long double” for maximum precision, and then rounds to the required precision at the end. This is not valid going from 64->32 bit floats and you could probably induce this bug if you tried hard enough. Anyway it is not sound to round twice in any numeric operation, and by doing the calculation at one precision and then rounding that to the required precision that is what gcc is doing. This is a bug in GCC, floating point is well defined, there isn’t randomness to it, and they’ve introduced additional rounding that is not valid.
- secondcoming 6y agoIf you add '-ffast-math' to the compiler parameters in his godbolt link you get the 'right' answer. Doing that means the division is done with the 'divsd' instruction, which clang uses even without the 'fast-math' flag. That doesn't mean it's necessarily the right thing to do though. Clang can be made generate the same code as the 'broken' gcc output by declaring the variables as 'long double' instead of 'double'. So I think this is where the issue lies? Edit: godbolt link requiring large monitor: https://godbolt.org/z/tFiZwW https://godbolt.org/z/tFiZwW
- saagarjha 6y ago-ffast-math should rebrand itself to -ffast-and-more-correct-math ;)
- garmaine 6y agoIt’s normally not more correct.
- saagarjha 6y ago;)
- all-fakes 6y agoI thought fast-math meant (among other things) that the compiler could algebraically simplify things using some rules for real numbers that don't apply to 32-/64-bit floating point numbers. I'd expect that reducing the number of operations would more often than not reduce the error due to rounding, which is what most people would want when they talk about "correctness". Which step in this (very naive) argument is wrong?
- saagarjha 6y agoThe code lacking -ffast-math doing math incorrectly, and -ffast-math ending up doing it correctly by accident.
- 6y ago
- sonofgod 6y agoAn example where this matters -- a new Mario 64 trick on the Virtual Console version because summing round-towards-zero numbers over a sine wave cycle produces a small delta which builds up over a couple of days to make a jump possible https://www.kotaku.com.au/2018/06/the-mario-64-trick-that-takes-three-days-to-complete/ https://www.kotaku.com.au/2018/06/the-mario-64-trick-that-ta...
- Konohamaru 6y ago> An example where this matters video games don't matter lol
- Spivak 6y agoAn example of real life software where this the difference has a noticeable impact.
- jmblpati 6y agoA small clarification: this trick is not new. It was discovered very early on in the investigation of the Wii VC release of Super Mario 64. It was rediscovered by a group of people who were interested in an application of it (which enables one to beat SM64 without the A button) a couple of years ago, and the discovery that the glitch is a result of inaccurate rounding mode emulation followed a couple weeks later.
- sonofgod 6y agoThank you for the clarification.
- dx87 6y agoFor anyone else interested in similar tricks, here is an in-depth video where someone talks about getting a star with half a button press. It starts out straightforward enough, but eventually gets really interesting when they start talking about parallell worlds that you can only access by abusing bugs to reach unrealistic speeds. https://www.youtube.com/watch?v=kpk2tdsPh0A&t=77s https://www.youtube.com/watch?v=kpk2tdsPh0A&t=77s
- Traster 6y agoIf you're interested in the accuracy of floating point operations, I recommend taking a look at the khronos openCL[1], it's not as clear cut as you would think. I think there are good reasons not to expect floating point ot be bit accurate. [1]:https://www.khronos.org/registry/OpenCL/specs/2.2/html/OpenCL_C.html#relative-error-as-ulps https://www.khronos.org/registry/OpenCL/specs/2.2/html/OpenC...
- alkonaut 6y agox87 is crazy. Especially since you can’t easily know when something is truncated from 80 bits to 64. Add a parameter to a function and an unrelated calculation behaves differently because it now ran out of 80bit registers. Congratulations, you now have a spurious test failure you will spend weeks trying to find. The only way of doing FP on x86 and keeping your sanity is by making sure you use 32/64 bit floats 100% of the time. That is: use SSE registers and never the 80-bit x87 arithmetic. I thought compilers for 64 bit programs these days did exactly this. I know the .NET Jit does for example.
- chrisseaton 6y agoI don't think any serious compiler for AMD64 is using x87 - it's ancient history. You would have to go deliberately far out of your way to encounter this problem.
- jabl 6y agoIf you use long double on x86-64 Linux, you're using x87. Though that has none of the problems discussed here, as there's no rounding when spilling regs to memory.
- chrisseaton 6y agoI think using long double counts as going 'deliberately far out of your way'.
- nwallin 6y agoThe article is about 32 bit gcc. 64 bit gcc defaults to SSE, although you can force x87 if you want to. (although... maybe don't do that)
- deleted 6y ago[deleted]
- brudgers 6y agoAs I understand it, the IEEE floating point standards do not require exact rounding for binary to decimal conversions. See What Every Computer Scientist Should Know About Floating-Point Arithmetic https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.h... The behavior appears to be consistent with the standard. Keep in mind that GCC needs to provide the same result irrespective of the architecture it runs on and one reason to use GCC for floating point is to take advantage of dedicated hardware. And the reason to use floating point is speed not precision (at least when using decimal floating point). If precision is critical then arbitrary precision or fixed point is the way to go. Normally it isn’t that critical.
- tijsvd 6y agoThe article is not about binary to decimal conversion. It's about rounding the result of a division, all in binary. A decimal equivalent would be like saying 200/3 = 66, so rounded towards zero instead of nearest. That is usually not problematic, but it can be in some situations, and in some of those, rounding to nearest can mitigate (for example, when you sum a large number of division results). The author points out that the behaviour should be configurable, but gcc doesn't seem to honor that.
- maxlybbert 6y agoLook for the heading, “The IEEE Standard” in “What Every Computer Scientist Should Know About Floating-Point Arithmetic,” then scroll down to the subheading “Operations.” That section begins, “The IEEE standard requires that the result of addition, subtraction, multiplication and division be exactly rounded. That is, the result must be computed exactly and then rounded to the nearest floating-point number (using round to even).” A few paragraphs later, it adds, “There is not complete agreement on what operations a floating-point standard should cover. In addition to the basic operations +, -, × and /, the IEEE standard also specifies that square root, remainder, and conversion between integer and floating-point be correctly rounded. It also requires that conversion between internal formats and decimal be correctly rounded (except for very large numbers).” I will admit that I worked as a programmer for several years without knowing that. I learned it from a former physics professor who has now moved on to become a quant with a hedge fund. He usually didn’t worry much about it, but when he needed to be sure a particular calculation was as accurate as possible, he’d start counting individual operations and stay away from particular function calls.
- upofadown 6y agoSo the complaint is that this loses half a bit of precision on average? Is that a lot?
- kevingadd 6y agoIf you lose half a bit on every operation and do multiple operations, that could be pretty bad. If it's half a bit at the end of a long chain of operations, that's tolerable and not worse than the precision specified by the docs for some execution environments (GPU shaders, for example) It also matters a lot whether the error is predictable and possible to compensate for.
- mjcohen 6y agoBack in the 90's I was working on some Ada code with some problems that were traced to conversion of strings to floating point at run time differing from conversion at compile time.
- known 6y agoDIY rounding function is given here; In unbiased rounding, .5 rounds to even, rather than always up, so 1.5 rounds to 2 but 4.5 rounds to 4 https://www.gnu.org/software/gawk/manual/html_node/Round-Function.html https://www.gnu.org/software/gawk/manual/html_node/Round-Fun...
- gcc_programmer 6y agoThe thread comments didn't dissapoint: > someone defended gcc > someone attacked it > someone who just started learning C or C++ made a comment > someone mentioned clang did it right > someone mentioned Rust
- gcc_programmer 6y agoand I got downvoted for not contributing meaningfully to the topic! The set is full
- helltone 6y agoI haven't seen anyone mentioning it, but I wonder if -ffloat-store fixes this?
- fancyfredbot 6y agoThe worst of this is that internal 80 bit x87 registers are truncated to 64 bit when written to memory. You have essentially no control over when this happens from the language level. It meant that floating point code was all rather undefined and you got different results on different compilers. As someone who has had to write multi platform numerical code I'm very thankful it's gone from AMD64.