One multiply the compiler removed

Everything in this field so far has needed a machine with more than one worker. This does not. It needs one processor, one thread, one core, and a compiler that was allowed to make a substitution every language standard permits.

a*b + c is two operations. A processor with a fused multiply-add computes it as one, forming the product to full width and rounding the sum once; a processor without one rounds the product first and the sum second. Both are conforming — IEEE-754 specifies what each operation returns, and says nothing about which operations a compiler emits — and C and Fortran both explicitly permit the contraction. No flag in the source records which one happened.

The difference is one rounding of one product: a relative 10⁻¹⁶. Here is what it decides.

From k = 27 the expression as written returns zero. Not a small number: zero, which has no sign, and the sign is the only thing anybody wants from a 2 × 2 determinant.

Extending the sweep of sizes shows how much of the range that zero occupies.

The failure point never moves and the failed region grows to fill the range. At ceilings of 28, 36, 40 and 45 the naive form is exact to k = 26 in every one — the same integer, not a trend — while the count of sizes returning exactly zero runs 2, 10, 14 and 19. Every size past 26 fails, so the zero region is not a patch of bad luck at particular k: it is everything above a threshold fixed by the format.

The fused column is the control, and it is exact at all four ceilings. That is what makes this a statement about the contraction rather than about the magnitude of the entries: the same entries, the same determinant, the same precision, and one form is right at every size while the other is right up to 26 and identically zero after. A difference of one rounding in one product, followed to its conclusion, is the difference between an exact answer at k = 45 and no answer at all.

Where the zero comes from #

The matrix is chosen so that its two products are consecutive integers. With x = 2ᵏ,

ad = (x+2)(x+2) = x² + 4x + 4,
bc = (x+1)(x+3) = x² + 4x + 3,

so the determinant is 1 exactly, at every k, and the two products differ by one part in x². Binary64 carries 53 significand bits, so once x² passes 2⁵³ — which is k = 27 — the spacing of the representable numbers near x² is larger than 1 and the two products round to the same float. Their difference is zero.

Nothing has gone wrong in any operation. Each product is correctly rounded, the subtraction of two equal numbers is exact, and the answer is zero. This is cancellation takes the answer in its purest form: the subtraction introduces no error at all and exposes error the operands were already carrying, and here the operands were carrying all of it.

What the fusion recovers #

Kahan’s algorithm for the same determinant is three lines and needs the fusion:

w = b·c                 rounded
e = fma(b, c, −w)       the exact residual of that rounding
f = fma(a, d, −w)       ad − w, with one rounding
det = f − e

The second line is where the fusion earns its keep. fma(b, c, −w) computes b·c − w with a single rounding, and since w is the rounded product, the true value of b·c − w is exactly representable — so e is the exact error of the rounded product, recovered by an instruction. Without the fusion that quantity needs Dekker’s splitting, six operations and a 26-bit split of each operand.

Then f is ad − w computed to one rounding, which is small and therefore exact, and the answer is f − e. The measured error is zero at every k in the sweep, not merely small.

That the sign of e is subtracted rather than added is the whole algorithm, and it is worth recording that getting it wrong returns 7 where the answer is 1 — and does so only in the regime where the naive route has already failed, so a test that checked the easy sizes would have passed.

Where this bites in practice #

A 2 × 2 determinant is not a corner of the subject. It is:

An orientation test. Whether three points turn left or right is the sign of exactly this determinant, and computational geometry has an entire literature about it — exact predicates, adaptive precision, filters — for the reason this figure shows. Numerical linear algebra mostly does not, and has the same expression in a dozen places.

A 2 × 2 pivot. The symmetric indefinite factorisations pick between a 1 × 1 and a 2 × 2 pivot, and the choice reads a determinant of exactly this shape. When symmetry is not enough is where this site measures that decision; the number it reads is one this essay says two builds may disagree about.

A Givens rotation’s discriminant, a discriminant of a quadratic, a cross product, an area. Every one of them is ad − bc.

And a definiteness test. A symmetric 2 × 2 is positive definite exactly when a > 0 and ad − b² > 0, so a determinant computed as zero turns a definite matrix into a semi-definite one. A matrix that is definite on one machine takes that all the way, with a census rather than an example.

Which of the two is right #

An awkward question, and the answer is not the one the figure suggests.

