How is fma() implemented
floating-point, fma, ieee-754, instruction-set
Solution
The actual implementation varies from platform to platform, but speaking very broadly:
If you tell your compiler to target a machine with hardware FMA instructions (PowerPC, ARM with VFPv4 or AArch64, Intel Haswell or AMD Bulldozer and onwards), the compiler may replace calls to `fma( )` by just dropping the appropriate instruction into your code. This is not guaranteed, but is generally good practice. Otherwise you will get a call to the math library, and:
When running on a processor that has hardware FMA, those instructions should be used to implement the function. However, if you have an older version of your operating system, or an older version of the math library, it may not take advantage of those instructions.
If you are running on a processor that does not have hardware FMA, or you are using an older (or just not very good) math library, then a software implementation of FMA will be used instead. This might be implemented using clever extended-precision floating-point tricks, or with integer arithmetic.
The result of the `fma( )` function should always be correctly rounded (i.e. a "real fma"). If it is not, that's a bug in your system's math library. Unfortunately, `fma( )` is one of the more difficult math library functions to implement correctly, so many implementations have bugs. Please report them to your library vendor so they get fixed!
Is there an intrinsic to make sure that a real FMA is used, when the precision is relied upon?
Given a good compiler, this shouldn't be necessary; it should suffice to use the `fma( )` function and tell the compiler what architecture you are targeting. However, compilers are not perfect, so you may need to use the `_mm_fmadd_sd( )` and related intrinsics on x86 (but report the bug to your compiler vendor!)
Problem
According to the documentation, there is a `fma()` function in `math.h`. That is very nice, and I know how FMA works and what to use it for. However, I am not so certain how this is implemented in practice? I'm mostly interested in the `x86` and `x86_64` architectures. Is there a floating-point (non-vector) instruction for FMA, perhaps as defined by IEEE-754 2008? Is FMA3 or FMA4 instruction used? Is there an intrinsic to make sure that a real FMA is used, when the precision is relied upon?