Floating-point arithmetic: an audit of what every number here runs on

Every number I publish is computed in the same arithmetic: 64-bit floating point, the format nearly every computer uses for decimals. It holds about sixteen significant digits and rounds everything else away. I had used it on every project and never once checked what it actually guarantees.

On August 30, 2026 I audited it from the bits up. The ground truth throughout was exact fractions and whole numbers of any size, which Python can do with no rounding at all. I never used a wider floating-point format as the reference, because a wider guess is still a guess.

The rounding error is not lost

When a computer adds two numbers and rounds, it looks as if the lost part is gone. It is not. The error of that one addition is itself exactly representable as a 64-bit number, and six more additions recover it exactly.

I checked that on ten million pairs, from the smallest numbers the format can hold to 10³⁰⁰, against exact arithmetic. There were no exceptions. The error of a multiplication is recoverable too, and two separate methods for recovering it agreed bit for bit on a million pairs.

This is a tool, not a curiosity. Pairing each number with its own rounding error gives about 32 significant digits with ordinary hardware. The library I built that way landed within 6 parts in 10³³ of the exact answer.

Cancellation decides almost everything

Adding numbers that mostly cancel is the hard case. The measure of it is the condition number: the sum of the sizes divided by the size of the sum. A condition number of 10¹⁰ means the answer is ten billion times smaller than the numbers that made it.

Plain left-to-right addition loses accuracy in proportion to the condition number. So does the pairwise method numpy uses, a constant factor better. So does Kahan’s compensated summation, which is the method most people reach for. It was only a median 6.5 times better than plain addition. I had predicted it would stay accurate up to enormous condition numbers. It held only to about 100.

Two methods stayed flat at the limit of the format across fourteen orders of magnitude of cancellation. One is Neumaier’s refinement of Kahan’s method. The other, math.fsum, returns the exactly rounded answer and ships with Python.

The random numbers live on a friendlier grid

One test came back with exactly zero disagreements over a million triples of random numbers, which looked like a bug. It was the grid.

Python’s random.uniform(-1, 1) returns whole-number multiples of 2^-52. Any two of them add to another number on that grid, so every sum is exact. Its random.random() gives exact sums three quarters of the time. Ordinary 64-bit numbers in the same range manage half. I derived each fraction from the bit patterns first and then measured them on two million pairs.

So anyone testing floating-point behaviour with Python’s random numbers is testing a kinder arithmetic than the real one.

One extra bit, thirty-eight correct digits

Rump’s example is a short expression, named for the mathematician Siegfried Rump, that floating point gets badly wrong. With its two inputs, 77617 and 33096, the exact answer is minus 54767/66192, about minus 0.827. I computed that with fractions.

64-bit floating point gives about minus 1.18 × 10²¹, wrong by twenty-one orders of magnitude. Twice the precision gives 1.1726, positive and still wrong. That 1.1726 is the value the example is famous for. Rump’s 1988 paper reports it, from all three precisions of an IBM System/370 mainframe, whose floating point is not today’s standard. In my runs it came out of every precision from 64 to 121 bits, and never out of standard 64-bit floating point.

Then the answer falls off a cliff. At 121 bits it has no correct digits. At 122 it has thirty-eight.

I had an explanation ready before I checked: the condition number of the expression is 2^124.85, two bits from the cliff. It was a coincidence. The real reason is that the inputs are whole numbers, so every term in the expression is exactly representable once the precision is wide enough. The largest term needs exactly 122 bits. At 122 bits the whole calculation becomes exact.

Nudging the inputs off whole numbers separates the two explanations. Then the standard error bound is tight, and the answer improves gradually rather than all at once. Past the crossover both versions gain 0.301 correct digits per extra bit, which is the theoretical log₁₀2 to three places.

Agreement is not evidence

Here is the finding I use most. Two different formulas for the variance of a dataset, run on data with heavy cancellation, gave the same answer to the last bit in all 400 datasets. Both answers were wrong, by about one part in a hundred thousand.

You cannot see that by comparing the two answers. You can see it by comparing each one’s error against the exact value, across many datasets, and asking whether the errors move together. When two methods share a step that throws precision away, their errors are almost perfectly correlated. When they share nothing, the correlation sits near zero.

Writing this up, I found the experiment was not reproducible. It seeded its random data from Python’s string hash, which Python changes on every run, so each run drew different datasets. The −0.026 I reported that night for two unrelated methods was one draw I could never get back. I fixed the seed and ran it at 20 different seeds. The two methods that share a step correlate at 0.996 to 0.998 every time. Methods that share nothing land between −0.10 and +0.12. The finding holds in every draw. The single number I printed did not mean much on its own.

Where I was wrong

Eight of nineteen predictions lost.

Two of my headline results were properties of my test data, not of arithmetic. I had built cancellation out of exact pairs of equal and opposite numbers. That made the true sum always representable, so every exact method scored a free zero. It also meant sorting the numbers placed each pair side by side, where they cancel perfectly, which is why one sorting result won by a factor of 34 million. With data that does not cancel exactly, the same comparison showed a 1.7% difference.

One prediction’s pass condition was that the answer has the right sign. The answer has the right sign at 24 bits while being wrong by thirty orders of magnitude. A test a wrong answer can pass is not a test.

My first scan for a curve’s crossover started at 100 bits and reported 100. The crossover is at 40. The search range had chosen the answer.

A text search that my tool did not support returned zero matches instead of an error, and zero was the answer I expected. The real count was 395.

And the figures had problems I only saw by looking. One axis label is meant to read sum|x| / |sum x|. The plotting font lacked the vertical bar, so it read “sum x / sum x”, which is always 1. Since that night the font has gained the bar, and the figure shows it. Other figures had curves running through their titles.

Re-run for this write-up, every program reproduced its saved results apart from timings and the unseeded experiment above.