On this family the fused form is exactly right at every size and the unfused form is exactly wrong, so the temptation is to conclude that fusing is the correct implementation and the contraction is a hazard to be turned off. That conclusion is wrong, and the measurement that shows it is one field essay further on: over two hundred Gram matrices whose definiteness is settled exactly in BigInt rationals, the two forms disagree twenty-six times, and the fused form is the one that is right nine of those times. The unfused form is right seventeen.

The difference is that this essay’s determinant is computed by an algorithm designed around the fusion. Kahan’s three lines use the fma to recover an exact quantity, which is a use of the instruction rather than an accident of compilation. Take the same instruction and apply it to an ordinary expression nobody designed for it, and it moves the answer by one rounding in a direction that has no reason to be the right one.

So the two claims to keep apart:

• an algorithm written for the fusion is better than the expression it replaces, measurably and by a lot — this figure, where one is exact and the other returns zero;
• an expression the compiler happened to contract is not better or worse, it is different, and which of the two lands closer to the truth is a property of the data.

The first is a reason to write fma explicitly. The second is a reason not to imagine that a compiler flag is a correctness setting.

Why no residual can see it #

The site’s standing instrument for is this number trustworthy is a residual: compute the answer, put it back into the problem, and measure what is left. It does not work here, and the reason is worth stating because it is the reason this failure is silent.

A residual is computed from the same expression the answer was. If the compiled form contracted the multiply and the add, it contracts them in the residual too — so the residual is computed in the arithmetic whose behaviour is in question, and it reports agreement with itself. The refusal this essay publishes is that statement mechanised: an assertion that a residual computed the same way detects the contraction is fed the case and must reject.

This is a genuinely different failure from the ones the site’s other eighteen fields handle. There, a residual is small and the answer is wrong because the problem was sensitive, and the two measurements together account for it — the identity three errors and one number sets out. Here the residual is not merely uninformative; it is computed by the thing being asked about.

The instruments that do work are the ones that share no arithmetic with the computation: an exact determinant in BigInt rationals, which is what this site’s exact.js provides and what the definiteness census uses for its verdict; or an interval computed in directed rounding, which is proving the answer is in the box’s machinery. Both are second routes in the site’s usual sense, and both cost more than the computation they are checking.

How large the difference between the two forms is #

Bounded, and small, which is what makes the determinant result surprising rather than obvious.

Over ten thousand random arguments, the fused and unfused forms of a*b + c differ by at most 1.36·10⁻¹⁶ relative to |ab| + |a*b+c| — 0.61 of a unit roundoff, which is one rounding of the product and one of the sum, as it must be.

That last clause names the bound and the number in front of it is not the bound; it is the largest value ten thousand draws happened to produce, and a maximum over a sample can only grow. Extending the sweep:

   draws          largest gap
   1,000            0.619 u
   10,000           0.659 u
   100,000          0.659 u
   1,000,000        0.663 u
   4,000,000        0.666 u

Still creeping at four million, and heading for about two thirds of a unit roundoff rather than settling at 0.61. The bound is one unit roundoff, which is what “one rounding of the product and one of the sum” says: the gap is p·δ₁ + (p + c)(δ₂ − δ₃) over |p| + |p + c|, each δ at most u/2, so u is the ceiling — and the sampled maximum is two thirds of the way to it and not obviously going further.

There is one case sampling cannot find, because the gap there is exactly zero. When c is negligible beside a·b, both forms return fl(ab) — the unfused one because the addition changes nothing, the fused one because a·b + c rounds to the same float. Two hundred thousand draws with c at 10⁻⁴⁰ give a maximum gap of 0, and no amount of sampling in the general case would have suggested that the two forms have a region of exact agreement.

None of that changes the essay’s conclusion, which is about a difference of one rounding deciding a sign. It changes what the number is: the gap is bounded by u, is two thirds of u on four million draws, and is exactly zero on a set the sampling never visits.

The fourth time this collection has quoted a sample maximum #

It is worth naming rather than fixing quietly, because the same mistake has now appeared in four essays on this site with a different subject each time: four ratios drawn from seven runs, a count of distinct answers over four hundred permutations, a ratio held fixed along the one axis that cannot move it, and this.

The shape is always the same. A quantity is sampled, its extreme value is reported, and the sentence around it turns the extreme into a property — there is no case above it, the ratio is 1.02, compensation is three times better. In every instance the underlying claim survives and its number does not, and in every instance the fix is one line: report the statistic the sample supports, or sample until it stops moving, or derive the bound.

