11 - Floating-Point Numbers: Why 0.1 + 0.2 Is Not 0.3

Python for data, part 4: decimal numbers in a computer (IEEE 754) and why 0.1 + 0.2 != 0.3

Machine Learning Engineer course · Module 01 · Python for data · Lesson 4

You write 0.1 + 0.2 == 0.3 and Python answers False. It is not a defect of the language and you do not fix it by changing library: it is the consequence of how every computer writes numbers with a decimal point, and in data work it shows up in forms far less innocent than a failing comparison. This lesson follows lists, tuples, dictionaries and sets and functions, files, exceptions and objects: from the bits you work out why 0.1 does not exist in binary, you compute the distance between two neighbouring numbers and you meet three cases of the trade in which floating-point arithmetic breaks an algorithm that is correct on paper, namely the variance with the "textbook" formula, the sum of many terms and the softmax.

The commonplace says floats are approximate but the error is negligible. In the variance computed with the convenient formula the error produces a negative result, a quantity that by definition cannot be: you will see it run further down, with the real output.

What it is, why it exists, when it is NOT enough

The problem

A computer has finite memory, the real numbers are infinite. You must choose a way of representing a subset of the reals in a fixed number of bits. The IEEE 754 standard (1985, revised in 2008 and 2019) is the universal answer: almost every processor implements it, so the basic operations (sum, product, division) give the same result on different machines.

The analogy: scientific notation with limited space

Imagine you can write numbers only with 3 significant digits and an exponent: 1.23 × 10⁴. You represent 12300 well, but not 12345 (it becomes 1.23 × 10⁴ = 12300) and you cannot tell 12300 from 12349. The bigger the number, the bigger the gap between two representable numbers. A float is this, in binary, with 53 significant digits.

When floats are NOT enough

• Money: use integers (cents) or decimal.Decimal.
• Exact equality between results of computations: use tolerances.
• Very large integers (beyond 2⁵³): a float loses them; use int, which has arbitrary precision.
• Exact symbolic computation (fractions): fractions.Fraction.

The common thread is that a float is a good choice for measured quantities, where the uncertainty of the measurement is already bigger than the rounding error, and a bad choice for quantities that are counted or compared exactly. For money the problem is twofold: cents are not representable in binary and the error grows with the amount, as exercise 2 below shows.

Setup

You only need NumPy; the rest is the standard library (struct, math, decimal, fractions). The examples were each run twice, with the same output, on Python 3.14.4 and NumPy 2.3.5. The install command below is a template to adapt to your environment (not run here, where NumPy was already present).

python -m pip install numpy

First example: the famous 0.1 + 0.2

import math

print(0.1 + 0.2)
print(0.1 + 0.2 == 0.3)
print(math.isclose(0.1 + 0.2, 0.3))

Output: 0.30000000000000004, False, True.

Python is not wrong: it prints the shortest string that identifies the float obtained from the computation, and that float is not 0.3. The third comparison, with math.isclose, is the remedy we develop later. To see why the result is exactly that, let us look at the bits.

Theory: the float64 format, derived

Structure (IEEE 754 binary64)

A Python float takes 64 bits, divided into three fields:

For "normal" numbers the value is

x = (−1)^s × (1 + f/2⁵²) × 2^(e − 1023)

The "1 +" in front of the fraction is the hidden bit: in normalised form the mantissa always has the shape 1.xxxx… in binary, so that leading 1 is not stored. This gives a free digit: 52 stored bits, 53 of precision.

Let us look at the bits of 0.1

import struct

def campi_float64(x):
    """Return (sign, exponent, fraction) as bit strings."""
    n = struct.unpack(">Q", struct.pack(">d", x))[0]
    b = f"{n:064b}"
    return b[0], b[1:12], b[12:]

s, e, f = campi_float64(0.1)
print(s, e, f)
print("exponent:", int(e, 2), "-> 2^", int(e, 2) - 1023)
mantissa = 1 + int(f, 2) / 2**52
print(mantissa * 2 ** (int(e, 2) - 1023))

Output: sign 0, exponent 01111111011, fraction 1001100110011001100110011001100110011001100110011010; then exponent: 1019 -> 2^ -4; then 0.1.

