Fused Multiply-Add, Single Rounding, and the Cost of Contraction

The expression a * b + c names two operations, and on most processors built before the mid-2000s it compiled to two instructions with a rounding after each. A fused multiply-add computes the same expression with the multiplication performed exactly and a single rounding applied to the final sum. The difference is one rounding out of two, which sounds like a minor economy until the expression sits inside a loop that runs a million times, and until the same source compiled for two machines starts producing different last bits.

Two roundings against one

Without fusion, a * b + c is evaluated as fl(fl(ab) + c). The product is rounded to p bits, and then the sum is rounded. The error has two contributions, and the standard bound is

$$ |\operatorname{fl}(\operatorname{fl}(ab) + c) - (ab + c)| \le u|ab| + u|ab + c| $$

which is bounded by roughly 2u · max(|ab|, |ab + c|). With fusion the product is not rounded at all: the exact value ab is added to c and the sum is rounded once, giving the standard model |δ| ≤ u relative to the exact ab + c.

The improvement is a factor of two in the worst case, and it is not merely a constant. In a dot product of length n, the unfused version rounds each product and each partial sum, so the accumulated bound is proportional to n·u; the fused version removes one of the two roundings per term. For a Horner evaluation of a polynomial, fusion halves the number of roundings per coefficient. Neither change alters the algorithm, the memory traffic, or the loop structure.

There is a caveat that belongs with any worst-case statement: a fused result is never worse than the unfused bound, but for a particular input the unfused computation can happen to land closer to the exact value. Fusion improves the guarantee, not every individual result.

What the unit computes

The arithmetic is not difficult to describe, and describing it explains why the operation costs what it does. A binary64 FMA takes two 53-bit significands, multiplies them to a 106-bit exact product, aligns that product against the 53-bit addend, adds, normalises, and rounds once.

The alignment is the expensive part. The addend can be many binades larger or smaller than the product, so the adder has to shift one operand by up to the width of the combined range. Implementations use a wide adder — wider than the 106-bit product, with extra bits for the carry and for the sticky information — and a normalization shifter after it. That is more area than a multiplier and an adder separately, which is why early processors did not include the operation and why some low-power cores still implement it in multiple instructions.

Two edge cases are handled by the same machinery. If the exact ab + c is tiny relative to the operands, the result may be subnormal, and the single rounding is applied at the subnormal spacing rather than at the normal precision. If the exact result is exactly zero, the sign follows the addition rules rather than the product’s sign, so fma(a, b, −ab) returns +0 in round-to-nearest even when ab is negative.

The operation is required by IEEE 754, under the name fusedMultiplyAdd, and C exposes it as fma in <math.h>. A call to fma is required to produce the correctly rounded fused result whether or not the hardware has the instruction; on a target without it, the library routine computes the same value more slowly, typically by splitting the operands and combining partial products.

Contraction, and why the same source gives two answers

A compiler is permitted to replace a separately rounded multiply followed by an add with a single fused instruction. The transformation is called contraction, and it is a source of last-bit differences that are invisible in the source code.

The C standard provides a pragma to control it: #pragma STDC FP_CONTRACT ON, OFF, or DEFAULT, scoped to the translation unit or to a block. The default is not fixed by the standard, and it differs between compilers and between standard-conformance modes; a compiler may also contract only within a single expression, or across statements when a more aggressive mode is selected. GCC and Clang expose the choice as -ffp-contract=off, on, or fast, where the last permits contraction across statement boundaries.

The practical consequence is that a numerical routine can produce different bits when compiled for a machine with FMA, for a machine without it, and for the same machine with a different optimisation level. Since contraction is an accuracy improvement, the difference usually appears in the last bit or two, which is exactly the difference that a regression test comparing bit patterns will catch and a tolerance-based test will not.

Contraction also interacts badly with code that depends on the rounding of an intermediate product. Dekker’s splitting technique, used to compute an exact product from two rounded ones, assumes that the product is rounded in a particular way and that the arithmetic obeys the format’s precision. A contracted expression inside such a routine can break the invariant that the algorithm relies on. The modern replacement, described below, uses the fusion deliberately rather than fighting it.

Recovering the exact product

Fusion makes the exact residual of a product available in two instructions, which is the foundation of the extended-precision techniques in Error-Free Transformations and Double-Double Arithmetic.

#include <math.h>

/* Exact product: p + e == a * b, with p the rounded product and e the residual.
   Requires a fused multiply-add. */
static void two_product(double a, double b, double *p, double *e)
{
    *p = a * b;
    *e = fma(a, b, -*p);   /* exact: a*b - fl(a*b) is representable */
}

The residual e is exactly representable, and the reason is the same one that makes the fused operation useful: a*b − fl(a*b) is the rounding error of a single correctly rounded operation, and that error is always expressible in the same format. Without FMA the same result needs the operands split into high and low halves, several multiplications, and a chain of exact subtractions — more instructions and more assumptions about the rounding mode.

A dot product can then be accumulated as a sum of exact products, with each residual folded into a compensation term. The result is a dot product whose error bound does not grow with n at first order, at roughly twice the arithmetic of the naive loop.

/* Dot product with a single compensation term. More accurate than the
   naive loop, less expensive than a full expansion. */
double dot_compensated(const double *x, const double *y, size_t n)
{
    double s = 0.0, c = 0.0;
    for (size_t i = 0; i < n; i++) {
        double p = x[i] * y[i];
        double e = fma(x[i], y[i], -p);   /* exact product residual */
        double t = s + p;
        c += (s - t) + p + e;             /* both residuals, accumulated */
        s = t;
    }
    return s + c;
}

The loop is serial: s and c are carried across iterations, so the compiler cannot vectorise it without changing the association. That is the same trade-off that compensated summation makes, and it is the reason high-performance libraries use blocked or per-lane variants that keep the compensation while exposing parallelism.

Where fusion changes more than the last bit

Three situations turn the choice from a rounding detail into a design decision.

The first is a code path that must produce identical results across machines. If a regression suite compares output files byte for byte, contraction has to be disabled, or the build has to guarantee FMA availability, and the choice has to be recorded with the results. Disabling contraction costs accuracy and, on some targets, throughput; enabling it costs reproducibility across a heterogeneous fleet.

The second is any algorithm that extracts the rounding error of a product. Dekker’s splitting and its descendants assume unfused semantics; code that uses them and is compiled with contraction enabled can violate its own invariants. The fix is to use the FMA-based formulation, which is both faster and immune to the transformation.

The third is the interaction with reassociation. Contraction alone preserves the expression tree; a compiler that also reassociates can turn a sum of products into a different sum of products, and the two transformations together can move a result much further than either alone. The distinction between “the compiler fused an operation” and “the compiler changed the arithmetic” is the difference between a one-bit difference and an unbounded one, and it is the subject of Reassociation, Excess Precision, and the Limits of Bitwise Reproducibility.

None of this makes fusion a hazard to be avoided. It is a strictly more accurate primitive, available on every mainstream processor since the 2010s, and the reason dense linear algebra routines reach their accuracy bounds. The hazard is that it is applied by the compiler rather than requested by the source, so the arithmetic a program performs is decided partly by build flags that are rarely read alongside the code.