Reassociation, Excess Precision, and the Limits of Bitwise Reproducibility

Floating-point addition is not associative. That single fact, combined with compilers that are allowed to exploit the algebraic identities of real arithmetic, means the arithmetic a program performs is only partly determined by its source text. Two builds of the same file can produce different last bits, and neither is wrong. Understanding which transformation is responsible is the difference between a reproducibility policy that works and a set of flags copied from a mailing list.

Three transformations, three mechanisms

The first is reassociation. A compiler may evaluate (a + b) + c as a + (b + c) when it is licensed to assume that floating-point addition behaves like real addition. The two orders generally give different results, because each rounds a different partial sum. The licence is not granted by default in a conforming C implementation, but it is granted by any optimisation level that enables the “unsafe” floating-point bundle, and it is granted implicitly to vectorising transformations that need to split a reduction into lanes.

The second is contraction, where a multiply and an add become one fused operation with a single rounding. It is a smaller change than reassociation, since the expression tree is unchanged and one intermediate rounding is removed, but it is applied by the compiler rather than requested by the source, so it changes results invisibly. The pragma and flag controls are described in Fused Multiply-Add, Single Rounding, and the Cost of Contraction.

The third is excess precision: the hardware or the language implementation evaluates an expression in a format wider than the operands, then rounds on assignment. This one is not an optimisation at all. It is a property of the evaluation model, and it can make the result depend on whether a value happened to stay in a register.

Excess precision and the register spill

The x87 floating-point unit, which was the only way to do floating-point on x86 before SSE2, computes in an 80-bit format with a 64-bit significand, a 15-bit exponent, and an explicit integer bit. A C expression involving double operands is evaluated in that wider format, and the result is rounded to binary64 only when it is stored to memory.

The consequence is that the value of an expression depends on register allocation. Consider x = a * b + c where the intermediate a * b stays in a register and is added at 64-bit precision. If the compiler runs out of registers and spills a * b to memory, the store rounds it to 53 bits, and the subsequent addition produces a different result. A change in the surrounding code, or in the compiler’s register allocator, changes the arithmetic of a line that was not touched.

C’s evaluation model has a name for this: FLT_EVAL_METHOD, defined in <float.h>. A value of 0 means expressions are evaluated in the format of their operands; 1 means float and double are evaluated as double; 2 means they are evaluated as long double. The x87 configuration reports 2. On x86-64, where SSE2 handles float and double natively, the value is normally 0, and the excess-precision problem disappears from that platform, which is why code that behaves consistently on x86-64 can still differ on a 32-bit build of the same program.

The language does provide a way to remove the excess: assignment to a variable, a cast, or a function call forces the value to the nominal type. A store to a volatile object is the blunt version, since the compiler cannot keep it in a register:

volatile double rounded = a * b;   /* forces a binary64 store even when the
                                      register holds an 80-bit value */

Compilers also offer a mode that implements the language’s rules exactly (GCC’s -fexcess-precision=standard) at some cost in optimisation freedom, because every assignment then becomes a rounding point that cannot be optimised away.

What the fast-math bundle licenses

The flags collected under names like -ffast-math are not one transformation but a bundle of assumptions, and each assumption changes which programs are equivalent.

Reassociation is the first: (a+b)+c = a+(b+c) holds only if the rounding of the intermediate sums is ignored. The assumption that no NaN or infinity is produced is the second: under it, x == x is always true, comparisons can be simplified, and code that tests for NaN is dead and can be removed. The assumption that the sign of zero does not matter is the third: x + 0.0 can be simplified to x even when x is -0.0, and 0.0 * y can be replaced by 0.0 even when y is negative. Reciprocal substitution is the fourth: a / b can become a * (1/b), adding a rounding.

Each is a real transformation with a real performance motive, and each breaks a program that depends on the corresponding detail. A solver that uses x != x to detect a failed iteration, a geometric routine that relies on the sign of a zero, a compensated summation that depends on an exact residual, all are silently miscompiled by the bundle, and none produces a diagnostic.