How to read the code. struct.pack(">d", x) converts the float into 8 bytes according to the IEEE 754 standard (the sign > fixes the byte order, d means "double", that is 64 bits); struct.unpack(">Q", ...) reinterprets the same 8 bytes as an unsigned 64-bit integer. The string f"{n:064b}" writes that integer in binary with 64 digits. Then the program cuts out the three fields of the figure: the first bit is the sign, the next 11 are the exponent, the last 52 the fraction. The last line rebuilds the value with the format's formula (mantissa × 2^(e − 1023)), and gives 0.1: the formula really is what the computer does.

Read: e = 1019, so the factor is 2^(1019−1023) = 2⁻⁴ = 1/16. The mantissa is about 1.6, and 1.6/16 = 0.1. But look at the fraction: 1001 1001 1001 … repeats. Why?

Derivation: 0.1 in binary never ends

To convert a decimal fraction to binary you repeatedly multiply by 2 and take the integer part as the next digit.

At step 6 the remainder is 0.2 again, already seen at step 2: the cycle 0011 repeats forever. Therefore

0.1₁₀ = 0.0001 1001 1001 1001 …₂

(the bits 1001 1001 … you saw in the fraction are this same cycle, read from the first 1 after shifting the point by four places for the exponent −4; the last group is 1010 instead of 1001 because the cut rounds up). A number has a finite binary expansion if and only if, written as an irreducible fraction p/q, the denominator q is a power of 2. 0.1 = 1/10 = 1/(2·5): it contains a factor 5, so the expansion is infinite. Exactly like 1/3 = 0.3333… in decimal, because 3 does not divide 10.

The computer must cut after 53 significant digits and round. The number it gets is slightly different from 0.1. By how much?

from decimal import Decimal
from fractions import Fraction

print(Decimal(0.1))
print(Fraction(0.1))

Output: 0.1000000000000000055511151231257827021181583404541015625 and 3602879701896397/36028797018963968. The float 0.1 is actually this fraction, whose denominator is 36028797018963968 = 2⁵⁵: a power of 2, as it must be. This is why 0.1 + 0.2 does not give exactly 0.3: you add two numbers that have already been rounded, and the sum is rounded again.

Not associative: the order of the sums matters

A surprising fact that follows from rounding: addition between floats is not associative. In mathematics (a + b) + c = a + (b + c); in floating point it is not, because every sum rounds the intermediate result.

import math
import numpy as np

a, b, c = 1e16, -1e16, 1.0
print("(a + b) + c =", (a + b) + c)
print("a + (b + c) =", a + (b + c))

rng = np.random.default_rng(0)
x = rng.random(1_000_000)
ordinamenti = {
    "original order": x,
    "ascending": np.sort(x),
    "descending": np.sort(x)[::-1],
}
esatta = math.fsum(x)                                    # accurate sum, a reliable reference
for nome, v in ordinamenti.items():
    s = float(np.sum(v))
    print(f"{nome:<18}: {s!r}  difference from the accurate sum: {s - esatta:+.2e}")

Output obtained:

(a + b) + c = 1.0
a + (b + c) = 0.0
original order    : 500159.2564636844  difference from the accurate sum: -5.82e-11
ascending         : 500159.25646368443  difference from the accurate sum: +0.00e+00
descending        : 500159.25646368443  difference from the accurate sum: +0.00e+00

How to read the code. The first part shows the extreme case: 1e16 + (-1e16) is 0 and then 1 is added (result 1), but adding first -1e16 + 1 gives -1e16 again (the one is lost because floats around 10¹⁶ are 2 apart) and the final result is 0. Same numbers, same set of operations, different results.

The second part sums a million random numbers in the three orders; math.fsum provides a reference obtained with a procedure that keeps track of the partial errors (it returns the float closest to the exact sum). The differences are tiny: the original order departs from the accurate sum by about 6·10⁻¹¹, while the two sorted orders coincide with it; with another array, or another seed, the situation would be different. What matters is that the result depends on the order, even if by very little. This is why two runs of the same training on a GPU, which sum in different orders, can give results that differ in the last digits, and why an exact-equality test between two "equivalent" computations is never reliable.

