The first program is stupid enough to look broken:

#include <cstdio>

int main() {
    float x = 16777216.0f;
    float y = 1.0f;
    std::printf("%.1f\n", x + y);
}

Output:

16777216.0

Not 16777217.0.

The number sixteen million seven hundred seventy-two thousand two hundred sixteen survives the addition. The one disappears.

Why is 16777216.0f + 1.0f still 16777216.0f? And why does the same shape of thing happen when summing gradients, computing variance, doing a dot product, reducing across devices, or running a model in lower precision?

This is a numerics post, but the way in is C and C++ because those languages let us look at the actual bits without much ceremony. Assumed background: beginner C or C++, beginner arithmetic intuition, and enough patience to let one impossible-looking print statement pull us down several layers.

Local machine:

Apple M4
Apple clang version 15.0.0 (clang-1500.3.9.4)
Target: arm64-apple-darwin24.6.0
NumPy 2.2.5, OpenBLAS 0.3.29 (Homebrew python3.13)

Unless a block says otherwise, the outputs below come from one probe file, numerics_probe.cpp, which prints a labeled line for every experiment. The inline snippets show the shape of each experiment; the output blocks keep the probe’s labels. The probe source is listed with the rest of the artifacts at the end.

The final thread is an experiment I didn’t run yet: a Metal kernel (Metal is Apple’s GPU API) that accumulates in a chosen precision until it stalls.

Act I: Integers are not simple

1. A byte is eight switches, and 255 plus 1 becomes 0

Start with something that looks more primitive than floats:

#include <cstdint>
#include <cstdio>

int main() {
    uint8_t x = 255;
    uint8_t y = x + uint8_t{1};
    std::printf("%u\n", (unsigned)y);
}

Output from my probe:

uint8_t 255 + 1 -> 0

uint8_t is an unsigned 8-bit integer. Eight bits gives 256 different bit patterns:

00000000 -> 0
00000001 -> 1
...
11111110 -> 254
11111111 -> 255

There’s no 256th value. If unsigned arithmetic is done modulo 256, then:

255 + 1 = 256
256 mod 256 = 0

For unsigned integer types in C and C++, wraparound is defined. The arithmetic is modulo 2^N, where N is the width of the type.

For uint8_t, N is 8:

mod 2^8 = mod 256

For uint32_t, N is 32:

4294967295 + 1 -> 0

For uint64_t, N is 64:

18446744073709551615 + 1 -> 0

Bit manipulation, hashes, checksums, ring buffers, and low-level protocols all lean on unsigned types because the wraparound is part of the contract.

It’s also why bugs can look perfectly legal:

size_t n = 0;
std::printf("%zu\n", n - 1);

Output:

size_t(0) - 1 -> 18446744073709551615

size_t is unsigned on this machine. 0 - 1 wraps to the largest representable size_t.

That one line explains a whole family of bounds-check bugs:

for (size_t i = v.size() - 1; i >= 0; --i) {
    ...
}

If v.size() is zero, v.size() - 1 is enormous. And i >= 0 is always true for unsigned i. The loop does not count down to safety. It wanders off with a clean conscience.

The clean way to write a reverse loop with size_t is shaped differently:

for (size_t i = v.size(); i-- > 0;) {
    use(v[i]);
}

Read it as: start at one-past-the-end, test the old value against zero, then decrement.

Or avoid the index with reverse iterators, which walk the vector back to front:

for (auto it = v.rbegin(); it != v.rend(); ++it) {
    use(*it);
}

The machine has a simple rule, the language has a rule, and your source code has to respect both. Unsigned wraparound is a precise modulo arithmetic system, and the bug is forgetting you entered it.

2. Negative numbers had three designs, and two’s complement won

Unsigned is easy. Every bit pattern means a nonnegative integer.

Signed integers need negative values too. Historically, machines tried several encodings.

Sign-magnitude is the obvious design:

00000101 -> +5
10000101 -> -5

The top bit is the sign. The remaining bits store the magnitude.

Problem: there are two zeros.

00000000 -> +0
10000000 -> -0

Also, arithmetic hardware becomes annoying. Adding numbers with different signs needs special cases.

One’s complement negates by flipping all bits:

00000101 -> +5
11111010 -> -5

Still two zeros:

00000000 -> +0
11111111 -> -0

Arithmetic improves, but end-around carry and signed zero still make life unpleasant.

Two’s complement negates by flipping all bits and adding one:

00000101 -> +5
11111010 -> flip bits
11111011 -> add one -> -5

For 8-bit signed integers:

00000000 -> 0
00000001 -> 1
...
01111111 -> 127
10000000 -> -128
10000001 -> -127
...
11111111 -> -1

One zero. Good.

Subtraction becomes addition of a negated value. Also good.

Example:

5 - 3 = 5 + (-3)

  00000101   5
+ 11111101  -3
-----------
1 00000010

Drop the carry beyond 8 bits:

00000010 -> 2

Two’s complement won because it makes hardware simple. Addition and subtraction can use the same adder. There is one zero. Sign extension works cleanly. Comparisons and overflow detection have regular patterns.

Sign extension is one of those details that makes two’s complement feel inevitable.

If an 8-bit two’s-complement value is -5:

11111011

and you want the same value in 16 bits, copy the sign bit to the new high bits:

11111111 11111011

That’s still -5.

For +5:

00000101

extend with zeros:

00000000 00000101

The same operation works for all values: replicate the top bit. That single rule is what makes widening a signed integer cheap in hardware.

Two’s-complement also makes -1 visually special:

8-bit:   11111111
16-bit:  11111111 11111111
32-bit:  11111111 11111111 11111111 11111111

All bits set. Masks spell that ~0u in unsigned code, and it’s what an accidental conversion of -1 to unsigned hands you: the largest unsigned value.

The C and C++ standards historically allowed other signed representations, even though mainstream hardware had already converged on two’s complement. For the machines you are likely to touch, two’s complement is the real representation. The language rules around overflow, though, are still not “signed integers wrap like unsigned.”

That distinction matters immediately.

3. Signed overflow is undefined behavior, not wraparound

This function looks like it detects overflow:

extern "C" int overflow_check(int x) {
    if (x + 1 < x) return 1;
    return 0;
}

If x is INT_MAX (the largest value an int can hold), then on a two’s-complement machine the hardware result of x + 1 would wrap to INT_MIN, which is less than x. The extern "C" just keeps the symbol name unmangled so the function is easy to spot in the assembly.

Compile at -O2, a normal optimization level:

$ clang++ -std=c++20 -O2 -S numerics_probe.cpp -o numerics_probe_O2.s

The assembly:

_overflow_check:
	mov	w0, #0
	ret

mov w0, #0 puts zero in the return register; ret returns. The check vanished. The function always returns 0.

Now compile with -fwrapv, which tells Clang to treat signed overflow as two’s-complement wraparound:

$ clang++ -std=c++20 -O2 -fwrapv -S numerics_probe.cpp -o numerics_probe_wrapv.s

Assembly:

_overflow_check:
	mov	w8, #2147483647
	cmp	w0, w8
	cset	w0, eq
	ret

cmp compares x against 2147483647, and cset w0, eq sets the return register to 1 exactly when they were equal. Now the function checks whether x == INT_MAX.

The optimizer didn’t “miss” the overflow in the first build. It used the C++ rule. Signed overflow is undefined behavior, the standard’s term for operations a valid program is assumed never to perform, which frees the compiler to reason as if they can’t happen. In any valid C++ program, x + 1 cannot overflow. Therefore x + 1 < x is impossible. Therefore return 0.

This is where integer representation and language semantics separate. The adder on this machine wraps two’s-complement bits; that’s a hardware fact. The C++ rule says signed overflow is undefined. So the optimizer may assume signed overflow never happens.

This is the source of the INT_MAX + 1 bug class. A program tries to check whether an allocation size, index, length, or offset overflowed after doing signed arithmetic. The compiler proves the check is impossible under the language rules and deletes it.

The security shape is usually:

int bytes = count * sizeof(Item);
if (bytes < count) {
    return error;
}
buffer = malloc(bytes);

or:

int end = offset + length;
if (end < offset) {
    return error;
}

If the arithmetic overflows before the check, the program has already stepped into undefined behavior. The optimizer can reason from “valid programs don’t do that” and rewrite the check away or move code around it.

The shape that survives review checks before the operation, or uses an operation that reports overflow:

int out;
if (__builtin_add_overflow(offset, length, &out)) {
    return error;
}

For multiplication:

size_t bytes;
if (__builtin_mul_overflow(count, sizeof(Item), &bytes)) {
    return error;
}

Those builtins lower to efficient overflow-checking machine code. They also tell the compiler exactly what you mean. Numerics bugs often come from writing a mathematical intention in a form the language doesn’t define.

Unsigned arithmetic wraps by definition. Signed arithmetic doesn’t.

If you need checked signed arithmetic, use a checked helper, compiler builtin, wider type, or explicit bounds check before the operation (std::numeric_limits<int>::max() below is the standard library’s spelling of INT_MAX):

if (x == std::numeric_limits<int>::max()) {
    // handle overflow
}
int y = x + 1;

Do not ask undefined behavior whether it happened. It is allowed to answer by erasing the question.

4. Integer promotions make small types less small than they look

From the probe:

