The Floor Beneath the Instability

A sidebar from Divide by Zero, by Roberto Mas

On the morning of 19 October 2026, two runs of the same program began from fields that were identical to the last bit. The program simulated the Gierer–Meinhardt equations, Alan Turing's activator and inhibitor in their most studied form, on a grid of 64 by 64 cells. Every cell held the same values, the uniform state, except for a small smooth ripple at the center, 0.000036 high.

The mathematics was not in doubt. A linear stability calculation, the kind Turing did by hand in 1952, said the uniform state was unstable: any ripple of the right width should grow by a factor of e every 26 time units until it became a field of spots.

The first run did what the theory said. By time 244 the spots were there.

The second run differed in one setting. It advanced the clock in steps a quarter as long, which is what a careful modeler does to make a simulation more accurate. At time 1,500 its field was still flat, and the ripple at the center was exactly the height it had started at, to every stored digit. By the theory it should have grown by a factor of about 10^24. Nothing flagged it. The program reported a calm field, in the same format it reported everything else.

Half a step

The program stored each value as a 32-bit float, with 24 binary digits of precision. Near the uniform activator level, 1.525, the gap between one storable number and the next is 2^-23, about 0.00000012. That gap is called a unit in the last place, or ULP.

At each step the program computes the ripple's growth, adds it, and rounds the result to the nearest storable number. The growth per step is the ripple's height times the growth rate times the length of the step. When that product is less than half a ULP, rounding returns the old value. The ripple grows by nothing, and it does so again on the next step.

So there is a floor. Below a certain height, a ripple the theory says must grow cannot grow at all:

δ* = ULP ÷ (2 · σ_max · Δt)

Here σ_max is the fastest growth rate from the stability calculation and Δt is the length of a step. For the first run the floor was 133 ULP. For the second, with steps a quarter as long, it was 530 ULP. The ripple was 300 ULP high, above one floor and below the other.

The formula is a half-step argument, not a theorem, and the measured edge sits a little above it. Across three step lengths, ripples at 1.2 times the formula's floor stayed frozen and ripples at twice the floor grew, at 85 to 91 percent of the predicted rate. The floor moved exactly as the formula said it should, doubling each time the step was halved.

More careful, less hearing

The ordinary instinct runs the wrong way here. Shorter steps cut the error of approximating a smooth change by a staircase, and textbooks recommend them. They also cut the growth added at each step, while the rounding gap stays the same size. Each halving of the step doubles the height a disturbance must reach before the program can register it.

Wider numbers lower the floor without removing it. A 64-bit float has a gap about 540 million times smaller, so the same field in double precision would hear ripples hundreds of millions of times fainter. It would still be deaf below its own floor.

None of this is new in kind. Small increments lost against a large running total are the reason William Kahan published a compensated summation method in 1965, the year before the report this book opened with. What is particular here is the consequence for a system poised to break its symmetry. The arithmetic does not merely blur the answer; below the floor it turns an unstable equilibrium into a stable one. [I have not searched the literature for a statement of this threshold. Until that search is done, this is an observation, not a priority claim. †]

Which machine this is

Of the Coda's three machines, this is the Pentium. It does not refuse. Asked whether the field will form a pattern, it answers no, plainly, and the answer looks like every other answer it gives. No NaN appears, because nothing invalid has happened; every operation was performed exactly as the standard specifies.

The check that caught it was Nicely's. A second run, built differently and asked the same question, disagreed. Here the second machine did not need to be a second chip: the same program with a different step length was enough, because the floor moves with the step while the physics does not. A result that changes when only the step changes was never a result about the equations.

The program can also be made to hear its own floor. It knows σ_max, the step length and the value it is storing, so it can compute δ* and compare it with the height of whatever disturbance it is carrying. When the theory predicts growth and the field sits below the floor, it can say so instead of reporting calm. The simulation this sidebar describes now does that.

The standard, applied

What is unknown. Nothing about chemistry or morphogenesis. What is hidden is whether a given calm field is calm because of the equations or because of the arithmetic.

Which framework makes it unknown. Fixed-width floating point, rounded to nearest, combined with an explicit step that adds growth in small increments. Change either one and the floor moves.

What would change the judgment. A second run at a different step length or a wider format; a disturbance raised above twice the floor; or the program's own comparison of δ* with the field it holds.

Status, 19 October 2026. Reproduced at three step lengths in one implementation, running on one JavaScript engine. Not yet tested on other hardware, with implicit integrators, or with fused multiply-add arithmetic, any of which could move the edge. [The general phenomenon is long known; whether this threshold has been stated for Turing instabilities is unchecked. †]

Source note (for Notes on the Sources)

Runs made on 19 October 2026 with Turing Bench, a browser reaction-diffusion simulator, and its headless twin under Node 22 (V8 engine). Model: Gierer–Meinhardt with production ρ = 0.1, activator decay μ = 0.08, inhibitor decay ν = 0.12, basal source 0.002, no saturation; diffusion 0.08 (activator) and 3.2 (inhibitor); 64 × 64 periodic grid; nine-point Laplacian; explicit Euler; fields stored as 32-bit floats. Linear stability gives σ_max = 0.0379 per time unit and a fastest wavelength of 13.3 cells. Disturbances were Gaussian ripples of width 3 cells added to the uniform state and compared against an unperturbed twin run.

The 300-ULP ripple of the opening: spots by time 244 at Δt = 0.0996; unchanged to the last bit at time 1,500 at Δt = 0.0249. Scripts: exp4.js and exp6.js in the bench workspace.