Machine epsilon and spacing

Between two consecutive powers of 2, 2^k and 2^(k+1), there are exactly 2⁵² equally spaced float numbers (the possible values of the fraction). So the distance between two neighbouring floats in that interval is

ulp(x) = 2^(k−52) for 2^k ≤ x < 2^(k+1)

(ulp = unit in the last place). In particular, above 1.0: ulp(1) = 2⁻⁵² ≈ 2.22 × 10⁻¹⁶. This number is called the machine epsilon.

Rounding to nearest, the maximum error is half the spacing, so the relative error of a single operation is at most u = 2⁻⁵³ ≈ 1.11 × 10⁻¹⁶: about 16 decimal digits of precision.

import math
import sys

print(sys.float_info.epsilon, 2**-52, math.ulp(1.0))
print(math.ulp(1e16), math.ulp(0.0))
print(2.0**53 == 2.0**53 + 1)

Output: three times 2.220446049250313e-16; then 2.0 and 5e-324; then True.

Interpret the last three lines:

• math.ulp(1e16) is 2.0: around 10¹⁶ floats are 2 apart. So 1e16 + 1 is not representable and is rounded.
• math.ulp(0.0) is 5e-324, the smallest positive float (denormalised).
• 2.0**53 == 2.0**53 + 1 is True: from 2⁵³ onwards floats no longer represent all the integers. This is why large integers (IDs, timestamps in nanoseconds) must never be converted to float.

The relative spacing stays roughly constant (between 2⁻⁵³ and 2⁻⁵²); the absolute one grows with the number. Floats are very precise near zero and less and less precise far from it.

A second, contrasting example: float32 and float16

float64 gives about 16 decimal digits. Neural networks often use float32 (about 7 digits) or float16 (about 3). Let us see what is lost, with two experiments.

import numpy as np

# 1) accumulate 20 million ones, in float32
uni = np.ones(20_000_000, dtype=np.float32)
print("cumsum in float32 (sequential):", float(np.cumsum(uni)[-1]))
print("np.sum in float32 (pairwise)  :", float(uni.sum()))
print("np.sum with float64 accumulator:", float(uni.sum(dtype=np.float64)))

# 2) the limits of float16
with np.errstate(over="ignore"):
    print("float16: 60000 * 2 =", np.float16(60000) * np.float16(2))
print("float16: 0.1 becomes", float(np.float16(0.1)), "| float32: 0.1 becomes", float(np.float32(0.1)))
print("float32: 16777216 + 1 =", float(np.float32(16777216) + np.float32(1)))

Output obtained:

cumsum in float32 (sequential): 16777216.0
np.sum in float32 (pairwise)  : 20000000.0
np.sum with float64 accumulator: 20000000.0
float16: 60000 * 2 = inf
float16: 0.1 becomes 0.0999755859375 | float32: 0.1 becomes 0.10000000149011612
float32: 16777216 + 1 = 16777216.0

How to read the code. In the first part, np.cumsum sums in sequence: when the partial total reaches 2²⁴ = 16,777,216, in float32 the spacing between consecutive numbers becomes 2, so adding 1 no longer changes the total, which stays stuck at 16,777,216 instead of reaching 20 million. np.sum, on the other hand, internally uses a pairwise sum (it sums small groups and then combines them), which keeps the addends of similar size and produces the correct result of 20 million: the summation algorithm matters as much as the precision of the type. With dtype=np.float64 the accumulator has an enormous margin. In the second part, float16 cannot represent 120000 (the maximum is 65504) and returns inf; it represents 0.1 as 0.09997... (an error of 0.02%), while float32 as 0.10000000149.... The practical lesson for deep learning: reduced precisions speed things up and save memory but must be used with more precise accumulators and with attention to the limits of the representable range (hence loss scaling, defined in the list further down).

Three cases where the algorithm breaks

Case 1: catastrophic cancellation

Subtracting two nearly equal numbers wipes out the significant digits they share and leaves only the "noise" of the least significant bits.

print((1e16 + 1) - 1e16)
print((1e16 + 2) - 1e16)

