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