uint8_t 200 * 2 has type-size 4 value 400; narrowed to uint8_t 144

The code:

uint8_t a = 200, b = 2;
auto product = a * b;
uint8_t narrowed = product;

product isn’t a uint8_t. On this machine, sizeof(product) is 4. The operands were promoted before multiplication, so the multiplication produced an int with value 400. Only when assigning back to uint8_t did it narrow modulo 256:

400 mod 256 = 144

C and C++ do integer promotions on small integer types like char, signed char, unsigned char, short, and sometimes enums. If all values fit in int, they promote to int.

So “I used uint8_t, so the arithmetic is 8-bit” is often wrong: storage is 8-bit, but the expression may be int.

Then there is the classic signed/unsigned comparison:

std::printf("%s\n", (-1 > 0u) ? "true" : "false");

Output:

-1 > 0u -> true

Because 0u is unsigned, -1 is converted to unsigned for the comparison. On a 32-bit unsigned comparison, -1 becomes 4294967295.

The expression becomes:

4294967295 > 0

True.

This shows up all the time when comparing signed indexes against sizes:

int i = -1;
if (i < vec.size()) {
    // vec.size() is size_t, usually unsigned
}

i may be converted to a huge unsigned value. Compilers warn about this if you enable the right warnings. Listen to them.

Width traps have the same flavor:

uint32_t bytes = width * height * 4;

If width and height are 32-bit and large, the multiplication can overflow before being assigned to a wider type. Writing:

uint64_t bytes = width * height * 4;

doesn’t necessarily fix it if the multiplication still happens in 32-bit first. Force a wider operand:

uint64_t bytes = uint64_t(width) * height * 4;

Numerics bugs are often type bugs wearing arithmetic clothes.

One more trap in the same family:

char c = 200;

Whether plain char is signed or unsigned is implementation-defined. On one target, that value may behave like 200. On another, it may behave like -56. If the sign matters, write signed char, unsigned char, uint8_t, or std::byte, depending on what you mean.

std::byte is especially useful when you mean “raw bits, not a tiny integer.” It refuses accidental arithmetic:

std::byte b{0xff};
// b + 1; // no

Choosing the type that matches the intent lets the compiler catch these mistakes for you.

5. Masks, shifts, and popcount are the useful bit mechanics

We need a little bit manipulation before floats, hashing, and quantization make sense.

A mask selects bits (0xff is hexadecimal for 255, the low eight bits, all set):

uint32_t low8 = x & 0xff;

A shift moves bits:

uint32_t high8 = x >> 24;
uint32_t bit10 = (x >> 10) & 1;

Set a bit:

x |= (1u << k);

Clear a bit:

x &= ~(1u << k);

Toggle a bit:

x ^= (1u << k);

Count set bits:

int n = std::popcount(x); // C++20, <bit>

Float decoding is mask and shift:

sign     = bits >> 31
exponent = (bits >> 23) & 0xff
fraction = bits & 0x7fffff

Bloom filters are mask-like thinking at a data-structure scale: hash an item to a few bit positions and set those bits. Hash tables use bit mixing and masks constantly. Quantization packs low-bit integers into machine words, then unpacks with shifts and masks inside kernels.

Floats are just the next encoding built from the same masks and shifts.

Act II: Why floating point had to be invented

6. Fixed-point works until dynamic range breaks it

Integers can’t store fractions directly. The simplest trick is fixed-point.

Pick a scale and store an integer:

stored = real_value * scale
real_value = stored / scale

If scale is 100, then:

12.34 -> 1234

Binary fixed-point does the same thing with powers of two. Q notation names how many integer and fractional bits you allocate.

Q8.8 in a signed 16-bit integer means:

8 integer-ish bits
8 fractional bits
scale = 2^8 = 256

The stored integer 3200 means:

3200 / 256 = 12.5

Addition is easy if the scales match:

1.5  -> 384
2.25 -> 576
sum  -> 960
960 / 256 = 3.75

Multiplication needs a correction:

1.5 * 2.25 = 3.375

Stored:

384 * 576 = 221184

But both inputs were already scaled by 256, so the product is scaled by:

256 * 256

To get back to Q8.8, shift right by 8:

221184 >> 8 = 864
864 / 256 = 3.375

That shift discards low bits. You have rounding questions again. Fixed-point doesn’t avoid numerical analysis; it makes the scale explicit.

This is great for domains with known range and precision: audio samples, embedded control, old game engines, financial cents if you are careful with decimal rules.

But fixed-point has a hard trade:

more fractional bits -> more precision near zero, less range
fewer fractional bits -> more range, less precision

One fixed scale can’t represent:

0.000000001
1000000000

with the same relative usefulness in a small number of bits.

Suppose you keep Q8.8 in 16 bits. The step size is:

1 / 256 = 0.00390625

The range is roughly:

-128 to +127.996

Switch to Q4.12:

step = 1 / 4096 = 0.000244140625
range roughly -8 to +7.999

More fractional precision, much less range.

Switch to Q12.4:

step = 1 / 16 = 0.0625
range roughly -2048 to +2047.9375

More range, less fractional precision.

Scientific notation solves this by moving the scale:

1.23 x 10^8
1.23 x 10^-8

Floating point is that trick, binary scientific notation packed into bits, with the scale stored per value.

7. Binary scientific notation gives us sign, exponent, and significand

In decimal scientific notation:

12345 = 1.2345 x 10^4

In binary:

13 = 1101b = 1.101b x 2^3

The useful pieces are the sign, the exponent (where the binary point moves), and the significand (the meaningful digits).

IEEE-754, the 1985 standard that defines what float means everywhere, calls the 32-bit format binary32 or fp32. It uses 32 bits as:

1 sign bit
8 exponent bits
23 fraction bits

Layout:

S EEEEEEEE FFFFFFFFFFFFFFFFFFFFFFF

For normal nonzero finite numbers, the value is:

(-1)^sign * 1.fraction * 2^(exponent - 127)

The exponent is stored with a bias of 127. That means stored exponent 127 represents actual exponent 0. Stored exponent 151 represents actual exponent 24.

The significand has an implicit leading 1 for normal numbers. The stored fraction has 23 bits, but the precision is 24 bits:

1.xxxxxxxxxxxxxxxxxxxxxxx

The exponent is stored biased instead of as a signed integer for a practical reason: unsigned exponent fields sort and compare more conveniently in hardware. Smaller stored exponent means smaller magnitude for positive normal numbers. The bias maps negative exponents into small positive stored values:

actual exponent -3 -> stored 124
actual exponent  0 -> stored 127
actual exponent 24 -> stored 151

Exponent field zero is reserved for zero and subnormals; all-ones (255) is reserved for infinities and NaNs. So normal fp32 exponents are stored from 1 through 254, representing actual exponents -126 through +127.

8. Crack a float open

Here is the decoding helper:

static uint32_t bits_of(float x) {
    uint32_t u;
    std::memcpy(&u, &x, sizeof u);
    return u;
}

Use memcpy, not pointer casting. Type-punning (reading the same bytes as a different type) through incompatible pointers can violate aliasing rules, the compiler’s assumptions about which pointer types may refer to the same memory. memcpy is the boring correct way; compilers optimize it away.

Output from the probe:

1.0f         1 0x3f800000 0 01111111 00000000000000000000000 sign=0 exponent=127 fraction=0x000000
0.1f         0.1000000015 0x3dcccccd 0 01111011 10011001100110011001101 sign=0 exponent=123 fraction=0x4ccccd
2^24         16777216 0x4b800000 0 10010111 00000000000000000000000 sign=0 exponent=151 fraction=0x000000
2^24+1       16777216 0x4b800000 0 10010111 00000000000000000000000 sign=0 exponent=151 fraction=0x000000

1.0f:

sign = 0
stored exponent = 127
actual exponent = 127 - 127 = 0
fraction = 0
value = +1.0 * 2^0 = 1.0

0.1f:

printed value = 0.1000000015
bits = 0x3dcccccd

Decimal 0.1 isn’t exactly representable in binary floating point. Just like 1/3 is 0.33333... in decimal, 1/10 repeats in binary:

0.1 decimal = 0.0001100110011001100... binary

The float stores the nearest representable value. That value is a little larger than 0.1.

That’s where the classic Python or C++ output comes from:

0.1 + 0.2 -> 0.30000000000000004

0.1 and 0.2 were never exact binary floating-point values to begin with. The addition is performed on nearby binary values, then rounded again, then printed in decimal.

Decimal fractions with denominators containing primes other than 2 repeat in binary. Since:

10 = 2 * 5

most human decimal fractions are not finite binary fractions.

Fractions like these are exact:

0.5    = 1/2
0.25   = 1/4
0.125  = 1/8
0.75   = 3/4

Fractions like these are not:

0.1 = 1/10
0.2 = 1/5
0.3 = 3/10

None of this is the float type being sloppy. The representation is binary, and some decimal values just don’t fit it.

The probe prints 0.1f as 0.1000000015, which is already a hint that the stored value isn’t 0.1, but I wanted the whole number, not a rounded preview. The bits 0x3dcccccd decode to an exact rational. The significand with its implicit leading 1 is 13421773, scaled by 2^23, and the stored exponent puts the binary point so that the exact value the float holds is:

