std::fma, std::fmaf, std::fmal

Parameters

Return value

If successful, returns the value of x * y + z

If a range error due to overflow occurs, ±HUGE_VAL, ±HUGE_VALF, or ±HUGE_VALL is returned.

If a range error due to underflow occurs, the correct value (after rounding) is returned.

Error handling

Errors are reported as specified in math_errhandling.

If the implementation supports IEEE floating-point arithmetic (IEC 60559),

• If x is zero and y is infinite or if x is infinite and y is zero, and
  • if z is not a NaN, then NaN is returned and FE_INVALID
  • if z is a NaN, then NaN is returned and FE_INVALID
• If x * y is an exact infinity and z is an infinity with the opposite sign, NaN is returned and FE_INVALID
• If x or y are NaN, NaN is returned.
• If z is NaN, and x * y is not 0 * Inf or Inf * 0, then NaN is returned (without FE_INVALID

Notes

This operation is commonly implemented in hardware as fused multiply-add CPU instruction. If supported by hardware, the appropriate FP_FAST_FMA?

POSIX (fma, fmaf, fmal) additionally specifies that the situations specified to return FE_INVALID

Due to its infinite intermediate precision, std::fma is a common building block of other correctly-rounded mathematical operations, such as std::sqrt or even the division (where not provided by the CPU, e.g. Itanium

As with all floating-point expressions, the expression x * y + z may be compiled as a fused multiply-add unless the #pragma STDC FP_CONTRACT

The additional overloads are not required to be provided exactly as (A). They only need to be sufficient to ensure that for their first argument num1, second argument num2 and third argument num3

Example

#include <cfenv>
#include <cmath>
#include <iomanip>
#include <iostream>
 
#ifndef __GNUC__
#pragma STDC FENV_ACCESS ON
#endif
 
int main()
{
    // demo the difference between fma and built-in operators
    const double in = 0.1;
    std::cout << "0.1 double is " << std::setprecision(23) << in
              << " (" << std::hexfloat << in << std::defaultfloat << ")\n"
              << "0.1*10 is 1.0000000000000000555112 (0x8.0000000000002p-3), "
              << "or 1.0 if rounded to double\n";
 
    const double expr_result = 0.1 * 10 - 1;
    const double fma_result = std::fma(0.1, 10, -1);
    std::cout << "0.1 * 10 - 1 = " << expr_result
              << " : 1 subtracted after intermediate rounding\n"
              << "fma(0.1, 10, -1) = " << std::setprecision(6) << fma_result << " ("
              << std::hexfloat << fma_result << std::defaultfloat << ")\n\n";
 
    // fma is used in double-double arithmetic
    const double high = 0.1 * 10;
    const double low = std::fma(0.1, 10, -high);
    std::cout << "in double-double arithmetic, 0.1 * 10 is representable as "
              << high << " + " << low << "\n\n";
 
    // error handling 
    std::feclearexcept(FE_ALL_EXCEPT);
    std::cout << "fma(+Inf, 10, -Inf) = " << std::fma(INFINITY, 10, -INFINITY) << '\n';
    if (std::fetestexcept(FE_INVALID))
        std::cout << "    FE_INVALID raised\n";
}

0.1 double is 0.10000000000000000555112 (0x1.999999999999ap-4)
0.1*10 is 1.0000000000000000555112 (0x8.0000000000002p-3), or 1.0 if rounded to double
0.1 * 10 - 1 = 0 : 1 subtracted after intermediate rounding
fma(0.1, 10, -1) = 5.55112e-17 (0x1p-54)
 
in double-double arithmetic, 0.1 * 10 is representable as 1 + 5.55112e-17
 
fma(+Inf, 10, -Inf) = -nan
    FE_INVALID raised

See also