Here the bound was available and was written in the same sentence. u, from one rounding of the product and one of the sum, is a proof; 0.61·u is an observation about ten thousand draws. The essay had both and quoted the second, which is the general form of this mistake — the derivation and the measurement are two routes to one number, and this site’s own habit is to run both and compare them.

So the entire content of this essay is a difference of one rounding, and it decides a sign. The mechanism is the one the whole field runs on: a rounding committed at the scale of the products, followed by a subtraction that removes everything except the rounding. The products are 10¹⁶ and the answer is 1, so the summation condition number of this two-term sum is 3.6·10¹⁶ — the second highest entry in the vector that hides it’s census, in a sum with two terms in it.

The same instruction, in a longer computation #

A determinant is two products, so the contraction’s effect is a single rounding and the analysis is exact. What happens when the expression is longer is worth measuring rather than extrapolating, because both the size and the character of the difference change.

Take Horner’s rule, which is s = s*x + c in a loop and is therefore the contraction’s natural habitat — every step of it is a multiply-add. Evaluate (x − 1)⁶, expanded into x⁶ − 6x⁵ + 15x⁴ − 20x³ + 15x² − 6x + 1, at 401 points within 10⁻³ of x = 1. The true value is a sixth power and is never negative; it is under 10⁻¹⁸ across the whole window, which is below the rounding of the terms that produced it.

The unfused evaluation returns a negative value at 179 of the 401 points. The fused one returns a negative value at 196. And the two disagree about the sign at 98 — a quarter of the window.

The largest difference between the two evaluations anywhere in that window is 1.44·10⁻¹⁵: six roundings rather than one, exactly as the bound says. What has changed is not the size but the consequence, because there is no longer a single answer that is right — both traces are made entirely of rounding, and neither is more nearly the sixth power than the other. A square that evaluates negative is that measurement, and it is where the field stops being able to say which build is correct.

What a caller can do #

Four options, in ascending order of cost.

Use the fused form deliberately. Kahan’s three lines, with an explicit fma intrinsic rather than a a*d - b*c the compiler may or may not contract. The error is under two ulps of the true determinant whatever the cancellation, which is the strongest statement available about any expression on this site. This is free and it is the answer for a 2 × 2.

Forbid the contraction. -ffp-contract=off and its equivalents make the compiled form predictable. That buys reproducibility across builds and buys the less accurate of the two forms, which is the trade this field keeps producing.

Filter. Compute the naive determinant along with a bound on its error; if the bound does not exclude zero, fall back to something exact. This is what computational geometry’s adaptive predicates do and it is nearly free in the common case, because the common case is a determinant that is not near zero and the bound settles it in two extra operations. The structure is the one an estimate that can be fooled describes for condition estimation, with the sign reversed: there a cheap estimate can be wrong in the flattering direction, here a cheap bound is only ever allowed to return unknown.

Compute exactly. Two products in a twoProduct expansion and an exact sum of four terms, which is what this essay’s ground truth is. About twelve operations for a 2 × 2 and it settles the question completely — and the expansion is the same object the reproducible summation policies of the sum that cannot be wrong are built from, so a code base that has one has the other.

The point of listing them is that all four exist, none is difficult, and the default — write a*d - b*c and let the compiler decide — is the one option that gives an answer whose accuracy depends on a flag nobody set.

What the sweep says about the format #

One reading of the hero figure that is easy to miss: the failure point is not a property of the matrix, and moving to a wider float does not remove it.

The naive determinant fails at k = 27 because 2²ᵏ crosses 2⁵³. In binary32, with 24 significand bits, the same family fails at k = 12; in a hypothetical 113-bit format it fails at k = 57. The failure moves and it does not go away, because the family is parameterised by exactly the quantity the format is: how many bits separate the products from their difference.

Every real problem has a fixed k, of course, and a wider float may put it on the safe side. That is the ordinary and correct use of precision and this site has a whole knob for it — precision as the knob is the arithmetic field’s founding observation. What the figure adds is that the shape is unchanged: there is a threshold, it is set by the format, and on either side of it the naive expression is either exactly right or exactly wrong, with nothing in between. A determinant near that threshold is not approximately computable. It is computable or it is zero.

And the fused route is exact on both sides of it, at every size in the sweep, which is why the recommendation is an algorithm rather than a format.

One line #

A compiler is allowed to replace two roundings with one, from source that says nothing about it, and on an expression whose subtraction removes everything except the rounding that is the whole of the answer.