0.1f exact fraction = 13421773/134217728
0.1f exact decimal  = 0.100000001490116119384765625

That’s not 0.1 with some display fuzz on the end. It is a specific rational number, a hair above 0.1, and the size of the gap is exact too:

0.1f - 0.1 = 1.490116119384765625e-9  (above)

The float landed about 1.49e-9 north of a tenth and stayed there. The same thing happens one size up, in double, which is where the famous printout actually comes from:

0.1(f64) exact       = 0.1000000000000000055511151231257827021181583404541015625
0.2(f64) exact       = 0.200000000000000011102230246251565404236316680908203125
(0.1+0.2)(f64) exact = 0.3000000000000000444089209850062616169452667236328125
repr(0.1+0.2)        = 0.30000000000000004

0.1 stored slightly high, 0.2 stored slightly high, their sum stored slightly high, and the ...00004 you see printed is just the shortest decimal string that round-trips back to that exact binary value. Python isn’t rounding wrong when it prints 0.30000000000000004. It’s reporting the truth about a number that was never 0.3 to begin with. (These exact expansions came from the Homebrew Python with NumPy 2.2.5, the same interpreter as the rest of the post.)

2^24:

stored exponent = 151
actual exponent = 151 - 127 = 24
fraction = 0
value = 1.0 * 2^24 = 16777216

2^24 + 1 gives the exact same bits as 2^24. The value is the same float; printing has nothing to do with it.

9. The implicit leading 1 is the whole game

Normal fp32 numbers have 24 bits of precision:

1 hidden bit + 23 stored fraction bits

That means all integers up to 2^24 are exactly representable:

0
1
2
...
16777216

At 2^24, the spacing between adjacent representable floats becomes 2, because the significand has 24 bits. Around 2^24, the representation is:

1.00000000000000000000000 x 2^24

The last significand bit is worth:

2^(24 - 23) = 2

So the next representable float after 16777216 is 16777218, not 16777217.

Another way to say the same thing:

fp32 significand precision = 24 binary digits
16777216 in binary = 1 followed by 24 zeros

To represent 16777217, you would need:

1000000000000000000000001

That’s 25 significant bits. fp32 has room for 24.

16777217 is exactly halfway between:

16777216 and 16777218

The default rounding mode is round to nearest, ties to even. 16777216 wins because its least significant significand bit is even.

The printout at the top is no longer mysterious:

16777216.0f + 1.0f -> exact real result 16777217
nearest fp32 value -> 16777216

I wanted to see the edge instead of just asserting it, so I cast a run of integers to float and back to a wider integer and checked which ones came home unchanged:

16777214 -> float 16777214.0 -> back 16777214  exact
16777215 -> float 16777215.0 -> back 16777215  exact
16777216 -> float 16777216.0 -> back 16777216  exact
16777217 -> float 16777216.0 -> back 16777216  SNAPPED
16777218 -> float 16777218.0 -> back 16777218  exact
16777219 -> float 16777220.0 -> back 16777220  SNAPPED
16777220 -> float 16777220.0 -> back 16777220  exact

Below 2^24 every integer is exact. At 2^24 the even ones survive and the odd ones snap to a neighbor. 16777217 snaps down to 16777216; 16777219 snaps up to 16777220. Half the integers are simply gone, and which half is decided by the low bit, which is the ties-to-even rule showing up as a fact about which numbers exist.

One power of two further up, the spacing doubles again to 4, and three out of every four integers vanish:

33554428 -> float 33554428.0 -> back 33554428  exact
33554429 -> float 33554428.0 -> back 33554428  SNAPPED
33554430 -> float 33554430.0 -> back 33554430  exact
33554431 -> float 33554432.0 -> back 33554432  SNAPPED
33554432 -> float 33554432.0 -> back 33554432  exact
33554433 -> float 33554432.0 -> back 33554432  SNAPPED
33554434 -> float 33554432.0 -> back 33554432  SNAPPED
33554435 -> float 33554436.0 -> back 33554436  SNAPPED
33554436 -> float 33554436.0 -> back 33554436  exact

33554430 is even but still snaps, because at 2^25 being even isn’t enough; you have to be a multiple of 4. 33554434 is even and still snaps. The survivors are ...428, ...432, ...436, the multiples of 4. The integer number line doesn’t degrade gently past 2^24; it thins out in exact powers of two, and the 1.0f that vanished at the top of the post is just the first casualty of a process that eventually deletes almost everything.

Act III: The IEEE-754 edges everyone skips

10. Subnormals and the performance cliff I didn’t get on M4

Normal floats use:

1.fraction * 2^(exponent - bias)

For fp32, the smallest positive normal is:

1.0 * 2^-126

The number line doesn’t stop there. IEEE-754 also has subnormal numbers. When the stored exponent is zero and the fraction is nonzero, the leading bit is not implicit 1. It is 0:

0.fraction * 2^-126

This creates gradual underflow. Instead of falling straight from the smallest normal number to zero, the number line continues with reduced precision down to:

2^-149

That’s std::numeric_limits<float>::denorm_min().

Subnormals exist because sudden underflow to zero is nasty. Gradual underflow preserves tiny differences and makes underflow less discontinuous.

The price is precision. Normal numbers have the hidden leading 1:

1.fraction

Subnormals have:

0.fraction

As you get closer to zero, the leading zeros eat significant bits. The smallest subnormal has only one meaningful bit. It is not a precise tiny number; it is the last rung before zero.

Still, that rung matters. Without subnormals, this can happen:

x and y are different tiny normal-ish values
x - y underflows straight to 0

With gradual underflow, the difference may survive as a subnormal. Algorithms that depend on tiny residuals get a softer landing.

People disable them because on some hardware, subnormal arithmetic is much slower. FTZ means flush-to-zero: treat subnormal results as zero. DAZ means denormals-are-zero: treat subnormal inputs as zero. These modes trade IEEE gradual underflow for speed.

I expected to measure a cliff. On this Apple M4 scalar probe, I did not.

Cleaner timing probe:

normal add             sink=1.333333e+00 time=72.594 ms
subnormal add          sink=2.802597e-45 time=46.419 ms
normal multiply        sink=3.333333e-01 time=44.946 ms
subnormal multiply     sink=0.000000e+00 time=45.674 ms

Those are single-run numbers from subnormal_probe.cpp. A re-run while editing gave 70.7 ms for the normal add and left the rest of the shape alone. The subnormal multiply is basically the same as normal multiply here. The subnormal add is faster in this particular microbenchmark, which is a sign that the benchmark is measuring pipeline details and dependency shape, not a universal law. On this M4 scalar probe, no subnormal slowdown cliff appeared.

That’s still useful. Subnormals are a special IEEE-754 region; some hardware handles them slowly; some environments flush them; you measure before assuming.

If you write numerics code that depends on gradual underflow, FTZ/DAZ settings can change results. If you write performance code that accidentally wanders into subnormals on hardware with a slow path, one tiny value can wreck a kernel (a kernel here is a hot numerical loop, no relation to the OS kernel).

“Performance cliff” is a hardware claim. IEEE-754 defines the values and operations; how fast each region runs is up to the microarchitecture. Some chips do subnormals in hardware. Some trap to microcode. Some flush under fast-math modes. Some vector paths and scalar paths differ.

The advice is boring: know whether your environment flushes, measure on the hardware you ship, and don’t assume your laptop and your accelerator agree.

11. Signed zero, infinities, and NaNs are values

The probe output:

+0.0f        0 0x00000000 0 00000000 00000000000000000000000 sign=0 exponent=0 fraction=0x000000
-0.0f        -0 0x80000000 1 00000000 00000000000000000000000 sign=1 exponent=0 fraction=0x000000
+inf         inf 0x7f800000 0 11111111 00000000000000000000000 sign=0 exponent=255 fraction=0x000000
-inf         -inf 0xff800000 1 11111111 00000000000000000000000 sign=1 exponent=255 fraction=0x000000
NaN          nan 0x7fc00000 0 11111111 10000000000000000000000 sign=0 exponent=255 fraction=0x400000
NaN == NaN -> false

Positive and negative zero compare equal:

0.0f == -0.0f // true

But they are not identical bit patterns. The sign can matter:

1.0f / 0.0f  -> +inf
1.0f / -0.0f -> -inf

Signed zero also preserves directional information near limits. Approaching zero from below and from above can matter in complex arithmetic and transcendental functions. A plain real-number classroom usually erases this detail. IEEE-754 keeps it because computers need to implement numerical functions consistently near edges.

Infinities are produced by overflow or division by zero in floating point. They are not exceptions in the C++ exception sense. They are values that propagate through arithmetic:

large * large -> +inf
1 / +inf -> +0

NaN means Not a Number. It appears from invalid operations:

0.0 / 0.0
sqrt(-1.0)
inf - inf

NaN doesn’t equal itself:

nan == nan // false

This lets you test with:

std::isnan(x)

not:

x == NAN

There are quiet NaNs and signaling NaNs. Quiet NaNs propagate without trapping (halting into an error handler) in the usual environment. Signaling NaNs are intended to raise an invalid-operation exception when used, though language and hardware behavior around them is subtle and often tamed by compilers or runtime settings.

NaNs also have payload bits. In the fp32 encoding, exponent all ones plus nonzero fraction means NaN; the fraction bits can carry diagnostic information. Many systems don’t preserve payloads carefully through all operations, but the space exists.