The practical rule is that fast-math is a statement about the program’s numerical requirements, not a performance setting. A code that can tolerate it is a code whose results are already approximate to a tolerance much larger than u, and whose control flow does not branch on floating-point special values.

Ordering as the contract

When bitwise reproducibility is required, the only durable contract is a fixed order of operations, and it has to be fixed at every level: the expression, the loop, the vector width, and the parallel decomposition.

At the expression level, contraction and reassociation must be disabled, and the excess-precision mode must be pinned. At the loop level, a reduction must be written so that vectorisation cannot change its association, or the vector width must be part of the build’s identity. A loop that sums an array with s += a[i] and is auto-vectorised into four lanes produces a different result on a machine with eight-wide vectors, because the tree of partial sums has a different shape. The error consequences of that change of shape, and the algorithm that makes the order irrelevant, are covered in Compensated Summation and How Summation Order Bounds Error.

At the parallel level, the reduction order in a thread pool or an MPI collective is usually unspecified. An OpenMP reduction clause leaves the combination order to the implementation, so the result varies with thread count and with scheduling. An MPI MPI_Reduce with a user-defined operation combines in an order the standard does not fix. Neither is a defect; both mean that the reduction is not reproducible unless the program imposes the order itself, for example by reducing into a fixed tree of ranks.

There is a cost to all of this, and it is not small. Disabling vectorisation of reductions can cost a factor of several on memory-bound loops; disabling contraction costs one rounding per term; forcing stores for excess precision costs register pressure. A program that needs bitwise reproducibility pays these costs everywhere, including in the parts that do not need it, unless the strict settings are applied per translation unit.

Reproducibility by exactness

The alternative to fixing the order is to make the result independent of it. An exact accumulator computes the correctly rounded exact sum of its inputs, and the exact sum does not depend on the order in which the terms are added. Two programs that sum the same values with different trees, different vector widths, and different thread counts produce the same bits, because they are all rounding the same exact quantity once at the end.

The mechanism is binned summation. Each input is decomposed by exponent into a fixed set of bins, each bin accumulates into an integer accumulator wide enough to hold the maximum possible contribution, and the bins are summed from the largest down at the end. Because the bins are integer, the accumulation is exact and order-independent; because the number of bins is fixed, the memory cost is bounded and independent of n.

The price is arithmetic per element — several operations to extract the bin index and add — and the technique is most attractive where the reduction dominates the runtime or where the reproducibility requirement is absolute. Libraries implementing reproducible BLAS use this approach to guarantee identical results across thread counts and machines, which is a stronger property than identical results on one machine with a fixed build.

Exactness also settles the cases that ordering cannot. A compensated sum with a fixed order is reproducible but still not exact; a change in the input data changes the result in ways that a tolerance must absorb. An exact sum is reproducible and correctly rounded, so the only remaining variation is in the inputs themselves.

What remains irreproducible

Even with strict flags, a fixed order, and an exact accumulator, three sources of variation remain outside the program’s control.

Library functions are the largest. The C standard does not require sin, exp, or pow to be correctly rounded, and implementations differ in the last bit and sometimes in more. Two platforms running the same strict code can produce different results from the same transcendental call, and a library upgrade can change them. The functions that are correctly rounded on a given platform are a property of that platform’s math library, not of the language.

The second is the hardware’s own guarantees. IEEE 754 requires correctly rounded results for the basic operations, and mainstream processors meet that requirement, but fused operations, reciprocal estimates, and conversions may have implementation-defined behaviour at the edges, particularly in the subnormal range.

The third is the data. A reproducible computation on irreproducible inputs — a file whose values were written with six digits, a sensor whose readings vary, a random seed drawn from the clock — produces different results no matter how strict the arithmetic. Bitwise reproducibility is a property of the function from inputs to outputs; it says nothing about whether the inputs are the same, and a process that needs identical outputs has to pin the inputs as carefully as it pins the compiler.