Scientific Computing · Floating Point

Subnormals, Signed Zero, and the Cost of Gradual Underflow

Why the exponent field needed two reserved encodings, what a subnormal costs in precision, and how signed zero and NaN payloads change the meaning of equality.

The normal-number encoding described by the sign, exponent, and fraction fields covers a range that stops abruptly at 2^emin. Below that point the fraction field has nothing left to give: the exponent cannot go lower, and the significand is already at its minimum of 1.0. The obvious response is to declare everything smaller to be zero. IEEE 754 does not do that, and the three encodings it invented instead — subnormals, signed zero, and NaN — are the source of most of the surprising behaviour in floating-point programs.

The gap below the smallest normal

For binary32 the smallest positive normal is 2^-126, about 1.17549435e-38. A subtraction such as 1.5e-38 - 1.4e-38 has an exact result of 1e-39, which is smaller than that. Without a special encoding the only options are to return zero, losing the entire result, or to raise an exception and stop.

Gradual underflow keeps the value by reinterpreting the exponent field’s all-zeros pattern. Instead of meaning “exponent emin with a leading 1”, it means “exponent emin with a leading 0”, so the significand becomes f/2^t in [0, 1). The spacing in this range is fixed at 2^(emin−t): 2^-149 for binary32, about 1.4e-45, and 2^-1074 for binary64, about 4.9e-324. The smallest positive subnormal is exactly that spacing, and the largest is (1 − 2^-t) × 2^emin, immediately below the smallest normal.

The range below 2^emin therefore holds 2^t − 1 values per sign — 8,388,607 of them for binary32 — and they are evenly spaced rather than logarithmically spaced. That change of spacing is the whole cost of the mechanism. A normal number’s relative precision is constant at p bits; a subnormal’s absolute precision is constant at 2^(emin−t) while its relative precision collapses as the value shrinks. The smallest positive subnormal has a rounding error of up to half its own value, so its relative error bound is 0.5 rather than 2^-24.

What gradual underflow preserves

The design question is whether that loss is worth the trouble compared with flushing to zero. Kahan’s argument, made in the context of complex elementary functions, is that abrupt underflow breaks algebraic identities that programs depend on, and it breaks them silently.

Consider x = (a * b) / c with a and b small and c smaller still. The product a * b may fall below 2^emin while the final quotient is perfectly representable. With flush-to-zero, the intermediate becomes zero and the result is zero; the same expression rewritten as a * (b / c) gives a correct answer. Under gradual underflow the product is merely imprecise, and the division recovers most of the magnitude. The two orderings still disagree, but they disagree in the last bits rather than in every bit.

The second property is exactness of subtraction near cancellation. Sterbenz’s lemma says that if y/2 ≤ x ≤ 2y, then x − y is exactly representable. In the subnormal range that remains true, because the difference of two nearby subnormals is a multiple of the subnormal spacing and is small enough to be encoded. Flush-to-zero destroys it: two distinct values in the subnormal range both become zero, and the difference becomes zero, so a test for equality reports them as equal when they are not.

The price is paid by hardware. Subnormal operands and results need the leading-zero count and a variable shift in the significand path, and on several microarchitectures they are handled by a slow assist path rather than by the main arithmetic datapath. The penalty is microarchitecture-specific and can be large, which is why performance-sensitive code sometimes enables flush-to-zero deliberately. On x86 the control is a pair of bits in the SSE control and status register, one for flushing results and one for treating subnormal inputs as zero; compilers expose them through flags rather than through a portable interface, and the aggressive optimisation levels enable them as a side effect.

Two zeros

The sign bit applies to zero as well, so the encoding has +0 and −0. They compare equal under ==, and they behave identically in addition and subtraction of nonzero values. They differ wherever the sign of a zero carries information.

The most visible case is division: 1.0 / +0.0 is +∞ and 1.0 / −0.0 is −∞, both raising the divide-by-zero exception. The sign of an infinity produced by overflow comes from the sign of the operands, and if the overflow happened through an underflowed quantity, the zero’s sign is the only record of which direction the value was heading. The same reasoning applies to log(−0), which is −∞ rather than a NaN, and to atan2, whose result depends on the signs of both arguments and is defined for zero arguments.