Quiet versus signaling NaN is encoded in the fraction field. The exact bit convention is specified by the platform’s IEEE behavior; in the common encoding, the top fraction bit distinguishes quiet NaNs. My generated NaN:

0x7fc00000

has the quiet bit set.

In everyday C++ you mostly see quiet NaNs because language and library operations generally avoid turning every invalid operation into a synchronous trap. But the floating-point environment can record exception flags:

invalid
divide-by-zero
overflow
underflow
inexact

These are not C++ exceptions. They are sticky floating-point status flags unless traps are enabled in an environment that supports them. Most application code never looks at them. Numerical libraries sometimes do.

IEEE-754 didn’t add these values for fun. They let numeric programs keep going while carrying information:

overflow -> infinity
invalid -> NaN
underflow -> subnormal or zero
signed zero -> directional limit information

That can be better than crashing or silently wrapping.

12. Rounding is a hardware protocol, not a print setting

IEEE-754 defines rounding modes. The ones most programmers meet are:

roundTiesToEven
roundTowardZero
roundTowardPositive
roundTowardNegative
roundTiesToAway

C exposes the common runtime modes through <cfenv>:

FE_TONEAREST
FE_TOWARDZERO
FE_UPWARD
FE_DOWNWARD

The default is nearest-even.

Probe around 1 + 2^-24:

nearest    1+2^-24 -> 1 bits=0x3f800000
upward     1+2^-24 -> 1.00000012 bits=0x3f800001
downward   1+2^-24 -> 1 bits=0x3f800000
towardzero 1+2^-24 -> 1 bits=0x3f800000

At 1.0, the next fp32 value is:

1 + 2^-23

Halfway is:

1 + 2^-24

Nearest-even rounds halfway to the value whose last significand bit is even. Here, 1.0 wins. Upward rounds toward +infinity, so it chooses the next float.

Ties go to even because always rounding halves upward creates a statistical bias. Ties-to-even splits exact halfway cases between lower and higher representable values depending on the low bit. It isn’t magic unbiasedness for all real data, but it avoids a simple systematic shove.

Hardware doesn’t keep infinitely many bits and then philosophize. Floating-point units compute with a few extra bits beyond the stored precision. The classic names are guard, round, and sticky:

guard:
  bit just beyond the kept precision

round:
  next bit after guard

sticky:
  OR of all remaining discarded bits

Those bits tell the unit whether the exact result is below halfway, exactly halfway, or above halfway relative to the representable candidates.

For the 1 + 2^-24 case:

1.00000000000000000000000
+0.000000000000000000000001

The extra 1 lands just beyond the last kept bit. It is the guard bit. There are no lower bits, so round and sticky are zero. That means exactly halfway.

Nearest-even asks whether the kept significand is odd or even. For 1.0, the kept low bit is zero, so it stays 1.0. If the lower candidate had an odd low bit, the tie would round upward.

Sticky is the “was there anything else?” bit. If any discarded bit beyond guard/round is 1, sticky becomes 1. This lets hardware avoid storing an unbounded tail while still distinguishing:

exactly halfway
slightly above halfway

That one sticky bit is a tiny summary of infinitely many possible discarded bits.

Rounding is why the real arithmetic result and the floating-point result are different objects. Every operation is rounded back into the finite set of representable numbers.

The 1 + 2^-24 probe shows the tie resolving down at 1.0, but the even-versus-odd half of the rule is easier to watch where the spacing is exactly 1. At 2^23 (8388608), adjacent floats are one apart, so adding 0.5f lands you exactly on the fence between two representable integers, a real halfway tie with nothing beyond it:

8388608 + 0.5 -> 8388608.0 (low bit even, ties down)
8388609 + 0.5 -> 8388610.0 (low bit odd,  ties up)

Same operation, same 0.5f, opposite outcomes. 8388608 is even, so the tie rounds back down to it and the 0.5 disappears. 8388609 is odd, so the tie rounds up, and it rounds all the way to 8388610, skipping 8388609 because that value doesn’t get to be the answer to a tie. The rounding didn’t pick the nearest integer; there wasn’t a nearest, both candidates were exactly 0.5 away. It picked the one with a zero in the last bit. If the hardware always rounded halves upward instead, a long run of ties like these would drift steadily up. Ties-to-even makes the drift cancel on average by sending half the ties each way, sorted by a bit that’s close to random across real data.

13. Machine epsilon, ULP, and the floating number line

Machine epsilon for fp32 is the distance between 1.0 and the next representable float:

epsilon = 2^-23 = 1.1920929e-7

ULP means unit in the last place. Around a given value, one ULP is the spacing between adjacent representable floats.

The spacing isn’t constant:

around 1.0:      spacing = 2^-23
around 2.0:      spacing = 2^-22
around 4.0:      spacing = 2^-21
...
around 2^24:     spacing = 2
around 2^25:     spacing = 4

Every time the exponent increases by one, spacing doubles.

That means representable floats are dense near zero and sparse far away, in relative terms. Between powers of two, there are the same number of significand steps. The absolute gap changes with scale.

This is the number-line picture:

[1, 2):       many floats spaced 2^-23 apart
[2, 4):       same count, spaced 2^-22 apart
[4, 8):       same count, spaced 2^-21 apart
...
[2^24,2^25):  same count, spaced 2 apart

So adding a tiny number to a big number can do nothing, the tiny number falls between representable slots.

ULP is also a better error unit than decimal digits when debugging floats. Saying:

the answer differs by 3 ULPs

means the result is three representable floats away from a reference at that scale. Saying:

the answer differs by 0.000001

is meaningless without knowing whether the true value is around 1, around 1e9, or around 1e-30.

Relative error and ULP error answer different questions: relative error is good for scale-free mathematical comparison, ULP error for representation-level comparison. Near zero, relative error can explode and ULP analysis gets tricky because subnormals have fixed spacing. No one metric is perfect.

14. 16777217

The opening print, one more time:

16777216.0f + 1.0f printed with %0.1f -> 16777216.0

16777216 is 2^24.

At that scale, fp32 spacing is 2.

The exact result 16777217 is halfway between adjacent floats:

16777216
16777218

Nearest-even chooses 16777216.

There’s a second way to hit the same wall:

float s = 0.0f;
for (;;) {
    float old = s;
    s += 1.0f;
    if (s == old) break;
}

Probe output:

stalled after attempt 16777217: sum stayed 16777216.0

The sum of ones stops increasing at the same value. Once the running sum reaches 2^24, adding one is exactly the doomed operation from the opening.

The mechanism is precise: fp32 has 24 bits of precision, at 2^24 the spacing is 2, adding 1 lands halfway, and nearest-even rounds back to 2^24.

That stall isn’t the end of the number, though. It’s the end of the number for that increment. If the sum freezes because 1.0f is smaller than half the local spacing, then a bigger increment should thaw it. So I let the loop stall, doubled the increment, kept going, and did that over and over to see how far a running sum of steps could climb before it ran out of exponent entirely:

increment 1 stalls at 16777216.0 (bits 0x4b800000)
increment 2 stalls at 33554432.0 (bits 0x4c000000)
increment 4 stalls at 67108864.0 (bits 0x4c800000)
increment 8 stalls at 134217728.0 (bits 0x4d000000)
increment 16 stalls at 268435456.0 (bits 0x4d800000)
...
increment 1.01412e+31 stalls at 170141183460469231731687303715884105728.0 (bits 0x7f000000)
increment 2.02824e+31 stalls at inf (bits 0x7f800000)

Every stall is an exact power of two, and every stall value is exactly double the one before it. Read the exponent field in the bit column: 0x4b800000, then 0x4c000000, then 0x4c800000. The fraction stays all zeros the whole way up; only the exponent climbs, one notch per stall. The increment that unsticks each floor is exactly half the spacing at that floor, so doubling the increment buys exactly one more power of two before it, too, becomes smaller than half an ULP.

The last finite stall is 0x7f000000, which is 2^127, the largest power of two an fp32 can hold. One more doubling of the increment and the running sum leaves the finite range: 0x7f800000 is +inf. The sum that couldn’t get past sixteen million, given a big enough step, walks the entire exponent field and falls off the top into infinity. The stall at 2^24 and the stall at infinity are the same event photographed at two scales, an addend that’s smaller than half of one ULP, and the space between them is just eight bits of exponent counting up.

Act IV: When arithmetic betrays you

15. Absorption is big plus tiny losing the tiny

The summation probe used:

1e8
then 2,000,000 copies of 1
then -1e8

The exact real sum is:

2,000,000

Local output:

oracle f64 2000000.0
naive f32  0.0
kahan f32  2000000.0
neum f32   2000000.0
pair f32   1999744.0

Why did naive fp32 give zero?

At 1e8, the fp32 spacing is 8. Adding 1.0f is too small to move the sum. Each one gets absorbed:

100000000 + 1 -> 100000000

Do that two million times, then subtract 1e8, and you get:

0

The small increments were never stored.

This is absorption, also called swamping. A small addend disappears when the running sum is too large relative to it.

The exact threshold comes from ULP spacing. Around 1e8, fp32 spacing is 8. 1.0f is less than half an ULP there, so it rounds away. If you add small values to each other first, they can become large enough to survive meeting the big value. That’s the intuition behind blocked, pairwise, and tree reductions.

