`fma` not defined for `BFloat16`

[dzhang314]

Does anyone know the correct way to define fma for BFloat16? Based on my understanding of the double-rounding theorems, I think it should be correct to simply cast to Float32 and back:

@inline Base.fma(x::BFloat16, y::BFloat16, z::BFloat16) =
    BFloat16(fma(Float32(x), Float32(y), Float32(z)))

But I've observed that this sometimes returns different results from BFloat16(Float32(x) * Float32(y) + Float32(z)), which puzzles me because I think they ought to be equivalent in round-to-nearest-even (provided that Float32 has more than twice plus two the precision of BFloat16, which is true).

[dzhang314]

Scratch that -- I see my mistake, and BFloat16(Float32(x) * Float32(y) + Float32(z)) is definitely not a correct implementation. The following implementations should all be correct:

• BFloat16(fma(Float32(x), Float32(y), Float32(z)))
• BFloat16(fma(Float64(x), Float64(y), Float64(z)))
• BFloat16(Float64(x) * Float64(y) + Float64(z))

You're applying the double-rounding bound twice, so if you round-to-nearest-even after every operation, you need an intermediate format with roughly 4x the final desired precision.

[giordano]

Recently there was a related issue in llvm (apparently llvm gets it wrong, also for float16): llvm/llvm-project#131531