4 ms·
Okay I got nerdsniped by this, reimplemented and reproduced it. As others mentioned, the article does not show a flaw in the M1's FMA. The results are the same
by thxg 5y ago
Okay I got nerdsniped by this, reimplemented and reproduced it. As others mentioned, the article does not show a flaw in the M1's FMA. The results are the same on x86_64. It is just due to how this particular formula interacts with floating-point rounding, and it all depends on the order of magnitude of the numbers.
The basic operation here is
x = (((a * b) + c) - c) / b
which would yield (x == a) in exact arithmetic.
Let's denote by fp(y) the floating-point number closest to y. If you don't use FMA, then the middle of the formula is essentially fp(fp(M + c) - c) for some M. When you pick the natural choice of initial a, b and c having the same order of magnitude, this gives exactly M in 99.99% of cases. Then you're just measuring the difference between a and (a * b) / b. Instead, FMA gives you x = fp(fp(fp(a * b + c) - c) / b), which is admittedly more complex to study.
But as you pick larger and larger c values, you can observe that FMA becomes more accurate compared to the other version. When c reaches about 0.5e+6 times larger than a and b, FMA starts winning. Note that a is incremented 1M times as the proposed algorithm iterates, so this corresponds to (a * b) and c having the same orders of magnitude on average, as the loop progresses.
You can probably go much more in-depth in studying this, and explain it better. But essentially this is a non-issue, just floating-point being floating-point.
If you want to play with the code, you can find it here:
https://godbolt.org/z/WWqdcrzdh https://godbolt.org/z/WWqdcrzdh
- thxg 5y agoAlso, the article concludes that FMA hurts accuracy, and this is at best misleading (FMA is better in the worst case). But the example given does highlight something important: Depending on the case, the accuracy gains can be tiny (so much so that in this example the FMA version does indeed return larger errors in most runs), and the performance gains can be small too. On the other hand, the price to pay for letting the compiler freely decide whether/when to fuse * and + is huge: reproducibility. Since x87 was deprecated, having floating-point results be bit-for-bit consistent across platforms, compilers and optimization flags has been amazing for debugging. Losing that is a step back and a recipe for headaches. I was surprised that gcc enables -fexcess-precision=fast by default (in GNU mode), hence enabling transparent FMA and potential headaches. Specifying the standard (e.g. -stdc=c17) fixes that.
- dnautics 5y agoIirc x87 was even worse on some processors depending on the processor state, the processor could cache the CPU state on context switch and revert the x87 register to a f64, leading to irreproducibility at the runtime level (randomly fails if you get an unlucky thread thrash)
- rualca 5y ago> When c reaches about 0.5e+6 times larger than a and b, FMA starts winning. It should be noted that 1e+6 sticks out as right on the rounding limit of single precision floating point numbers, which is about 1e+7. Thus it suggests that the residue amplification from b ceases to be a problem once b simply falls out of the picture as it gets close to a roundoff error. Therefore, this test isn't evaluating FMA as much as it is showcasing how floating point division amplifies numerical residues.
- oceanghost 5y agoTo complement what you're saying... It took me a long time to understand: The number line for floats is non-linear. It has holes in it.