Output: 0.0 and 2.0. The first should be 1: 1e16 + 1 was rounded to 1e16, and the difference is zero. The second is 2 because 1e16 + 2 is representable.

Application: the "textbook" variance

The convenient formula Var(X) = E[X²] − E[X]² is algebraically correct (you derive it from the properties of the expected value). In floating point, if the data have a large mean relative to their spread, the two terms are enormous and almost equal, and the difference is dominated by noise.

import numpy as np

rng = np.random.default_rng(0)
x = rng.normal(0, 1, 1000) + 1e9          # spread ~1, mean ~1 billion

var_ingenua = np.mean(x**2) - np.mean(x) ** 2
var_due_passi = np.mean((x - x.mean()) ** 2)
print(var_ingenua, var_due_passi, np.var(x))

Output: -128.0, 0.9540463497672959, 0.9540463497672959. The naive formula returns a negative variance: E[X²] is about 10¹⁸ and at that scale floats are 128 apart (math.ulp(1e18)), so the difference is dominated by rounding. The two-pass method (first subtract the mean, then square) subtracts small numbers and works.

Engineering lesson: two formulas equal on paper are not equal on the computer. NumPy uses the two-pass method.

Case 2: summing many numbers

import math

import numpy as np

def somma_a_mano(valori):
    s = 0.0
    for v in valori:
        s += v
    return s

print(somma_a_mano([0.1] * 10))
print(somma_a_mano([0.1] * 1_000_000))
print(sum([0.1] * 1_000_000), math.fsum([0.1] * 1_000_000))
print(np.sum(np.full(1_000_000, 0.1)))

Output: 0.9999999999999999, 100000.00000133288, 100000.0 100000.0, 100000.00000000003.

The naive loop accumulates error: every sum rounds, and the errors add up (100000.00000133288 instead of 100000). There are three alternatives. math.fsum sums with extended precision and returns the correctly rounded result: it costs more, but it is the referee when you need a reference. Python's sum() since version 3.12 uses a compensated summation (Neumaier's algorithm) for floats, as the 3.12 release notes say, and this is where the 100000.0 comes from; before that it behaves like the hand-written loop.

np.sum uses a partial pairwise summation: according to the documentation of numpy.sum it is always applied when you do not indicate an axis, while along an axis it depends on which one. It sums in blocks recursively and reduces the growth of the error, in the worst case, from proportional to n to proportional to log n. Here it does not give exactly 100000, but the error (about 3e-11) is more than four orders of magnitude smaller than that of the loop (about 1.3e-6).

Kahan's algorithm (1965) keeps a compensation variable that collects the rounding error of each step and reintroduces it at the next:

def kahan(valori):
    somma, compensazione = 0.0, 0.0
    for v in valori:
        y = v - compensazione
        t = somma + y
        compensazione = (t - somma) - y     # what was lost in the addition
        somma = t
    return somma

print(kahan([0.1] * 1_000_000))

Output: 100000.0. Careful: an over-"clever" compiler can algebraically simplify (t − somma) − y to zero and disable the compensation; in Python it does not happen.

Trade-off. Fsum and Kahan are slower than the simple sum. For np.sum on ordinary arrays pairwise summation is already an excellent compromise. Special care is needed only with very many terms, values of very different scale, or where the exact result matters.

Case 3: softmax

The softmax function turns a vector of scores z into probabilities, and sits at the output of almost every multi-class classifier:

softmax(z)ᵢ = e^(zᵢ) / Σⱼ e^(zⱼ)

With large scores, e^z overflows: np.exp(710) is inf (the maximum float64 is about 1.8 × 10³⁰⁸, and e^709 ≈ 8.2 × 10³⁰⁷ is still representable). Then inf/inf = nan.

Remedy and proof of correctness. For any constant c:

e^(zᵢ − c) / Σⱼ e^(zⱼ − c) = ( e^(zᵢ) · e^(−c) ) / ( e^(−c) · Σⱼ e^(zⱼ) ) = e^(zᵢ) / Σⱼ e^(zⱼ)

because e^(−c) is a common factor and cancels. The result does not change. Choosing c = max(z) the largest exponent becomes e⁰ = 1 and no term exceeds 1: no overflow.

import numpy as np

