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

[llvmbot]

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

[llvmbot]

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

[LebedevRI]

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

[llvmbot]

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.

[Muon]

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.

[Muon]

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.

[Muon]

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.

[hvdijk]

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.

[Muon]

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

[hvdijk]

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.

[Muon]

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

[hvdijk]

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.

[RalfJung]

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