Wrong optimization: instability of x87 floating-point results leads to nonsense

@topperc


Bugzilla Link	44873
Version	trunk
OS	Linux
Reporter	LLVM Bugzilla Contributor
CC	@topperc,@LebedevRI,@RKSimon,@rotateright,@vinc17fr

Extended Description

x87 floating-point results are effectively unstable due to possible excess precision. Without extra care, this instability can taint everything around and lead to nonsense.

Instability is not limited to FP numbers, it extends to integers too:

#include <stdio.h>

__attribute__((noipa,optnone)) // imagine it in a separate TU
static int opaque(int i) { return i; }

int main()
{
    int z = opaque(1) + 0x1p-60 == 1;

    printf("z = %d\n", z);
    if (z) 
        puts("z is one");
}

$ clang -std=c11 -pedantic -Wall -Wextra -Wno-unknown-attributes -m32 -march=i686 -O3 test.c && ./a.out
z = 0
z is one

----------------------------------------------------------------------
clang x86-64 version: clang version 11.0.0 (https://github.com/llvm/llvm-project 14ecbd7b8ded18af6c95f6a9957da541d1ec0e80)
----------------------------------------------------------------------

And instability of integers then easily taints surrounding code:

#include <stdio.h>

__attribute__((optnone)) // imagine it in a separate TU
static int opaque(int i) { return i; }

static void f(int a)
{
    int z = opaque(0);

    printf("z = %d\n", z);
    if (z == a && a == 1)
        printf("z = %d\n", z);
}

int main()
{
    f(opaque(1) + 0x1p-60 == 1);
}

$ clang -std=c11 -pedantic -Wall -Wextra -m32 -march=i686 -O3 test.c && ./a.out
z = 0
z = 1

----------------------------------------------------------------------
clang x86-64 version: clang version 11.0.0 (https://github.com/llvm/llvm-project 14ecbd7b8ded18af6c95f6a9957da541d1ec0e80)
----------------------------------------------------------------------

gcc bug -- https://gcc.gnu.org/bugzilla/show_bug.cgi?id=93681

int z = opaque(1) + 0x1p-60 == 1;

0x1p-60 there is just a double
https://godbolt.org/z/mJS3RN
It's as expected if it is long double: https://godbolt.org/z/NVvJht

int z = opaque(1) + 0x1p-60 == 1;
0x1p-60 there is just a double
https://godbolt.org/z/mJS3RN
It's as expected if it is long double: https://godbolt.org/z/NVvJht

Yeah, 1 + 0x1p-60 is the simplest expression I could think of where both operands fit into double exactly but the result doesn't.
x87 evaluates it as long double and gets the exact result. But if the result is spilled into memory (as double) it's rounded to just 1.

If you force the expression to honest long double (e.g., by casting any of the operands to it) then everything becomes stable -- the result both in a register and in memory has the same type and the same value. The same happens with SSE -- the result is stored as double everywhere and there are no discrepancies.

This issue in particular should be resolvable by either always storing the result to memory before performing a comparison, or by setting the x87 FPU's precision to the destination's precision.

@d0k

I'm trying to fix the contradictory control flow here by promoting f32 and f64 ops to intel_fp80 on the x87. This wouldn't make the numerical results conform to IEEE 754, but it would at least make them consistent. However, I'm running into trouble because it appears that LLVM doesn't actually properly support promotion for floating-point operations. I tried removing the f32/f64 registers from the x87, setting the operations to Promote, and then fiddling with this part of TargetLoweringBase.cpp:

llvm-project/llvm/lib/CodeGen/TargetLoweringBase.cpp

Lines 1414 to 1419 in 14bc11a

    
           if (!isTypeLegal(MVT::f32)) { 
        
             NumRegistersForVT[MVT::f32] = NumRegistersForVT[MVT::i32]; 
        
             RegisterTypeForVT[MVT::f32] = RegisterTypeForVT[MVT::i32]; 
        
             TransformToType[MVT::f32] = MVT::i32; 
        
             ValueTypeActions.setTypeAction(MVT::f32, TypeSoftenFloat); 
        
           }

But now I'm getting errors because it's trying to promote constants to intel_fp80, and I don't know how to tell it that the x87 can load those from memory as a widening conversion.

Can anyone point me in the right direction? Pinging @d0k @topperc since they were the last to touch that part of the code.

It is sufficient to read from a global variable to reproduce this (godbolt):

#include <stdio.h>

float c = 1.0f;

static void f(int z) {
    printf("z = %d\n", z);
    if (z) 
        printf("z = %d\n", z);
}

int main(void) {
    f(c + 0x1p-52f == c);
}

This version should also reproduce correctly on Windows as well.

I'm trying to fix the contradictory control flow here by promoting f32 and f64 ops to intel_fp80 on the x87.

I'm looking at this and realising that that can avoid the instability, but cannot generate correct results. For f32, it can work but for f64, it cannot. Because intel_fp80 has less than double the precision of f64, it is possible to get double rounding errors where adding two f64 values gives a different result than promoting both to intel_fp80, performing the addition there, and then rounding back to f64. It would still be an improvement, but the work needed to get that done might be work that would need to be thrown out again later.

Yes, I'm only trying to fix the control flow issue for now. Actually getting correct results will involve a performance hit, and my previous inquiries suggest that this would only be merged if it was behind a compiler flag.

Also, are you thinking of products? Twice the precision is enough to represent products exactly. However, I'm almost certain you need much more than twice the precision to make double rounding give correct results for addition, if it's possible at all. (I don't have a specific counterexample on hand, but I can probably come up with something if I fiddle with the exponents and rounding breakpoints...)

Yes, I'm only trying to fix the control flow issue for now. Actually getting correct results will involve a performance hit, and my previous inquiries suggest that this would only be merged if it was behind a compiler flag.

I do not see how it's possible to fix even that without a performance hit, in general, but I can imagine that there is a possibility that this performance hit will be lower than that of any alternative.

Also, are you thinking of products? Twice the precision is enough to represent products exactly. However, I'm almost certain you need much more than twice the precision to make double rounding give correct results for addition, if it's possible at all. (I don't have a specific counterexample on hand, but I can probably come up with something if I fiddle with the exponents and rounding breakpoints...)

For addition, first consider round to nearest, because it is the most complicated rounding mode. For round to nearest, in order for double rounding to happen, and to produce wrong results, you need the round to extended precision to be inexact, and produce a value which is exactly at a halfway point for the reduced precision which then gets rounded the wrong way in the second round operation. But producing a value which is exactly at a halfway point requires the most significant extended precision bit to be set after the round to extended precision, which requires either that or the next bit to be set in the infinitely precise original result. Meanwhile, having the round to extended precision be inexact requires a bit after the extended precision bits to be set in the infinitely precise result. These cannot both be true when performing an addition if the extended precision adds sufficient bits. I'm not sure how many extra bits of precision are needed exactly, but it should be little more than twice the original precision. For the other rounding modes, it's much simpler, it's either round to zero or round away from zero (round up and round down are one of those, depending on the sign), round to zero is just discarding the extra precision bits, and round away from zero, if inexact, ensures that at least one bit will be set in the extra precision bits so that the second conversion rounds correctly. I don't believe double rounding should ever be an issue with fp32 additions promoted to fp80, in any rounding mode supported by x87. Happy to be corrected if wrong, of course.

That said, this is kind of distracting from the point, which is not how much precision is enough for addition to be safe in extended precision, it's that fp80 is definitely not precise enough for fp64 operations to be safely promoted to fp80.

The only way I think we can get correct results is if we let LLVM alter the x87 control word to set the precision to that needed by operations. And I actually suspect that if done right, that will have a lower impact on performance in code which does not mix different FP types, because it means LLVM can just set the precision at the start of a function, perform all the operations it needs, and then restore the precision at the end, whereas leaving x87 in its default precision would require extra instructions to ensure correct rounding after every operation. Code which does mix different FP types would need the control word to be changed throughout the function though.

I agree that a true fix would be immensely preferable. However, I want to add that it's not sufficient to set the control word to the right precision, because that doesn't modify the exponent range of the x87 registers. You still need to do extra ops to ensure correct handling of underflow, subnormals, and overflow.

Also, yes, there is a performance cost to promoting all ops to fp80. Since changing the control word plays havoc with pipelining, I suspect this is cheaper, but only profiling will tell.

(Also, indeed, double rounding is never a problem for directed rounding modes.)

However, I want to add that it's not sufficient to set the control word to the right precision, because that doesn't modify the exponent range of the x87 registers.

Thanks for that addition. You are correct, we need that as well.

I'm looking at this and realising that that can avoid the instability, but cannot generate correct results. For f32, it can work but for f64, it cannot. Because intel_fp80 has less than double the precision of f64, it is possible to get double rounding errors where adding two f64 values gives a different result than promoting both to intel_fp80, performing the addition there, and then rounding back to f64. It would still be an improvement, but the work needed to get that done might be work that would need to be thrown out again later.

Yeah a naive translation will give wrong results. But getting the exact results seems to be possible; Java is supposedly precisely emulating fp64 on the x87. See here and here for details.

Also, are you thinking of products? Twice the precision is enough to represent products exactly. However, I'm almost certain you need much more than twice the precision to make double rounding give correct results for addition, if it's possible at all. (I don't have a specific counterexample on hand, but I can probably come up with something if I fiddle with the exponents and rounding breakpoints...)

There's a general theorem that twice the precision is good enough, AFAIK. Here's the paper for that: https://hal.science/hal-01091186/document

This issue in particular should be resolvable by either always storing the result to memory before performing a comparison, or by setting the x87 FPU's precision to the destination's precision.

I don't think it will work in any case. With float32, there are no problems setting the destination precision as a 64 bit mantissa is more than twice the size. With float64 however, it seems to be very tricky to resolve. Mainly the double rounding issue seems to occur whether you set the precision mode or not. So float64 at 80-bit precision has the double rounding issue already mentioned.

But float64 at 64-bit precision control has double rounding for subnormal results. At issue is the x87 does not seem to properly handle subnormals - it only controls the 24/53 bit precision of the mantissa, but does not control the exponent. So subnormals in float64 will be normals in an 80-bit float and still suffer double rounding. For some reason I do not find the subnormal issue mentioned anywhere (though it is not hard to find some cases).

Here is an implementation of float64 multiplication on x87 that produces a correctly rounded result for all x87 rounding modes in LLVM IR:
https://godbolt.org/z/84KPes9ac

[IR] LangRef: state explicitly that floats generally behave according…

lang/rust: enable SSE2 by default on i386

Would it be workable to add a flag like -mfpmath=387-precise to go with -mfpmath=sse/-mfpmath=387 that save and restore the CW around each op and uses the Java implementations? This loses some efficiency (add/sub goes from 3 ops to ~6 and mul/div 3->~8) but at least gets a path for the right semantics and sound behavior off the ground that's more tolerable than soft float. In theory, an optimization could merge the save/restore around multiple ops if they aren't interspersed with something opaque.

I think it would be possible to do that combined with something similar to the TypeSoftPromoteHalf currently used for f16 (saving and restoring the CW shouldn't be needed around basic operations for f32, just f64). I don't think it should be an optional flag: the current behaviour is so broken (as it explicitly conflicts with the semantics of LLVM IR) that I don't think it's worth keeping.

I have an (admittedly currently mostly LLM vibe-fixed, as a PoC) branch that makes the LLVM and clang behavior align with GCC behavior with -fexcess-precision=standard -ffloat-store set. That would also at least align it with the stated behavior in the documentation:

Clang does not use excess precision arithmetic for most types or on most targets. For example, even on pre-SSE X86 targets where float and double computations must be performed in the 80-bit X87 format, Clang rounds all intermediate results correctly for their type.

Clang's behavior right now is pretty much equivalent to GCC with -fexcess-precision=fast (and without -ffloat-store). With those settings the same misbehavior (z = 0\nz is one) can even be observed in GCC.

Now I understand that this doesn't actually fix IEEE compliance fully (like the Java way would), but would there be appetite to bring that behavior in?
My main concern would be that it might make it actually harder to add the proper Java-like implementation later.

(saving and restoring the CW shouldn't be needed around basic operations for f32, just f64)

(wouldn't there still be double-rounding issues if the CW is set to 24bit for some reason?)

Would it be workable to add a flag like -mfpmath=387-precise to go with -mfpmath=sse/-mfpmath=387 that save and restore the CW around each op and uses the Java implementations? This loses some efficiency (add/sub goes from 3 ops to ~6 and mul/div 3->~8) but at least gets a path for the right semantics and sound behavior off the ground that's more tolerable than soft float. In theory, an optimization could merge the save/restore around multiple ops if they aren't interspersed with something opaque.

You could probably mitigate most of the performance issues just by making clang implement FLT_EVAL_METHOD==2 correctly (emit all intermediate operations as x86_fp80). But I'm not particularly sure that anyone is tracking 32-bit x87 performance numbers.

(wouldn't there still be double-rounding issues if the CW is set to 24bit for some reason?)

If the function is non-strictfp, it is UB to be in non-default FP env, and PC set to single-precision is definitely non-default FP env (Windows makes the default PC be double precision while SysV makes the default PC be extended precision).

I think it would be possible to do that combined with something similar to the TypeSoftPromoteHalf currently used for f16 (saving and restoring the CW shouldn't be needed around basic operations for f32, just f64). I don't think it should be an optional flag: the current behaviour is so broken (as it explicitly conflicts with the semantics of LLVM IR) that I don't think it's worth keeping.

What about TypeSoftPromoteHalf would be relevant - you mean just how the backend expands the ops? Good point about f32. Edit: though with strictfp I guess this needs to be saved+restored anyway

Edit: about TypeSoftPromoteHalf I guess this may just mean how each individual operation is wrapped in an extend/trunc similar to how the fadd/fsub/fmul/fdiv needs to be wrapped in some other instructions here

I'd of course prefer it to be non-optional as well but it seemed like there may have been some reluctance earlier in the thread. But...

You could probably mitigate most of the performance issues just by making clang implement FLT_EVAL_METHOD==2 correctly (emit all intermediate operations as x86_fp80). But I'm not particularly sure that anyone is tracking 32-bit x87 performance numbers.

Meaning users get the correct semantics (java ops) by default but could effectively opt into the old semantics by setting FLT_EVAL_METHOD? That sounds great to me.

Now I understand that this doesn't actually fix IEEE compliance fully (like the Java way would), but would there be appetite to bring that behavior in?

I am not a LLVM maintainer but my 2cents: if the fix is small and pretty easy, or if something from it is going to be needed with the switch to java semantics, then it's helpful enough. But if it's not going to be useful once the correct semantics are implemented then I think it would be better to just work on that instead because that's going to completely resolve more problems.

Meaning users get the correct semantics (java ops) by default but could effectively opt into the old semantics by setting FLT_EVAL_METHOD? That sounds great to me.

Users would get a correct semantics, but not necessarily the one they intend. FLT_EVAL_METHOD is set by the compiler, not the user, although there is a command-line flag that could be set to change its value. More likely, i386 compilation that requests the x87 FPU would want the default flag to be to use extended precision for everything, and users would have to opt into to full IEEE 754 conformance (much like how Java made you opt into strictfp) by resetting the flag to no extra precision.

The solution I would (still) like to see is a command-line flag which causes all of the float/double operations to call library functions instead of lowering directly to x87 instructions.

I don't really feel like it's really justified, at this point, to add complex lowerings for these operations in the compiler, especially given the code-size of the required sequences. Nor do I think it's really justified to do anything complex like mode change coalescing or anything like that.

I definitely think we ought to offer the mode which gets this right, but realistically, the userbase for 32-bit x86 without SSE2 is practically nil, so we should do it in the easiest way possible.

@beetrees

What is the concern with backend complexity? To me this doesn't seem much worse than f16 ops, which typically need a conversion to/from f32, or something like minimum that's ~14 instructions on x86. And if the correct behavior were to be opt-in as you suggested, changes in binary size wouldn't be much of a concern?

I don't really think optimizing the patterns are required but if it's something that may come nearly for free at some point (think this has been talked about in the context of strictfp?) then any drawbacks of switching the default pretty much go away.

Regardless, I do agree with @beetrees that it's worth changing the default anyway - the miscompilations are nasty. I really can't imagine that anybody would complain about a few extra instructions if they get to stop special casing i586 in math libraries and unskip a number of broken tests. That's just not a sane tradeoff if you're writing cross-platform float code and the alternative is handwritten x87 assembly. (Of course, no decision is needed here before such a flag even exists.)

Alright, I posted #215104, which should fix the instability issue, and provide a FLT_EVAL_METHOD=0 and FLT_EVAL_METHOD=2 with the ~same caveats as GCC (i.e. memory roundtrips for rounding; doesn't fix the issues that need the Java-style solution).

[IR] LangRef: state explicitly that floats generally behave according…

Navigation Menu

Repository navigation

Wrong optimization: instability of x87 floating-point results leads to nonsense #44218

Extended Description

12 remaining items

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

	if (!isTypeLegal(MVT::f32)) {
	NumRegistersForVT[MVT::f32] = NumRegistersForVT[MVT::i32];
	RegisterTypeForVT[MVT::f32] = RegisterTypeForVT[MVT::i32];
	TransformToType[MVT::f32] = MVT::i32;
	ValueTypeActions.setTypeAction(MVT::f32, TypeSoftenFloat);
	}

Wrong optimization: instability of x87 floating-point results leads to nonsense #44218

Description

Extended Description

Activity

llvmbot commented on Feb 11, 2020

LebedevRI commented on Feb 11, 2020

llvmbot commented on Feb 11, 2020

Muon commented on Feb 16, 2023

Muon commented on Oct 21, 2023

Muon commented on Nov 23, 2023

hvdijk commented on Nov 28, 2023

Muon commented on Nov 28, 2023

hvdijk commented on Nov 29, 2023

Muon commented on Nov 29, 2023

hvdijk commented on Nov 29, 2023

RalfJung commented on Nov 29, 2023

12 remaining items

GregoryMorse commented on Aug 2, 2025

johnplatts commented on Nov 12, 2025

tgross35 commented on Aug 6, 2026

beetrees commented on Aug 6, 2026

seritools commented on Aug 6, 2026

jcranmer-intel commented on Aug 6, 2026

tgross35 commented on Aug 6, 2026

jcranmer-intel commented on Aug 7, 2026

jyknight commented on Aug 7, 2026

tgross35 commented on Aug 9, 2026

seritools commented on Aug 9, 2026

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions