APFloat's addOrSubtractSignificand still throws assert in subtler FMA cases.

[eddyb]

#include <math.h>

float foo() {
    return fmaf(
        -0.000000000000000000000000000000000000014728589,
        0.0000037105144,
        0.000000000000000000000000000000000000000000055
    );
}

results in (e.g. via godbolt):

llvm/lib/Support/APFloat.cpp:1805:
  llvm::lostFraction llvm::detail::IEEEFloat::addOrSubtractSignificand(const llvm::detail::IEEEFloat&, bool):
    Assertion `!carry' failed.

However, that's not how the bug was found. I noticed 8-bit formats (like Float8E5M2 and Float8E4M3FN) were added to APFloat, and decided to try brute-forcing all possible inputs for a few common operations (it only takes 8 seconds, and ~98% of that is FMA, because it uniquely has 3 inputs).

So the first example I found was for Float8E4M3FN, namely: FMA(0.0254, 0.0781, -0.00195) (the encoded byte values being 0x0d, 0x1a, 0x81).

Then @solson helped me turn some formula sketches I made looking at the code, into Z3, which confirmed that the example I found was the only one across all possible 8-bit floats, but 16-bit and 32-bit have a lot more cases. (I'm only not including all that here because it's not really set up to generate ready to use examples, it's really finicky as-is)

---

My understanding is that e62555c fixed most of the previously problematic cases, but not the ones where both significands are equal before subtraction, which will only work if lost_fraction == lfExactlyZero.

(cc @ekatz)

But if there is a lost fraction, neither direction will work to subtract equal significands - the code seems to rely on being able to assume that only equal exponents can result in equal significands, but we know that's not true w/ FMA.

(also, the lost fraction is always subtracted, regardless of the subtraction direction, but before e62555c the source of the lost_fraction was tied to reverse, i.e. a.subtractSignificand(b, lost_fraction != lfExactlyZero) was always performed with lost_fraction coming from lost_fraction = b.shiftSignificandRight(...);, regardless of how a and b mapped to *this and temp_rhs, respectively - this seems significant as well, but I'm not sure how to tell or test this)

[eddyb]

Just noticed LLVM does have exhaustive testing for 8-bit floats, neat!

llvm-project/llvm/unittests/ADT/APFloatTest.cpp

Lines 5073 to 5080 in 83f3920

	TEST(APFloatTest, Float8ExhaustivePair) {
	// Test each pair of 8-bit floats with non-standard semantics
	for (APFloat::Semantics Sem :
	{APFloat::S_Float8E4M3FN, APFloat::S_Float8E5M2FNUZ,
	APFloat::S_Float8E4M3FNUZ, APFloat::S_Float8E4M3B11FNUZ}) {
	const llvm::fltSemantics &S = APFloat::EnumToSemantics(Sem);
	for (int i = 0; i < 256; i++) {
	for (int j = 0; j < 256; j++) {

However, AFAICT that doesn't cover FMA - which is entirely understandable, as it would increase the test duration ~50x compared to all binary ops combined.

Sadly, FMA is also more likely to hit weird edge cases (some of which have asserts), thanks to its intermediary multiplicative result's expanded dynamic range. So maybe there should at least be some disabled-by-default long-running test that could be manually run once in a while? (I would expect it to take a few seconds per format)

cc @d0k @reedwm @krzysz00 @majnemer (based on git blame on the existing exhaustive tests)

[eddyb]

Something I should've checked is the non-assertion output (try on godbolt):

Compiler	printf("%.50f\n", foo()); output
Clang 16.0	-0.00000000000000000000000000000000077037197775489434
GCC 13.1	-0.00000000000000000000000000000000000000000000000000

Which makes sense, the assert is guarding overflow (i.e. it's wrapping around within the significand APInt).

[eddyb]

I've let a bruteforce run for f16 (aka IEEEhalf) run overnight and got 12432 samples out of it.
(All of them confirmed to trigger the assert(!carry), in the current state of APFloat, i.e. b0abd48)

I've also went ahead added parallelism to my bruteforce tool:

• fuzz: allow bruteforce-ing wider formats and add parallelism. rust-lang/rustc_apfloat#5

and as such, I can now generate such examples much faster, if that's of any use.

I'm attempting a run for f32 (aka IEEEsingle), and thanks to Z3 I know to expect some results for FMA(_, 5, 1) and FMA(_, 9, 1) (bit patterns, not float values - i.e. some of the first few subnormals), which should hopefully be reached within mere hours.
_{(as latter inputs come from higher order iteration bits, so in FMA(a, b, c), it effectively acts like for(c : ...) for(b : ...) for(a : ...) test(a, b, c);)}