Floating-point multiplication is not associative

My research code had one physics formula written out by hand in seven places. Not seven related formulas. One formula, a pairwise interaction between particles, written seven ways: in the routines that build a system, in the GPU kernels that simulate it, and in two slow reference implementations whose whole job is to check the fast ones. During a GPU kernel refactoring campaign I decided to collapse the copies. Write the physics once, call it everywhere. A textbook deduplication.

I set the acceptance gate deliberately high: every stored matrix and every computed energy had to be bit-identical before and after the change. Not close. The same bits.

Under that gate, “extract the formula into one function” is impossible as stated. Not hard, impossible. Floating-point multiplication is not associative, and each of the seven copies multiplied its pieces in its own order.

Rounding happens after every operation

On paper, multiplication is associative: (a × b) × c and a × (b × c) are the same number, always. Floating-point numbers follow IEEE 754, the arithmetic standard almost every CPU and GPU implements, and it defines each operation as: compute the exact result, then round it to the nearest number the format can store. Round after every single operation. So (a*b)*c rounds a*b first, a*(b*c) rounds b*c first, and the two intermediate roundings need not land the same way.

In Julia, with ordinary 64-bit floats:

(0.1 * 0.2) * 0.3   # 0.006000000000000001
0.1 * (0.2 * 0.3)   # 0.006

Neither answer is the “true” 0.006, because 0.006 has no exact binary representation. Both are within a rounding whisker of the truth. They are simply not each other. In floating point, the grouping is part of the result.

Algebraically equal, bitwise incompatible

Before touching anything I audited the seven copies. They were algebraically identical, as hoped. Almost none of them were bit-compatible with each other. The disagreements came in three shapes.

The chains grouped differently. Each site’s string of multiplications had its own bracketing, usually by accident of how it was typed years apart. That is the class of the demo above.

The prefactors were factored differently. One site computed a prefactor from two particle diameters as (di*dj)^3. Others multiplied precomputed cubes. Identical in the algebra you learned at school. In single precision:

di = 1.1f0; dj = 0.9f0   # Float32
(di * dj)^3     # 0.97029907
di^3 * dj^3     # 0.97029895

Five parts in 10^8 apart, a step or two in the grid of numbers a Float32 can store (each step is called a ULP, a unit in the last place). That is the entire disagreement. It is also infinitely more than zero, which is what a bit-identical gate permits.

The projections divided at different moments. Some sites kept stored unit vectors and dotted against those. Others dotted against the raw displacement and divided by the distance afterwards. With a a stored unit vector and d a displacement, both Float32:

a = (-0.25539315f0, 0.29757085f0, -0.91990536f0)   # unit vector
d = (0.70521975f0, -0.1381614f0, -1.0842874f0)
# invr = 1/|d|, the inverse length of d

(a[1]*d[1] + a[2]*d[2] + a[3]*d[3]) * invr               # 0.59672177
a[1]*(d[1]*invr) + a[2]*(d[2]*invr) + a[3]*(d[3]*invr)   # 0.5967217

The whole audit compresses to three rows:

Expression pair	The two values	Where the class bites
(0.1*0.2)*0.3 vs 0.1*(0.2*0.3)	0.006000000000000001 vs 0.006	any chain of multiplications, bracketed by typing habit
(di*dj)^3 vs di^3*dj^3	0.97029907 vs 0.97029895	prefactors built from particle properties
dot then scale vs scale then dot	0.59672177 vs 0.5967217	projections divided by a distance

Five tries

That last pair deserves its own section because of how I found it. I drew random vectors and evaluated both expressions, expecting an instant counterexample. The first pair agreed, bit for bit. So did the second. And the third, and the fourth. The difference above is try number five.

Four clean agreements in a row is a spot check’s idea of proof. If I had eyeballed a couple of cases and moved on, which is exactly what a spot check is, I would have merged the two spellings and quietly changed the output of some small fraction of the billions of formula evaluations a long simulation performs. Most inputs agree. That is what makes this bug class so insidious: the reassurance is part of the trap.

Why the gate demands identical bits

A tolerance gate would have been far easier to pass. That is the argument against it. For a pure refactor, bit-identity is the cheapest absolute proof there is: if every output is the same bits, the change altered nothing observable, and no further discussion is needed. “Within 1e-12” opens a negotiation about how close is close enough, and that negotiation has no natural end.

There is a physics reason too. A long Monte Carlo simulation is chaotic in the practical sense: a last-bit difference in one intermediate value eventually flips one accept-or-reject decision, and from there the two runs part ways completely. The new run is not wrong physics. It is a perfectly valid alternative history of the same system. But it can no longer be compared against the pre-change reference, and that comparison was the whole point of the exercise. Bit-identity is also what lets you re-run last year’s simulation and get last year’s numbers back, exactly, which matters more with every chapter that depends on them.

Prove the gate before trusting it

One cheap trick I will now use every time: before changing any code, I ran the reference capture twice and compared it with itself. Same bits. That one boring pass proves the gate itself is deterministic, with no hidden run-to-run wobble, so any mismatch it reports later can only have come from my change. It cost minutes, and it upgrades “the gate passed” from a hopeful sentence to a meaningful one. I have been burnt before by a check that passed without proving anything; this is the vaccine for one strain of that disease.

One source of truth with two doors

The refactor that landed is smaller than the one I planned, and better. The physics of the formula, which boils down to a sign and a factor of 3, is now written once, in a tiny combine function. On top of it sit two thin arithmetic entry points, one per projection convention: one door for call sites that keep stored unit vectors, one for call sites that scale after the dot. And every call site keeps its own prefactor grouping and its own guard logic, frozen exactly as it was.

So “one function” became “one source of truth with two doors”. The gate held. Every matrix, every energy, bit-identical.

One deliberate exception: the two slow reference implementations were not rewired through the shared function, though that would have been the tidy thing to do. A validator that calls the code it is validating can never catch a bug in that code; agreement between a function and itself is not evidence. They stay textually independent, and their agreement with the shared arithmetic is pinned by tests instead.

Your parentheses are a contract, until fast-math

None of this would be controllable if compilers were free to reorder multiplications. Mostly, they are not. Because IEEE 754 makes reassociation an observable change, an ordinary compile preserves your source-level grouping: the brackets you write are the brackets you get. That is what made “freeze each call site’s grouping” a real engineering option rather than wishful thinking.

Fast-math is the flag that dissolves the contract. GCC and Clang’s -ffast-math, Julia’s @fastmath, and their relatives license the compiler to pretend floating-point algebra is real algebra and reassociate as it pleases. Sometimes that is a good trade for speed. But it is a trade, and you should know whether anything in your stack has made it on your behalf before you promise anyone bit-identical results.

Next time you deduplicate a formula

When a refactor of numerical code promises “no behaviour change”, decide out loud which promise that is: equal within tolerance, or equal to the bit. They are different projects. If the promise is bit-equal, then multiplication order is not a style choice. It is behaviour, and it is untouchable.

And before merging two spellings of the same formula, run them side by side over random inputs and compare the bits, not the values. == cheerfully calls -0.0 and 0.0 equal and refuses to call a NaN equal to itself, so in Julia use === or compare the raw bit patterns. Expect to need many tries before the first disagreement. The agreement that comes first is not evidence of safety; mine agreed four times before it didn’t.

This exact gotcha had come up in a conversation with my supervisor months ago, and I nodded along, because I knew it. Knowing it did not shorten the audit by a single minute. Some facts you only believe after they have cost you an afternoon.