Error grows with N because every addition is a rounding event. For same-sign numbers, a rough worst-case forward-error bound grows like:

O(n * epsilon)

until the bound becomes too pessimistic to be useful. The shape is enough for intuition: summing a million floats one by one is a million chances to round.

Order matters even when all values are positive:

ascending:
  small + small + small ... then big

descending:
  big + small + small ...

Ascending often loses less because the small values get a chance to accumulate before they are added to a large partial sum. Sorting before every sum is usually too expensive, but reduction algorithms exploit the same idea locally.

16. Catastrophic cancellation is losing the good bits you had

Absorption loses a small input before it matters.

Cancellation loses significant bits when subtracting nearly equal numbers.

The variance formula many people learn first is:

variance = mean(x^2) - mean(x)^2

Mathematically fine. Numerically dangerous.

Probe data:

1,000,000 values alternating between 9999.75 and 10000.25

The real variance is:

0.0625

Output:

one-pass f32 -2116246.75
two-pass f32 17266.3574
welford f32  0.0625
oracle f64   0.0625

The one-pass formula subtracts two large, close numbers:

mean(x^2) ~ 100000000
mean(x)^2 ~ 100000000

The small variance is in the low bits of the difference. But fp32 has already rounded the large terms. The subtraction exposes the rounding error. It can even produce a negative variance.

My “two-pass f32” also failed because I computed the mean with naive fp32 summation. Two-pass variance only helps if the mean is itself computed accurately enough. Welford’s online algorithm fixed this case because it updates mean and squared deviations in a more stable recurrence.

Welford’s update is:

count += 1
delta = x - mean
mean += delta / count
delta2 = x - mean
M2 += delta * delta2
variance = M2 / count

It avoids subtracting two huge nearly equal quantities at the end. It keeps the running mean and the running squared deviations together.

My first variance attempt was even uglier. I used values around 100000000.0f that differed by 1. But near 1e8, fp32 spacing is 8, so the inputs collapsed before the variance algorithm even ran. The probe’s record of that attempt:

one-pass f32 variance 2.7256422e+08
same data f64 variance 0

Converting those rounded floats to double gave a double-precision variance of 0. More precision later is not time travel.

The same pattern hits the quadratic formula:

x = (-b +/- sqrt(b^2 - 4ac)) / (2a)

When b and sqrt(b^2 - 4ac) are close, one branch subtracts nearly equal numbers. The stable version computes one root with the non-canceling sign and gets the other from:

x1 * x2 = c / a

Algebraically equivalent does not mean numerically equivalent.

17. Conditioning and stability are not the same thing

Two ideas get mixed together. Conditioning is how sensitive the mathematical problem is to input perturbations. Stability is how much extra error the algorithm introduces while solving it.

The problem can be ill-conditioned even with a good algorithm. Or the problem can be well-conditioned and your algorithm can still be terrible.

A condition number measures sensitivity. For a function f, the relative condition number is roughly:

relative change in output / relative change in input

If a tiny relative input change can cause a huge relative output change, the problem is ill-conditioned.

Example:

f(x, y) = x - y

If x and y are close, the result is tiny. Small input errors can dominate the output. Subtracting nearly equal measured quantities is ill-conditioned.

For subtraction, a useful relative sensitivity shape is:

(|x| + |y|) / |x - y|

If x - y is tiny, that denominator makes the condition number large. You can use a stable algorithm and still be limited by the fact that the input data doesn’t contain enough reliable information to determine the tiny difference.

Stability asks whether the algorithm behaves like it solved a nearby problem accurately. A stable algorithm doesn’t add much damage beyond what the conditioning forces. An unstable algorithm makes things worse.

Variance from mean(x^2) - mean(x)^2 is unstable for data with large mean and small variance. Welford’s algorithm is much more stable for the same mathematical problem.

This language matters because it prevents the wrong fix. An ill-conditioned problem calls for better data, reformulation, rescaling, or more precision. An unstable algorithm calls for a better algorithm. No algorithm removes the sensitivity inherent in the problem, but a stable one stops adding to it.

18. Compensated summation: Kahan, Neumaier, pairwise

Naive summation:

float s = 0.0f;
for (float x : xs) {
    s += x;
}

Kahan keeps a compensation term:

float sum = 0.0f;
float c = 0.0f;
for (float x : xs) {
    float y = x - c;
    float t = sum + y;
    c = (t - sum) - y;
    sum = t;
}

The compensation tries to remember the low-order part lost in the previous addition and feed it back into the next one.

For the absorption example, imagine:

sum = 100000000
x = 1

Naive fp32:

sum + x -> 100000000

Kahan sees that the addition failed to move the stored sum:

y = x - c
t = sum + y
c = (t - sum) - y

If t rounded back to sum, then (t - sum) is zero, so c becomes -1. On the next iteration, x - c is 2, feeding the lost unit back into the stream. Kahan is not mystical. It is a small memory of the rounding error.

Neumaier is a related variant that handles some cases better by compensating based on which term is larger:

float t = sum + x;
if (fabs(sum) >= fabs(x)) {
    c += (sum - t) + x;
} else {
    c += (x - t) + sum;
}
sum = t;

Pairwise summation recursively splits the array:

sum(left half) + sum(right half)

This reduces error growth from linear-ish in N to logarithmic-ish in N for many inputs, and it maps naturally to SIMD (single instruction, multiple data, one instruction operating on several values at once) and parallel reductions.

Local result on the absorption test:

naive f32  0.0
kahan f32  2000000.0
neum f32   2000000.0
pair f32   1999744.0

My simple pairwise implementation used naive blocks of 256 before combining. It improved the shape but still lost 256 ones in this adversarial order. Kahan and Neumaier recovered the exact result for this case.

NumPy’s docs say np.sum often uses a numerically better partial pairwise summation, especially when no axis is given, and may use it along the fast axis for axis-specific sums. My local NumPy result:

np.sum float32 default 1999984.0
np.sum float64 dtype 2000000.0

NumPy’s float32 sum is far better than my naive 0.0 but still not exact on this adversarial array. Asking for dtype=np.float64 (the dtype is the array’s element type) gives the oracle result here.

This “better but not exact” point matters. Pairwise summation improves error growth; it doesn’t make fp32 exact. Kahan and Neumaier improve this adversarial case; they add serial dependencies that are less convenient for wide SIMD and parallel reductions. Superaccumulators can be more exact; they cost more. There is no one summation algorithm that is best for every workload.

A library chooses a point in this space: naive is the fastest simple loop with the worst error profile, pairwise is a good parallel-friendly general reduction, compensated buys accuracy with more operations and dependencies, and a wider dtype is often the simplest safe choice when it’s available.

If the result matters, choose the accumulation strategy and dtype explicitly instead of assuming “sum” is one thing.

19. Non-associativity changes reductions

Real addition is associative:

(a + b) + c = a + (b + c)

Floating-point addition is not.

Probe:

float x = 1.0e20f, y = -1.0e20f, z = 3.14f;

Output:

(x+y)+z = 3.1400001
x+(y+z) = 0

In the first order:

x + y -> 0
0 + z -> z

In the second:

y + z -> y
x + y -> 0

z disappears when added to -1e20.

This is why parallel reductions can change answers. A serial loop has one order. A SIMD reduction has another. A tree reduction across threads has another. A distributed all-reduce across GPUs (a collective where every device ends up holding the combined result) has another, and atomics (updates the hardware applies one at a time, in whatever order they arrive) can make the order nondeterministic.

The difference is usually small. Sometimes it isn’t. If your training run differs in the fourth decimal across devices, reduction order is one suspect.

The annoying part is that performance pushes toward changed order:

serial loop:
  (((a0 + a1) + a2) + a3) ...

SIMD:
  lane-local partial sums, then horizontal reduction

threads:
  per-thread chunks, then combine

GPU:
  warp reductions, block reductions, atomics or trees

distributed:
  ring or tree all-reduce across devices

(Warps and blocks are the GPU’s thread-group sizes.) Every shape has a different rounding history. Deterministic reductions are possible, but you have to ask for them and pay for them.

Act V: Doing it fast and wide

20. Dot product is slow and wrong in the naive form

A dot product:

float dot(const float* a, const float* b, size_t n) {
    float s = 0.0f;
    for (size_t i = 0; i < n; ++i) {
        s += a[i] * b[i];
    }
    return s;
}

It’s slow and it’s wrong in a specific way: one scalar multiply-add chain doesn’t use the full machine, and one fp32 accumulator rounds every addition in sequence.

Modern CPUs want wide work. On Apple Silicon, NEON/ASIMD (Arm’s SIMD instruction set) vector registers can process multiple floats per instruction. A good dot product unrolls, vectorizes, uses multiple accumulators to break dependency chains, and reduces partial sums at the end.

But wider and faster doesn’t automatically mean more accurate. Multiple accumulators change the reduction order. FMA changes rounding. Accumulating in fp32 versus fp64 changes the result. BLAS libraries choose kernels with specific tradeoffs.

In the dot product, the speed choices and the accuracy choices are the same choices.

A scalar dependency chain is especially limiting:

s0 = 0
s1 = round(s0 + a0*b0)
s2 = round(s1 + a1*b1)
s3 = round(s2 + a2*b2)
...