Arithmetic on zeros follows rules that look arbitrary until the sign is seen as data. (+0) + (−0) is +0 in round-to-nearest, because the exact sum is zero and the standard picks the positive sign for an exact zero sum. In round-toward-negative the same expression yields −0, since the exact result is representable and the mode prefers the lower neighbour. x − x is +0 in round-to-nearest for the same reason, and −0 under downward rounding. Multiplication and division are simpler: the sign of a zero result is the exclusive-or of the operand signs, with no mode dependence. The reason the sign is chosen for an exact zero sum is the tie rule applied to the infinitely precise result, described in Round-to-Nearest-Even and the Mechanics of Correct Rounding.

A program that uses the sign of zero as a flag is relying on behaviour the standard guarantees but that compilers may not preserve. An expression rewritten as x - x by an optimiser can lose the sign; a comparison x == 0.0 cannot distinguish the two zeros at all, and the portable test is signbit(x), which the C standard defines to work on zeros as well as on other values.

NaN as an encoding, not a value

The all-ones exponent field with a nonzero fraction encodes NaN, and the fraction field is a payload. The most significant payload bit distinguishes quiet from signalling: clear means signalling, which raises the invalid-operation exception when consumed. That leaves 22 payload bits in binary32 and 51 in binary64, so the encoding has 16,777,214 distinct NaN patterns in binary32 counting both signs.

Almost nothing about the payload is portable. The standard requires that an operation on a quiet NaN return a NaN, but it does not require that the payload survive, and implementations differ: some propagate the first NaN operand’s payload through SSE-style instructions, some return a canonical value, and ARM’s floating-point units have a default-NaN mode that discards payloads entirely. The canonical quiet NaN produced by an invalid operation on x86 is the pattern 0xFFC00000 for binary32, sign bit set and only the quiet bit in the payload.

Comparisons are where NaN changes program logic rather than just results. Every ordered comparison involving NaN is false, including NaN == NaN, while the unordered != is true. That asymmetry is the basis of the idiom x != x for detecting NaN without calling a library function, though a signalling NaN operand can raise an exception during the comparison, and compilers under fast-math assumptions are permitted to fold the expression to false.

#include <math.h>
#include <stdio.h>

int main(void)
{
    double pz = 0.0, nz = -0.0;
    printf("%d %d\n", pz == nz, signbit(nz));      /* 1 1 */
    printf("%g %g\n", 1.0 / pz, 1.0 / nz);         /* inf -inf */

    double q = NAN;
    printf("%d %d\n", q == q, q != q);             /* 0 1 */
    printf("%d\n", isnan(q) != 0);                 /* 1 */
    return 0;
}

The standard also defines a total ordering that does distinguish what == conflates. It places −0 before +0, orders NaNs after every number, and orders negative NaNs before positive ones, with payloads compared as integers within a sign. That ordering is what a sort or a hash needs if it is to be a total order on bit patterns, and it is not what < implements.

Where the abstraction leaks

Three boundaries are worth knowing. The first is that the point at which a result is detected as tiny — before or after rounding — is implementation-defined, so the exact pattern of underflow flags raised near 2^emin can differ between machines that agree on every value.

The second is that the standard permits a flush-to-zero mode as an extension, and many compilers enable it under flags intended to relax floating-point semantics. Code that relies on subnormals for correctness, such as a solver that iterates to very small residuals, can fail under those flags while producing no diagnostic.

The third is that these special values interact with the error analysis that the rest of numerical computing rests on. The model fl(x op y) = (x op y)(1 + δ) with |δ| ≤ u holds only in the normal range; in the subnormal range the bound is absolute, and with a NaN operand there is no bound at all because there is no value. Distinguishing those regimes is the first step in any serious error analysis, which is the subject of Forward Error, Backward Error, and Condition Numbers.

References

  1. IEEE 754-2019: Standard for Floating-Point Arithmetic — IEEE
  2. What Every Computer Scientist Should Know About Floating-Point Arithmetic — David Goldberg, ACM Computing Surveys 23(1), 1991
  3. Branch Cuts for Complex Elementary Functions, or Much Ado About Nothing's Sign — W. Kahan, in The State of the Art in Numerical Analysis, Oxford, 1987