GNU GCC does not round floating-point divisions to the nearest value

https://lemire.me/blog/2020/06/26/gcc-not-nearest/

by guiambros • 6 years ago

158 78 6 years ago

wiml• 6 years ago

This 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• 6 years ago

But then why does python and JavaScript on the same machines get it right?

nitrogen• 6 years ago

Unsure: compiled to use SSE instead of x87 for fpmath?

stefan_• 6 years ago

Because 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.

nwallin4• 6 years ago

Because 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• 6 years ago

Python 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.

dpwm2• 6 years ago

My 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.

nwallin4• 6 years ago

clang 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

montalbano• 6 years ago

Learning 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• 6 years ago

rust (default settings, amd64) gets 0.501782303180000055

MayeulC• 6 years ago

Does Rust still uses LLVM? If so, I'd tend to expect a result similar to Clang.

torstenvl• 6 years ago

Clang with default settings seems to get it right:

https://imgur.com/a/d3LBlLK

zucker42• 6 years ago

Does it work with GCC on your computer? This might be 32-bit specific.

bonzini4• 6 years ago

It 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 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.

saagarjha5• 6 years ago

I 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/18176b6931d69bdabe48f137.... If anyone knows how to control floating point rounding in C reliably, I'd be interesting in hearing it :)

slavik81• 6 years ago

It'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/

josefx• 6 years 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• 6 years 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.

amluto3• 6 years ago

Perhaps the title should mention that this is 32-bit x86, which is rather obsolete. This should work okay on x86_64.

tbodt• 6 years ago

Are you sure? The FPU instruction set is the same.

dpwm2• 6 years ago

I 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

nwallin4• 6 years ago

Yes. 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

olliej5• 6 years ago

The 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.

chrisseaton3• 6 years ago

Who's going to use FPU on AMD64? Why on earth would you do that?

garmaine2• 6 years ago

Why would there be any difference?

olliej5• 6 years ago

My guess is that they’re using long double, and the abi definition of long double may be different on x86_64 (eg on windows long double is logically just double)

amluto3• 6 years ago

Because every x86_64 CPU supports SSE2, so compilers can assume it exists and the ABI passes floating point arguments in SSE2’s FP registers. This means that every sane compiler for x86_64 will use SSE2 and won’t use x87, and the problem won’t occur.

In contrast, the 32-bit C ABI passes floating point args in x87 registers, and there are 32-bit CPUs without SSE2, so x87 is used by default.

x87 is, in quite a few respects, awful.

secondcoming• 6 years ago

If 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

saagarjha5• 6 years ago

-ffast-math should rebrand itself to -ffast-and-more-correct-math ;)

garmaine2• 6 years ago

It’s normally not more correct.

sonofgod3• 6 years ago

An 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-ta...

jmblpati• 6 years ago

A 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.

sonofgod3• 6 years ago

Thank you for the clarification.

dx87• 6 years ago

For 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

sonofgod3• 6 years ago

Strongly recommended watching, in my opinion: very clearly describes some quite complicated manipulations, the circumstances in which they work and how to move between those circumstances, to achieve what should be impossible.

If you're going to watch one commentated "tool-assisted" run, this is it.

LambdaComplex• 6 years ago

Best followed by this video, which made me laugh my ass off for several minutes: https://www.youtube.com/watch?v=hSa-pwML4z4

(I grew up watching SpongeBob, if that matters)

caseymarquis• 6 years ago

Fun idea: If the universe is a simulation, then maybe dark energy is just a rounding error bug.

hyperman1• 6 years ago

According to Wikipedia:

   ..dark energy contributes 68% of the total energy..

so if anyone complains about the shoddy quality of x87 floats, point them to the universe and tell them it could be worse.

Traster• 6 years ago

If 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/OpenC...

alkonaut• 6 years ago

x87 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.

chrisseaton3• 6 years ago

I 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• 6 years ago

If 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.

nwallin4• 6 years ago

The 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)

brudgers2• 6 years ago

As 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.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• 6 years ago

The 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• 6 years ago

Look 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.

brudgers2• 6 years ago

Yes, you are correct. Thanks. This stuff is hard for me.

Reading down thread, I see that the behavior might be related to x87 since the build script references i386...which makes IEEE 754 conformance moot. The methods of achieving conformance when targeting i386 are described in the GCC documentation.

upofadown• 6 years ago

So the complaint is that this loses half a bit of precision on average? Is that a lot?

kevingadd• 6 years ago

If 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• 6 years ago

Back 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• 6 years ago

DIY 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-Fun...

gcc_programmer2• 6 years ago

The 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_programmer2• 6 years ago

and I got downvoted for not contributing meaningfully to the topic! The set is full

helltone• 6 years ago

I haven't seen anyone mentioning it, but I wonder if -ffloat-store fixes this?

fancyfredbot• 6 years ago

The 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.

GNU GCC does not round floating-point divisions to the nearest value

Related Stories