Each step depends on the previous one. The CPU cannot freely overlap all additions because the next addition needs the previous sum.

A faster kernel uses multiple accumulators:

s0 accumulates elements 0,4,8,...
s1 accumulates elements 1,5,9,...
s2 accumulates elements 2,6,10,...
s3 accumulates elements 3,7,11,...

Then it reduces s0+s1+s2+s3 at the end. This breaks the dependency chain and exposes more instruction-level parallelism (independent instructions the CPU can overlap). It also changes the rounding order.

SIMD does the same idea in registers. A vector register might hold four fp32 lanes:

[a0, a1, a2, a3] * [b0, b1, b2, b3]

producing four products at once. But those four lanes still need to become one scalar dot product eventually, so a horizontal reduction happens. That reduction order is part of the numerical result.

21. NEON and FMA: one rounding instead of two

FMA means fused multiply-add:

a * b + c

computed as one fused operation with one final rounding, not:

round(a*b)
round(product + c)

That’s both faster and often more accurate.

My arm64 assembly for the dot probe at -O3 contains:

fmadd	s0, s1, s2, s0

Even the plain loop:

s += a[i] * b[i];

uses fmadd in the scalar remainder path under this compiler configuration. The explicit builtin version:

s = __builtin_fmaf(a[i], b[i], s);

also emits:

fmadd	s0, s1, s2, s0

I expected to show a contrast. The compiler did not cooperate, and that is the actual finding. On this target, with these flags, Clang is happy to fuse the multiply-add.

If you need strict control over contraction, the knobs are -ffp-contract, the fast-math flags, pragmas, and explicit calls like std::fma. Numerical reproducibility often starts by pinning these choices. Otherwise two builds can differ because one fused and the other rounded twice.

FMA can improve accuracy, but it can also break tests that assumed two roundings. Example shape:

float r1 = (a * b) + c;       // may fuse
float r2 = std::fma(a, b, c); // requests fused behavior

If the compiler contracts the first expression, r1 and r2 match. If contraction is disabled, r1 may round after multiplication and then again after addition. The difference is usually one or a few ULPs. For a chaotic simulation or a long reduction, that can branch into visibly different output. “More accurate” and “bitwise same as before” are different requirements.

That -ffp-contract knob stopped being abstract when I tried to catch a single rounding in the act. There’s a classic move called an error-free transform: for a product p = x*y, the value fma(x, y, -p) is exactly the rounding error the multiply threw away, because the fused op computes x*y to full width and only rounds once, after subtracting p. If instead you compute (x*y) - p with two separate rounded operations, the x*y rounds to p first and the subtraction gives a clean zero. So the two spellings should disagree: one hands you the lost bits, the other hands you nothing.

I picked the literals x = 1.0000001f and y = 1.0000002f (which store as 1.00000012 and 1.00000024), whose exact product needs more than 24 bits, and printed both spellings:

p = x*y (rounded)      = 1.00000036 bits 0x3f800003
fma(x,y,-p) = rounding = 2.84217094e-14 bits 0x29000000
(x*y)-p                = 2.84217094e-14 bits 0x29000000

Both spellings gave the same nonzero error. The (x*y) - p line was supposed to be zero, and it wasn’t. The compiler had contracted (x*y) - p into a single fma too, so my “unfused” control was fused behind my back. This is the same wall the earlier fmadd probe hit, and it’s the default: at -O2, Clang runs -ffp-contract=fast, which lets it fuse a multiply and an add inside a single expression whenever it likes. To actually get two roundings I had to say so:

$ clang++ -std=c++20 -O2 -ffp-contract=off ladder_probe.cpp -o ladder_off
fma(x,y,-p) = rounding = 2.84217094e-14 bits 0x29000000
(x*y)-p                = 0 bits 0x00000000

Now the control is 0 and the fused version keeps the lost bits. And those lost bits are real information, not noise. Adding them back to p in double precision reconstructs the exact product:

double exact product   = 1.0000003576278971
p + fma_error (double) = 1.0000003576278971

The multiply threw away 2.84e-14, fma caught it, and the two pieces together are the answer the float couldn’t hold. This is the whole reason compensated algorithms and double-float tricks lean on FMA: it is the cheapest way to see your own rounding error. It’s also why the same flag that makes your code faster can quietly change which of these two programs you actually compiled.

22. What BLAS buys besides speed

BLAS is the Basic Linear Algebra Subprograms interface: vector and matrix kernels like dot, axpy (y = a*x + y), gemv (matrix times vector), and gemm (matrix times matrix). Libraries such as OpenBLAS, BLIS, MKL, Accelerate, and vendor GPU libraries implement those kernels with tuned assembly or intrinsics.

This NumPy build reports:

blas:
  name: openblas
  version: 0.3.29

SIMD:
  baseline: NEON, NEON_FP16, NEON_VFPV4, ASIMD
  found: ASIMDHP

So on this machine I should not claim NumPy dispatches to Accelerate (Apple’s built-in BLAS). It does not, at least for this Homebrew Python. It uses OpenBLAS.

BLAS buys speed: blocked memory access, cache reuse, SIMD kernels, unrolling, prefetching (fetching data before the loop asks for it), threading, and microarchitecture-specific scheduling. Sometimes it also buys accuracy: multiple accumulators, pairwise-ish reduction shapes, FMA, and wider internal accumulation for some low-precision kernels.

But BLAS doesn’t mean “mathematically exact.” It means “a highly engineered implementation of a specified operation.” The accumulation dtype and reproducibility guarantees depend on the routine, library, build flags, and hardware. Serious numerical code grows knobs for all of this, accumulation dtype, deterministic reductions, thread count, fast-math, BLAS backend, because the backend is part of the numerical environment.

Matrix multiply shows the memory reason. A naive triple loop:

for (int i = 0; i < M; ++i)
  for (int j = 0; j < N; ++j)
    for (int k = 0; k < K; ++k)
      C[i][j] += A[i][k] * B[k][j];

reloads data wastefully. A blocked GEMM works on tiles (small sub-matrices sized to fit in cache):

load a tile of A
load a tile of B
update a tile of C many times while data is hot

The numerical effect is that each output element is a dot product, but the implementation uses a carefully chosen accumulation order and register blocking. Vendor libraries tune tile sizes to cache, registers, SIMD width, and core count.

“NumPy is fast” usually means NumPy handed the hot loop to a BLAS library whose authors spent years on the memory hierarchy.

23. Low precision inputs often accumulate wider

The ML world relies on a crucial trick:

store/multiply in fewer bits
accumulate in more bits

I wanted to show the exact local proof:

bf16 sum of ones stalls at 256
bf16 matmul accumulates to 4096

This NumPy build has no np.bfloat16, and PyTorch/JAX are not installed. So I wrote a tiny software bfloat16 emulator with round-to-nearest-even.

Output:

255 bf16 serial sum 255.0
256 bf16 serial sum 256.0
257 bf16 serial sum 256.0
4096 bf16 serial sum 256.0
4096 bf16 inputs, fp32 accumulate, bf16 output 4096.0
bf16 bits for 256 0x4380
bf16 bits for 4096 0x4580

Why does serial bf16 sum stall at 256?

bfloat16 has:

1 sign bit
8 exponent bits
7 fraction bits

That’s 8 bits of precision including the hidden leading 1. At 256 (2^8), the spacing is:

2^(8 - 7) = 2

Adding 1 lands halfway between 256 and 258. Nearest-even chooses 256.

But if a matmul kernel reads bf16 inputs, multiplies them, and accumulates in fp32, then:

4096 copies of 1*1

can accumulate exactly in fp32 to 4096 before rounding the output back to bf16. 4096 is exactly representable in bf16, so the output is 4096.

This is the central low-precision matmul idea: low-precision storage saves bandwidth, wider accumulation saves the reduction. Without wider accumulation, low-precision training would be far more brittle.

The same principle appears in integer quantized inference. Int8 multiply may accumulate into int32:

int8 * int8 -> int16-ish product
sum many products -> int32 accumulator
requantize output

If you accumulated hundreds or thousands of int8 products back into int8, overflow would be immediate. The low-bit inputs only work because the accumulator is wider.

Act VI: Fewer bits on purpose

24. ML starves for bits because memory and bandwidth dominate

Some grounding first, for a reader who has never trained a model. A neural network’s weights are the numbers it learns; training nudges them, step by step, until the model’s outputs get less wrong. Each training step runs the model forward, measures how wrong the output was (the loss), computes a gradient for every weight, which direction to nudge it, and how hard, and lets an optimizer apply the nudges. The optimizer usually keeps running statistics per weight (optimizer state), and the intermediate values each layer produced on the way forward (activations) have to stick around for the backward pass.

A billion parameters at different precisions:

fp32: 4 GB
fp16: 2 GB
bf16: 2 GB
int8: 1 GB
int4: 0.5 GB

That’s just weights. Training adds gradients, activations, optimizer state, temporary buffers, and communication. Adam, a common optimizer, keeps two moment vectors, usually fp32. A model with 1 GB of weights can easily need many GB to train.