def softmax_ingenuo(z):
    e = np.exp(z)
    return e / e.sum()

def softmax_stabile(z):
    e = np.exp(z - z.max())
    return e / e.sum()

z = np.array([1000.0, 1001.0, 1002.0])
with np.errstate(all="ignore"):
    print(softmax_ingenuo(z), softmax_stabile(z))
piccolo = np.array([1.0, 2.0, 3.0])
print(softmax_ingenuo(piccolo), softmax_stabile(piccolo))

Output: [nan nan nan] versus [0.09003057 0.24472847 0.66524096]; on the small example the two coincide. This trick (and its relative for the logarithm, log-sum-exp) is in the code of every deep learning library. If one day a network of yours returns nan, suspect first of all an overflow or a division by zero.

How to defend yourself

Comparing with a tolerance

math.isclose(a, b, rel_tol=1e-9, abs_tol=0.0) returns true if

|a − b| ≤ max( rel_tol · max(|a|, |b|), abs_tol )

The relative tolerance is fine for numbers far from zero; near zero it is useless (rel_tol · 0 = 0), so an absolute tolerance is needed:

import math
import numpy as np

print(math.isclose(1e-12, 2e-12))                  # False: relatively, one is double the other
print(math.isclose(1e-12, 2e-12, abs_tol=1e-9))    # True: in absolute terms they are very close
print(np.allclose([1.0, 2.0], [1.0 + 1e-9, 2.0]))  # True

Output: False, True, True. For arrays: np.isclose and np.allclose. In tests: np.testing.assert_allclose.

Money and exact decimals

from decimal import Decimal

print(Decimal("0.1") + Decimal("0.2"), Decimal("0.1") + Decimal("0.2") == Decimal("0.3"))

Output: 0.3 True. Decimal reasons in base 10, so 0.1 is exact. Careful: build from a string, not from a float (Decimal(0.1) inherits the error, as seen above). In practice, for money, using integers in cents is often enough.

Rounding

Python's round() rounds .5 to the nearest even (round half to even), not always up: round(0.5) = 0, round(1.5) = 2, round(2.5) = 2, round(3.5) = 4. It is the standard's choice, which avoids a systematic drift in sums of many roundings. Also, round(2.675, 2) is 2.67: because the float near 2.675 is actually slightly below 2.675.

32- and 16-bit floats: machine learning lives here

(Values obtainable from np.finfo(type); np.finfo(np.float16).max is 65504.)

Neural networks are often trained in float32 or in 16-bit formats, because memory and speed matter more than precision: the weights are noisy anyway. Consequences to know:

• With float32, 16777216 + 1 == 16777216 (at 2²⁴ float32 can no longer tell consecutive integers apart). Accumulators are often kept at higher precision.
• In float16 the maximum is 65504: a larger value goes to inf. This is where the technique of loss scaling in mixed-precision training comes from.
• The behaviour changes with the order of the operations: training the same network twice on a GPU can give slightly different results because the parallel sums are in different orders (and floating-point addition is not associative: (a + b) + c ≠ a + (b + c) in general).

The three consequences have the same root: fewer bits mean a wider gap between numbers and a lower maximum. The remedy is not "always use float64", which costs memory and speed, but to choose where to keep the precision (in the accumulator and not in the weights) and where to protect the range of the values.

Special values. inf (overflow), -inf, nan (undefined result: 0.0 * inf). nan is not equal to anything, not even to itself: float("nan") == float("nan") is False. To test for it use math.isnan or np.isnan. A nan propagates: every operation with a nan gives nan, so a single corrupted value can pollute an entire metric.

A small numerical debugging protocol

When a computation gives a strange result (a nan, a negative variance, a metric that changes at every run), it pays to follow an order of investigation, from the cheapest to the most expensive. First: find where the first anomalous value appears (np.isnan(x).argmax() gives the position of the first nan in an array); the error is born there, everything else is propagation. Second: print the orders of magnitude of the quantities involved (minima, maxima, np.log10(np.abs(x).max())): a subtraction between two numbers with similar orders of magnitude, much larger than the result, is a candidate for cancellation.

Third: redo the same computation in a more precise type (float64 instead of float32, or Decimal/Fraction for small examples): if the result changes visibly, the problem is numerical and not logical. Fourth: rewrite the formula in a stable form (two passes, subtracting the maximum, log1p, fsum) and check that it matches the naive version on small data. Fifth: replace every == between floats with a declared tolerance. Always following the same order avoids losing hours hunting a logic bug that is not there.

Typical mistakes

Best practices

• Never == between computed floats. Always a declared tolerance, with a justification.
• Normalise before subtracting large values; prefer numerically stable formulas (two passes, log-sum-exp).
• Use library functions (np.var, np.sum, np.logaddexp) instead of rewriting formulas: they have been studied for stability.
• Declare the type (dtype) when it matters; do not trust the defaults for sensitive data.
• Check for nan (np.isnan(x).any()) at the input and after critical steps.
• Money: integers or Decimal.

The six rules descend from the same fact: a float is not the real number, it is its nearest representative. Library functions, declared tolerances and nan checks are three ways of not trusting that representative without having looked at it. Decide the policy once per project and write it in the README, so that it does not stay in the head of one person alone.

Mini-project

Build a small numerical laboratory lab_float.py:

1. A function ulp_attorno(x) that compares math.ulp(x) with the formula 2**(floor(log2(x)) - 52) for x in [1e-5, 1e20] and checks that they coincide.
2. A function somma_confronto(n) that returns the absolute errors (with respect to math.fsum) of the naive loop, np.sum, Kahan for n random values.
3. A comparison of the naive and stable variance as the offset 10**k varies, k = 0..12, with a chart of the relative error as a function of k.
4. A test showing that (a + b) + c != a + (b + c) for three values of your choice.

The first point checks the ulp formula against the library function, the second and third repeat the two cases of the lesson on your own data, the fourth forces you to find three numbers that violate associativity. For the chart in the third point you need Matplotlib, covered in visualisation with Matplotlib. There is no official solution: the criterion is that each function prints a number you had predicted before running it.

Challenge

• Easy. Find by hand the binary representation of 0.375 and check with struct that the fraction terminates.
• Medium. Write a stable log_somma_exp(z): log Σ e^(zᵢ) computed as c + log Σ e^(zᵢ − c) with c = max(z). Prove the identity.
• Hard. Implement Welford's algorithm for streaming variance and compare it with the naive one on the same data challenge with offset 10⁹.
• Expert. Show that floating-point addition is not associative with a minimal example found by random search. Then compute, for an array of n numbers, the upper bound of the error of recursive summation and of pairwise summation (Higham, ch. 4).
• Enterprise. Define a numerical policy for your team: when to use float32 and when float64, which tolerances in tests, how to detect a nan in production.

The five levels rise from a check by hand to a team decision. Easy and Medium are solved with what you have seen here; Hard and Expert ask you to consult the literature cited at the bottom; Enterprise has no code answer but a one-page document your colleagues can apply.

What you can do now

You can read the three fields of a float64 and rebuild the value with the format's formula, you can explain why 0.1 has no finite binary expansion (the denominator has a factor 5) and estimate with the ulp the distance between two neighbouring floats. You recognise the three points where an algorithm that is correct on paper breaks, namely cancellation, the sum of many terms and overflow in the exponential, and you know which remedy to use in each. The next step is NumPy: arrays, vectorisation and broadcasting, which brings all this to millions of numbers at a time.

The same series also includes: 13 - Pandas: Data Tables Without Surprises, 21 - Calculus 1: Limits and Derivatives, 25 - Matrix Calculus and Automatic Differentiation: From the Jacobian to an Autograd Written From Scratch, 27 - Entropy, Cross-Entropy and KL Divergence: Why Loss Functions Are All the Same Thing (which refer back to this lesson).

Sources

Bibliographic references cited in the text, not consulted online in this revision: IEEE Std 754-2019; D. Goldberg, What Every Computer Scientist Should Know About Floating-Point Arithmetic, ACM Computing Surveys 23(1), 1991; N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM, 2002; W. Kahan, Communications of the ACM 8(1), 1965. The measurements in this article (ulp, epsilon, outputs) were recomputed with Python 3.14.4 and NumPy 2.3.5.