Bandwidth matters as much as capacity. If a layer reads 1 GB of weights per token batch, int8 turns that into roughly 250 MB and int4 into roughly 125 MB, as long as the kernel keeps the weights compressed. If a kernel is memory-bound, halving bytes can nearly double throughput; if a distributed job all-reduces gradients, halving gradient bytes reduces network pressure. This is the roofline argument, performance is capped by the smaller of compute peak and memory bandwidth times arithmetic intensity, and low precision attacks the bandwidth side while often raising the hardware compute peak too.

The cost is numerical: less range, less precision, more rounding, more saturation, more need for scaling. And the scale metadata is real. One scale per tensor is tiny metadata and a worse fit; one scale per channel fits better and costs more; one scale per small group of values is the common compromise in large-language-model (LLM) quantization. Low precision is a data layout as much as a dtype.

25. The format zoo as consequences

Each format spends its bits differently.

fp16, IEEE binary16, spends its 16 bits as 1 sign, 5 exponent, 10 fraction, with a max finite value of 65504. Local NumPy:

np.float16(65504) * np.float16(2) -> inf
np.float16(1e-8) -> 0.0
np.float16(1e-5) -> 1e-05

fp16 has more mantissa bits than bf16 but far less exponent range. It can underflow and overflow easily in training.

bfloat16 spends the same 16 bits as 1 sign, 8 exponent, 7 fraction: the exponent width of fp32 with a much shorter fraction. It keeps range and spends precision, which suits deep learning because gradients and activations can have wide dynamic range.

tf32 (1 sign, 8 exponent, 10 fraction) is what NVIDIA introduced for tensor cores, the GPU’s dedicated matrix-multiply units: fp32-like range, reduced mantissa, fast matmul. It’s less a storage type than a compute mode for multiplying rounded inputs with wider accumulation.

The fp8 formats split one argument two ways: E4M3 spends 4 exponent bits and 3 mantissa bits (more precision, less range); E5M2 spends 5 and 2 (more range, less precision).

int8 and int4 aren’t floating point. They need scales, real_value ~= scale * integer_value, and sometimes zero-points for ranges that aren’t centered on zero: real_value ~= scale * (integer_value - zero_point).

The pattern across the zoo: range can live in exponent bits or in shared scale metadata, and precision can live in mantissa bits or be partly recovered by wider accumulation. Which format wins depends on the workload and on the failure mode you can tolerate.

Every format is a budget.

26. Softmax overflow and the max-subtraction trick

Naive softmax:

softmax(x_i) = exp(x_i) / sum_j exp(x_j)

Local NumPy, on three logits (the raw scores a model produces before they become probabilities):

logits = np.array([1000, 1001, 1002], dtype=np.float64)
np.exp(logits) / np.exp(logits).sum()

Output:

[nan nan nan]

exp(1002) overflows.

The fix subtracts m = max_j x_j from every logit before exponentiating:

softmax(x_i) = exp(x_i - m) / sum_j exp(x_j - m)

This is exact algebra, multiply top and bottom by exp(-m), and now the largest exponent is exp(0) = 1, with the rest below it.

Local output:

[0.09003057 0.24472847 0.66524096]

Logsumexp uses the same trick: log(sum exp(x_i)) = m + log(sum exp(x_i - m)).

Layernorm (a layer that recenters and rescales activations) has a different numerics trap: variance. If you compute variance as mean(x^2) - mean(x)^2, you’re back in cancellation land. Stable implementations use better reduction patterns and often accumulate in fp32 even when inputs are lower precision.

Attention, the transformer operation that scores every token against every other token, has both problems at once: big dot products (accumulation precision), a softmax over those scores (overflow without max subtraction), and another matmul-like accumulation after it. Transformer kernels are numerics code as much as GPU plumbing. FlashAttention is famous for memory movement, but its online softmax recurrence also preserves the max-subtraction stability while tiling the computation.

27. Mixed-precision training leaks precision on purpose

A real mixed-precision training step keeps activations and weights in fp16 or bf16 for the forward and backward passes, accumulates matmuls and convolutions in fp32, scales the loss to avoid gradient underflow, keeps fp32 master weights and fp32 optimizer state, and casts updated weights back to low precision after each step.

Loss scaling exists because fp16 gradients can underflow to zero. Multiply the loss by a scale S before backprop (backpropagation, the backward pass that computes all the gradients), and every gradient comes out scaled by S; divide them back before the optimizer step. If gradients overflow, reduce S. Dynamic loss scaling automates this.

Master weights stay in fp32 because a weight around 1.0 in fp16 has spacing of about 2^-10. A tiny optimizer update may be smaller than one fp16 ULP and vanish if applied directly. In fp32, small updates get to accumulate over time.

Optimizer state stays in fp32 for the same reason. Adam’s moments are long-running sums, exactly the kind of recurrence that suffers when precision is too low:

m_t = beta1 * m_{t-1} + (1-beta1) * g_t
v_t = beta2 * v_{t-1} + (1-beta2) * g_t^2

Stochastic rounding is another tool: round probabilistically based on distance to the neighboring representable values, so a tiny update below one ULP still moves the stored value sometimes, preserving it in expectation.

Mixed precision is a leak-management system. The leaks: activations can overflow or underflow, fp16 gradients can underflow, weight updates can fall below one low-precision ULP, optimizer recurrences accumulate rounding error, and reductions are non-associative and usually parallel. The patches, loss scaling, fp32 accumulators, fp32 master weights, fp32 optimizer state, stochastic rounding, careful normalization kernels, exist one per leak.

28. Quantization is scale management

Symmetric int8 quantization maps real ~= scale * q with q in [-127, 127], choosing scale = max(abs(x)) / 127 and q = round(x / scale). Asymmetric quantization adds a zero-point, real ~= scale * (q - zero_point), which helps when the range isn’t centered on zero, especially for activations.

Small worked example:

values = [-1.0, -0.5, 0.0, 0.5, 1.0]
int8 symmetric range = [-127, 127]
scale = 1.0 / 127
q = round(value / scale)

-1.0 -> -127
-0.5 -> -64
 0.0 -> 0
 0.5 -> 64
 1.0 -> 127

Dequantizing gives real ~= q * scale, and 64 / 127 is not exactly 0.5. Quantization error is now part of the model.

If one outlier is present:

[-1, -0.5, 0, 0.5, 100]

the scale becomes 100/127, and the small values collapse into very few integer levels. This is why outliers matter so much in LLM quantization.

Per-tensor quantization uses one scale for the whole tensor. Per-channel uses one scale per output channel or group, and usually preserves accuracy better because different channels can have different ranges.

GPTQ and AWQ are the post-training quantization methods you’ll meet first in the LLM world. GPTQ quantizes weights while accounting for approximate second-order sensitivity, layer by layer; AWQ picks scales that protect the weight channels activations care about most. I didn’t measure either one here, they’re names to chase, not results I’m reporting.

Inference kernels dequantize inside the kernel: load packed int4 weights, unpack nibbles, apply scales while feeding SIMD or tensor-core operations, accumulate wider, and maybe requantize the output. The dequantized values exist briefly in registers, never as an expanded matrix in RAM; materializing it would throw away the bandwidth win.

Quantization is numerics plus memory layout plus kernel engineering.

Act VII: The frontier

29. Determinism across hardware is not free

Two devices can disagree in the fourth decimal because they didn’t add in the same order.

Parallel reductions use trees. GPU reductions use warps, blocks, atomics, and sometimes nondeterministic scheduling. Distributed reductions depend on topology and algorithm. Atomics serialize in some order, but not necessarily the same order across runs. Since floating addition is non-associative, a different order gives a different rounded result. FMA availability changes results too, fused is one rounding, separate multiply and add is two. And hardware differs in internal precision, math-library approximations, flush-to-zero behavior, and BLAS kernels.

Deterministic ML modes exist, but they often cost speed. They constrain algorithm choices so reductions happen in reproducible orders and avoid nondeterministic kernels.

The useful question is which contract you asked for. Fastest lets the hardware and library choose aggressive orders and kernels. Reproducible constrains order, precision, and kernel choice, and pays for it in speed. Correctly rounded is a much stronger contract and often much more expensive.

30. Block floating point, MX formats, and posits

Block floating point shares one exponent across a block of numbers, leaving many low-bit mantissas per exponent. That works well when values in a block have similar scale, saves exponent bits, and can map nicely to hardware. It struggles when one block contains both huge and tiny values.

MX formats, the recent microscaling formats, push this idea into ML hardware design: small element formats plus shared scale metadata. Same pressure as everywhere else in this act: fewer bits, enough dynamic range, tensor cores kept fed. Posits are a different proposed number system with tapered precision, encoded with sign, regime, exponent, and fraction fields; they give more precision near 1 and different tradeoffs than IEEE floats.

I measured none of these, no hardware on my desk speaks MX or posits. They’re the reading list for this post, not the results.

The historical lesson from IEEE-754 still applies, though. A number format is an encoding plus rounding rules plus exception behavior plus hardware cost plus a software ecosystem, and the ecosystem part is how formats can be technically elegant and still lose.

31. The un-run experiment

The experiment I haven’t run yet:

write a Metal kernel
choose an input format
choose an accumulation format
sum ones until the accumulator stalls
repeat for fp32, fp16, bf16-emulated or hardware-supported paths
repeat for vector reductions and tree reductions
record the stall point and throughput

The CPU numbers say where each format should stall. The GPU numbers are the open thread.

Act VIII: the GPU numbers, months later

32. Update: PyTorch showed up, so I ran the un-run experiment

I wrote section 31 with no PyTorch on the machine. Months later, doing unrelated work, I found a second Python on this same M4, an Anaconda install carrying torch 2.9.1, and torch.backends.mps.is_available() returned True. MPS is PyTorch’s Metal backend; it runs tensor ops on the M4’s GPU through Metal, which is the exact door section 31 said was closed. This isn’t the hand-written Metal shader I sketched there, and it isn’t even the same interpreter as the rest of this post (this one has numpy 1.26.4, not the 2.2.5 in the header), so read it as a second environment sharing one piece of silicon. But it can sum ones on the GPU in a chosen precision until the accumulator stalls, which was the whole ask.

The most direct version first: accumulate one at a time, on the GPU, in fp16 and bf16, and stop when the sum stops moving.

fp16 serial GPU sum stalled at attempt 2049, sum stuck at 2048.0
bf16 serial GPU sum stalled at attempt 257, sum stuck at 256.0

bf16 stalls at 256 on the actual Metal GPU. That’s the number my software bf16 emulator predicted back in section 23, now confirmed by real hardware doing real bf16 adds one at a time. fp16 stalls at 2048, which is 2^11: fp16 keeps 10 stored fraction bits plus the hidden 1, so 11 bits of precision, so integers stay exact up to 2^11, and 2^11 + 1 is the doomed 16777216 + 1 shape one format down. Every format has its own sixteen-million, and for fp16 it is two thousand and forty-eight.

33. The reduction order rewrote the stall

Then I stopped hand-rolling the loop and asked the library to sum, which quietly changed the answer.

fp16 N=   256  sum(same dtype)=     256.0  sum(fp32 acc)=     256.0
fp16 N=  4096  sum(same dtype)=    4096.0  sum(fp32 acc)=    4096.0
fp16 N= 65536  sum(same dtype)=       inf  sum(fp32 acc)=   65536.0
bf16 N=   256  sum(same dtype)=     256.0  sum(fp32 acc)=     256.0
bf16 N=  4096  sum(same dtype)=    4096.0  sum(fp32 acc)=    4096.0
bf16 N= 65536  sum(same dtype)=   65536.0  sum(fp32 acc)=   65536.0

The serial bf16 loop stalled at 256. But torch.sum of 4096 bf16 ones, still in bf16, returns 4096 exactly. Same dtype, same data, and it went straight through the stall that trapped the serial version. The reason is the reduction shape. A tree reduction adds ones in pairs: 1+1 = 2, then 2+2 = 4, then 4+4 = 8, and every partial sum is a power of two, which bf16 represents exactly. The serial loop hits 256 + 1 and dies; the tree loop only ever computes 256 + 256 and sails on. Non-associativity from section 19 usually shows up as a small disagreement in a low bit. Here it is the difference between stalling at 256 and reaching 4096 in the same eight-bit-precision format. Order isn’t a rounding detail sitting on top of the answer; sometimes it is the answer.

The 65536 row has its own surprise. In fp16, torch.sum of 65536 ones returns inf. Not a stall, an overflow: fp16’s largest finite value is 65504, and a tree partial sum climbs past it (32768 + 32768 = 65536) and saturates to infinity. bf16 in the same spot returns 65536 cleanly, because bf16 spends its bits on fp32’s exponent range and can’t overflow until around 3e38. This is the range-versus-precision split from the format zoo turned into a live result: fp16 has more fraction bits and here dies of overflow, bf16 has fewer and survives. And asking for a wider accumulator (sum(dtype=torch.float32)) fixes the fp16 case, 65536 instead of inf, which is the low-precision-in, wide-accumulate-out pattern that makes mixed precision work at all.

Matmul tells the same story from the compute side:

fp16 K= 16384  (ones@ones)[0,0]=   16384.0  out dtype=torch.float16
bf16 K= 16384  (ones@ones)[0,0]=   16384.0  out dtype=torch.bfloat16

A 1xK times Kx1 matmul of ones is a dot product of K ones, and both low-precision matmuls return 16384 exactly, far past their serial stall points of 2048 and 256. The GPU’s matmul is not a naive same-precision accumulate; it is reducing wider, or in a tree, or both. I tried to pin down which by pushing K past fp16’s overflow ceiling, and that test came back inconclusive: at K = 65536 the fp16 matmul returns inf, and I can’t tell whether the accumulator overflowed or the fp16 output simply can’t hold 65536 (it can’t; 65536 rounds to inf in fp16 no matter how it was computed). So the exact internal accumulation dtype of the MPS matmul is something I did not manage to measure. I know it’s wider than a serial fp16 add, and past the output’s range I can’t see inside it.

34. I expected the two devices to disagree, and mostly they did

Section 29 asserts that two devices can disagree in the fourth decimal because they add in a different order. I’d never actually watched it happen on this machine, so I summed the same fp32 random vector on the CPU and on the GPU and compared bits:

N=2^12  cpu=    35.893959045410156  gpu=     35.89396667480469  DIFFER  (f64 ref 35.893965)
N=2^18  cpu=      397.376708984375  gpu=    397.37677001953125  DIFFER  (f64 ref 397.376754)
N=2^20  cpu=    -708.4118041992188  gpu=    -708.4116821289062  DIFFER  (f64 ref -708.411797)
N=2^22  cpu=     586.9429931640625  gpu=     586.9424438476562  DIFFER  (f64 ref 586.942748)
N=2^24  cpu=        -955.244140625  gpu=        -955.244140625  IDENTICAL  (f64 ref -955.244187)

There it is, measured. At N = 2^22 the CPU says 586.9430 and the GPU says 586.9424, and they split in the fourth significant figure, with the fp64 reference 586.942748 sitting between them. Neither device is wrong; they rounded a different sequence of partial sums. The first time I ran this, with a different random seed and only N = 2^22, the two happened to match to the bit, and I nearly wrote down that the disagreement was a myth. It isn’t; it’s data-dependent. The 2^24 row here matched too, and I don’t have a clean story for that one, maybe the two reduction trees line up at that size, maybe it’s luck. I’m not going to pretend I know.

The stranger result is the matmul:

matmul 512x512 fp32: max abs diff = 0.000e+00
bitwise identical = True

A 512-by-512 fp32 matmul is bit-for-bit identical between the CPU and the Metal GPU, even though the elementwise sum wasn’t. Two different pieces of silicon, two completely different kernels, and every one of the 262144 output elements agrees to the last bit. I don’t have a confident explanation. My guess is that both backends land on the same blocked reduction order for this shape, or round each fma the same way, but I didn’t disassemble either kernel, so that’s a hypothesis, not a finding. The disagreement is real for the reductions and absent for this matmul, and the honest reading is that “different hardware gives different bits” is true often enough to ruin a reproducibility test and false often enough that you can’t lean on it as a diagnostic.

The Metal shader I actually described in section 31, hand-written, one thread per accumulator, sweeping input and accumulation formats, is still unwritten. torch’s MPS backend answered the format-stall question and the cross-device question, but it hid the accumulator dtype behind its own kernels. The thing I still can’t see is what precision the GPU is really adding in when I don’t get to choose it, and the only way to find out is to stop borrowing torch’s kernels and write the reduction myself.

Sources and artifacts

Local scratch directory:

/private/tmp/numerics-notes

A copy of every probe lives in this blog’s repo under writing-notes/probes/numerics/, so the runs survive the next reboot.

Local tools and versions:

Apple M4
Apple clang version 15.0.0 (clang-1500.3.9.4)
Target: arm64-apple-darwin24.6.0
NumPy 2.2.5
OpenBLAS 0.3.29

The Act VIII GPU numbers come from a second interpreter on the same M4, an Anaconda Python with torch 2.9.1 and numpy 1.26.4, using the Metal (MPS) backend. It’s a different environment, flagged as such wherever its numbers appear.

Local artifacts used:

  • numerics_probe.cpp: integer wraps/promotions, float bit decoding, fp32 stall, special values, rounding modes, first subnormal timing, summation, non-associativity, cancellation
  • subnormal_probe.cpp: cleaner normal/subnormal timing on Apple M4
  • variance_probe.cpp: one-pass variance failure, naive two-pass failure, Welford fix
  • fma_probe.cpp: arm64 fmadd assembly for dot-product loops
  • ladder_probe.cpp: the exponent ladder from 2^24 to +inf, which integers survive past 2^24 and 2^25, ties-to-even at 2^23, and the FMA error-free transform with the -ffp-contract=off contrast
  • ml_probe.py: NumPy/OpenBLAS config, bfloat16 availability check, fp16 sum behavior, NumPy sum behavior. The fp16 overflow/underflow and softmax runs were separate one-liners against the same interpreter, not part of this file.
  • bf16_emulation.py: software bfloat16 rounding and accumulation-width proof
  • exact_probe.py: exact rational value of 0.1f and the exact decimals behind 0.1 + 0.2 (Homebrew Python, NumPy 2.2.5)
  • gpu_probe.py, gpu_probe2.py: serial and tree reductions, matmul accumulation, and CPU-versus-GPU reduction disagreement on the M4 GPU via PyTorch MPS (Anaconda Python, torch 2.9.1)

References I leaned on: