Floating point for scientists and engineers

In Julia 0.1 + 0.2 is 0.30000000000000004, so a test that checks it against 0.3 with == fails, and the formula (1 − cos x)/x², whose true value is very close to 0.5 for small x, returns 0 at x = 10⁻⁹. Both come from the way a Float64 stores numbers, and both have exact remedies.

Each section below is a run recorded on Julia 1.13.1: comparing results with isapprox and a tolerance that suits them, the grid of floats that eps measures, NaN, signed zeros, Inf and subnormal numbers, the rounding modes, and rewriting formulas and sums that lose digits. One rounding behaviour changed in 1.13, so its runs on 1.12 and on the 1.10 LTS are shown too.

Two floats agree within a tolerance, not exactly

A Float64 holds a binary fraction, and 0.1, 0.2 and 0.3 are not binary fractions, so each is stored as the nearest one and the sum lands a little off. Here is the comparison, then the same question asked with a tolerance:

x = 0.1 + 0.2
println("0.1 + 0.2 = ", x, ", minus 0.3: ", x - 0.3)
println("x == 0.3: ", x == 0.3, ", x ≈ 0.3: ", x ≈ 0.3)

# isapprox(x, y) asks whether |x - y| <= max(atol, rtol * max(|x|, |y|)).
println("default rtol for Float64, √eps: ", sqrt(eps(Float64)))
println("isapprox(1e6, 1e6 + 0.01): ", isapprox(1e6, 1e6 + 0.01))
println("isapprox(1.0, 1.01): ", isapprox(1.0, 1.01), ", with rtol = 0.02: ", isapprox(1.0, 1.01; rtol = 0.02))

# Next to zero a relative tolerance is almost nothing: give an absolute one.
println("isapprox(0.0, 1e-12): ", isapprox(0.0, 1e-12))
println("isapprox(0.0, 1e-12; atol = 1e-10): ", isapprox(0.0, 1e-12; atol = 1e-10))

from session import files, run

# == and ≈ on 0.1 + 0.2, then isapprox's default tolerance at three sizes and next to zero.
files("compare.jl")
run("""
julia compare.jl
""")

$ julia compare.jl
0.1 + 0.2 = 0.30000000000000004, minus 0.3: 5.551115123125783e-17
x == 0.3: false, x ≈ 0.3: true
default rtol for Float64, √eps: 1.4901161193847656e-8
isapprox(1e6, 1e6 + 0.01): true
isapprox(1.0, 1.01): false, with rtol = 0.02: true
isapprox(0.0, 1e-12): false
isapprox(0.0, 1e-12; atol = 1e-10): true

The sum missed 0.3 by 5.55e-17, so == said false and ≈, which is isapprox, said true. Its documentation gives the rule: x ≈ y when |x − y| ≤ max(atol, rtol × max(|x|, |y|)). The relative tolerancerelative tolerance: The largest difference allowed between two numbers as a share of the larger one: isapprox's rtol. Its default for Float64 is √eps, about 1.5e-8, so about half the digits must agree., rtol, defaults to √eps, 1.49e-8, which asks for about half of a Float64’s digits to agree; the absolute toleranceabsolute tolerance: The largest difference allowed between two numbers in their own units, whatever their size: isapprox's atol. It is 0 by default and is needed whenever a value may be zero., atol, defaults to 0. Step through the three comparisons of the run:

A relative tolerance scales with the numbers, which is right for a result whose size you do not know in advance: 1e6 and 1e6 + 0.01 agreed to within a part in 10⁸, so they were approximately equal, while 1.0 and 1.01 were not, until rtol = 0.02 allowed 2 %. Next to zero the scale is almost nothing, and the documentation says so plainly: with the default tolerances x ≈ 0 is equivalent to x == 0. A result that should be zero needs an atol in the units of the problem, and no value suits every problem: 1e-10 is tiny for a distance in metres between cities and enormous for one between atoms.

A float is a grid, and eps measures its spacing

A Float64 is a binary64 number of IEEE 754: a sign, 11 bits of exponent and 52 bits of fraction. Between two powers of two there are 2⁵² evenly spaced values, so the spacing doubles at every power of two. eps(x) gives the gap from x to the next larger value:

# Floats are a grid whose spacing grows with the size of the number: eps(x) is the gap above x.
for x in (1.0, 1000.0, 1.0e6, 1.0e16)
    println(rpad(x, 8), "  eps(x) = ", eps(x))
end
println("nextfloat(1.0) = ", nextfloat(1.0), ", prevfloat(1.0) = ", prevfloat(1.0))
println("1e16 + 1 == 1e16: ", 1e16 + 1 == 1e16, ", 2.0^53 + 1 == 2.0^53: ", 2.0^53 + 1 == 2.0^53)

# What is really stored for 0.1: the bits, and the exact value of those bits.
println("bitstring(0.1) = ", bitstring(0.1))
println("big(0.1)       = ", big(0.1))
println("eps(Float32) = ", eps(Float32), ", eps(Float16) = ", eps(Float16))

from session import files, run

# The spacing of Float64 values at four sizes, the neighbours of 1.0, and what is really stored for 0.1.
files("grid.jl")
run("""
julia grid.jl
""")

$ julia grid.jl
1.0       eps(x) = 2.220446049250313e-16
1000.0    eps(x) = 1.1368683772161603e-13
1.0e6     eps(x) = 1.1641532182693481e-10
1.0e16    eps(x) = 2.0
nextfloat(1.0) = 1.0000000000000002, prevfloat(1.0) = 0.9999999999999999
1e16 + 1 == 1e16: true, 2.0^53 + 1 == 2.0^53: true
bitstring(0.1) = 0011111110111001100110011001100110011001100110011001100110011010
big(0.1)       = 0.1000000000000000055511151231257827021181583404541015625
eps(Float32) = 1.1920929e-7, eps(Float16) = 0.000977

At 1.0 the gap is 2.22e-16, which is machine epsilonmachine epsilon: The gap between 1.0 and the next larger floating-point number: eps() is 2.220446049250313e-16 for Float64. eps(x) gives the gap above x, which doubles at every power of two., eps(); at 1000 it is 1.14e-13 and at 10¹⁶ it is 2.0, so 1e16 + 1 rounds straight back to 1e16. From 2⁵³ on the gap is 2 or more, and whole numbers start to go missing: 2.0^53 + 1 == 2.0^53 is true. Measured against the number itself, the gap stays between 1.1e-16 and 2.2e-16 at every magnitude, which is why a Float64 keeps about 15 to 16 significant decimal digits wherever it is.

The last lines show what 0.1 really is. Its 64 bits end in a repeating 1001 pattern, cut off and rounded, and big(0.1) prints the exact value of those bits, 0.1000000000000000055511151231257827…, a little more than 0.1. That tiny excess is where the 0.30000000000000004 of the first run came from. Float32 has 23 bits of fraction, so its epsilon is 1.19e-7, and Float16’s is 0.000977.

NaN, signed zero and Inf follow IEEE 754, and three equalities differ on them

IEEE 754 defines values for results that have no ordinary number: Inf for overflow and division of a non-zero number by zero, NaN, “not a number”, for undefined results such as 0/0, and two zeros, 0.0 and −0.0. Julia has three ways to ask whether two values are equal, and they disagree exactly on these values:

# Three kinds of equality on the values where they disagree.
println(rpad("x, y", 11), rpad("x == y", 8), rpad("isequal", 9), "x === y")
for (label, x, y) in (("NaN, NaN", NaN, NaN), ("-0.0, 0.0", -0.0, 0.0), ("1, 1.0", 1, 1.0))
    println(rpad(label, 11), rpad(x == y, 8), rpad(isequal(x, y), 9), x === y)
end
println("[1, NaN] == [1, NaN]: ", [1, NaN] == [1, NaN], ", isequal: ", isequal([1, NaN], [1, NaN]))

# Dictionaries and sorting use isequal and isless, so -0.0 and 0.0 are two keys, and NaN sorts last.
d = Dict(0.0 => "zero", -0.0 => "minus zero", NaN => "not a number")
println(length(d), " keys; d[NaN] = ", d[NaN])
println(sort([3.0, NaN, -0.0, 0.0, -Inf, 1.0]))

# Where Inf, NaN and the signed zeros come from, and the smallest numbers a Float64 holds.
println("1/0 = ", 1/0, ", -1/0 = ", -1/0, ", 0/0 = ", 0/0, ", Inf - Inf = ", Inf - Inf, ", 1/-0.0 = ", 1/-0.0)
println("floatmin() = ", floatmin(), ", nextfloat(0.0) = ", nextfloat(0.0))
println("issubnormal(1e-320) = ", issubnormal(1e-320), ", 1e-320 * 1e10 = ", 1e-320 * 1e10)

from session import files, run

# ==, isequal and === on NaN, the two zeros and 1 against 1.0; what a Dict and sort make of them; special values.
files("equality_table.jl")
run("""
julia equality_table.jl
""")

$ julia equality_table.jl
x, y       x == y  isequal  x === y
NaN, NaN   false   true     true
-0.0, 0.0  true    false    false
1, 1.0     true    true     false
[1, NaN] == [1, NaN]: false, isequal: true
3 keys; d[NaN] = not a number
[-Inf, -0.0, 0.0, 1.0, 3.0, NaN]
1/0 = Inf, -1/0 = -Inf, 0/0 = NaN, Inf - Inf = NaN, 1/-0.0 = -Inf
floatmin() = 2.2250738585072014e-308, nextfloat(0.0) = 5.0e-324
issubnormal(1e-320) = true, 1e-320 * 1e10 = 9.999888671827e-311

== follows IEEE 754: NaNNaN: Not a number: the floating-point result of an undefined operation such as 0/0 or Inf - Inf. NaN is not == to anything, itself included, so it is tested with isnan. is not equal to anything, itself included, and the two signed zerossigned zero: Either of the two zeros of IEEE 754 floating point, 0.0 and -0.0. They are == to each other but not isequal, and 1/-0.0 is -Inf while 1/0.0 is Inf. are equal. isequal treats every NaN as equal to every NaN and tells −0.0 from 0.0, and === asks whether no program could tell the two apart, which is false even for 1 and 1.0, an Int64 and a Float64. Dictionaries use isequal, so 0.0 and -0.0 became two keys and d[NaN] found its value; sorting uses isless, which put −0.0 before 0.0 and NaN last.

The last lines show where the special values come from and how small a Float64 can be. 1/0 is Inf and 1/-0.0 is −Inf, so the sign of a zero can decide a result. floatmin(), 2.2e-308, is the smallest number with all 53 bits of precision; below it come the subnormal numberssubnormal number: A floating-point number smaller in size than floatmin(), stored with fewer significant bits; for Float64, between 5.0e-324 and 2.2250738585072014e-308. issubnormal tests for one., down to 5.0e-324, which keep fewer bits as they shrink. 1e-320 is one, and multiplying it by 10¹⁰ gave 9.999888671827e-311 instead of 1e-310: most of its digits were already gone.

Rounding is a choice of mode, and ties go to even

round rounds to the nearest integer, and a value exactly halfway between two of them is a tie. Julia’s default mode, RoundNearest, sends ties to the even neighbour, as IEEE 754 does by default; other modes are one argument away:

# round sends ties to the even neighbour; other modes and the floor, ceil and trunc functions round as they say.
println(rpad("x", 6), rpad("round", 7), rpad("TiesAway", 10), rpad("floor", 7), rpad("ceil", 6), "trunc")
for x in (0.5, 1.5, 2.5, -2.5)
    println(rpad(x, 6), rpad(round(x), 7), rpad(round(x, RoundNearestTiesAway), 10),
            rpad(floor(x), 7), rpad(ceil(x), 6), trunc(x))
end
println("round(π; digits = 3) = ", round(π; digits = 3))
println("round(123456; sigdigits = 2) = ", round(123456; sigdigits = 2))
println("round(0.125; digits = 2) = ", round(0.125; digits = 2), ", round(1.15; digits = 1) = ", round(1.15; digits = 1))
println("big(1.15) = ", big(1.15))
println("round(Int, 2.5) = ", round(Int, 2.5), ", round(Int, 3.5) = ", round(Int, 3.5))

from session import files, run, use

# The rounding modes on ties, rounding to digits, then one behaviour that changed in 1.13: mod with an infinite divisor.
files("rounding.jl")
run("""
julia rounding.jl
julia -e 'println(VERSION, ": mod(5.5, -Inf) = ", mod(5.5, -Inf))'
""")
use("1.12.7")
run("""
julia -e 'println(VERSION, ": mod(5.5, -Inf) = ", mod(5.5, -Inf))'
""")
use("1.10.12")
run("""
julia -e 'println(VERSION, ": mod(5.5, -Inf) = ", mod(5.5, -Inf))'
""")

$ julia rounding.jl
x     round  TiesAway  floor  ceil  trunc
0.5   0.0    1.0       0.0    1.0   0.0
1.5   2.0    2.0       1.0    2.0   1.0
2.5   2.0    3.0       2.0    3.0   2.0
-2.5  -2.0   -3.0      -3.0   -2.0  -2.0
round(π; digits = 3) = 3.142
round(123456; sigdigits = 2) = 120000.0
round(0.125; digits = 2) = 0.12, round(1.15; digits = 1) = 1.2
big(1.15) = 1.149999999999999911182158029987476766109466552734375
round(Int, 2.5) = 2, round(Int, 3.5) = 4
$ julia -e 'println(VERSION, ": mod(5.5, -Inf) = ", mod(5.5, -Inf))'
1.13.1: mod(5.5, -Inf) = 5.5
$ julia -e 'println(VERSION, ": mod(5.5, -Inf) = ", mod(5.5, -Inf))'
1.12.7: mod(5.5, -Inf) = -Inf
$ julia -e 'println(VERSION, ": mod(5.5, -Inf) = ", mod(5.5, -Inf))'
1.10.12: mod(5.5, -Inf) = -Inf

round(1.5) and round(2.5) both gave 2.0, and round(0.5) gave 0.0: ties to eventies to even: The default rule for a value exactly halfway between two candidates: choose the even one, so round(2.5) is 2.0 and round(3.5) is 4.0. Julia calls the mode RoundNearest.. Rounding every tie up would push a long column of sums upwards on average; sending half of them down keeps the average where it was. RoundNearestTiesAway rounds ties away from zero, as C’s round does, so 2.5 gave 3.0 and −2.5 gave −3.0, and floor, ceil and trunc round down, up and toward zero whatever the fraction. round(Int, 2.5) applies the same rule and returns an integer, 2.

Rounding to decimal places meets the binary grid. round(0.125; digits = 2) gave 0.12, because 0.125 is exactly representable and is a true tie. round(1.15; digits = 1) gave 1.2, although big(1.15) shows that the stored value is a little less than 1.15; the documentation warns that rounding to digits in base 10 can be inexact for exactly this reason. Keep full precision while computing and round once, when you present a result.

Cancellation and summation order lose digits, and rewrites win them back

Subtracting two numbers that agree in most of their digits leaves only the digits where they differ, and those are the ones rounding has already touched: catastrophic cancellationcatastrophic cancellation: The loss of significant digits when two nearly equal numbers are subtracted: their leading digits cancel, and what is left is mostly the rounding error of each.. 1 − cos x does that for small x. This run compares the formula with an algebraically equal rewrite, using the identity 1 − cos x = 2 sin²(x/2), and measures both against a BigFloat of 256 bits; then it adds two sets of numbers in different orders:

using Printf
naive(x) = (1 - cos(x)) / x^2         # 1 - cos(x) subtracts two numbers that are almost equal
stable(x) = 2 * sin(x / 2)^2 / x^2    # the same function, with 1 - cos(x) written as 2 sin(x/2)^2

@printf("%-7s %-22s %-22s %s\n", "x", "naive", "stable", "errors: naive, stable")
for x in (1e-1, 1e-3, 1e-5, 1e-7, 1e-9)
    exact = naive(big(x))               # BigFloat, 256 bits by default: a reference to compare with
    err(v) = Float64(abs(v - exact) / exact)
    @printf("%-7.0e %-22.17g %-22.17g %.1e, %.1e\n", x, naive(x), stable(x), err(naive(x)), err(stable(x)))
end
println("precision(BigFloat) = ", precision(BigFloat))

# The same numbers added in different orders, against an exact BigFloat sum.
function left_to_right(xs)
    s = 0.0
    for x in xs
        s += x
    end
    return s
end
right_to_left(xs) = left_to_right(reverse(xs))

"Kahan's compensated sum: c keeps the low bits that each addition rounds away."
function kahan(xs)
    s = 0.0
    c = 0.0
    for x in xs
        y = x - c          # put back what the last addition lost
        t = s + y          # the big addition, which rounds
        c = (t - s) - y    # what that rounding lost, as a small number
        s = t
    end
    return s
end

small = [1.0; fill(1e-16, 10)]           # 1 followed by ten numbers below half the gap above 1
println("exact:         ", sum(big.(small)))
for (label, f) in (("left to right", left_to_right), ("right to left", right_to_left), ("Kahan", kahan))
    println(rpad(label * ":", 15), f(small))
end

tenths = fill(0.1, 10^6)                 # a million copies of the stored 0.1
exact = sum(big.(tenths))
println("exact - 100000: ", Float64(exact - 100000))
for (label, f) in (("left to right", left_to_right), ("sum (pairwise)", sum), ("Kahan", kahan))
    println(rpad(label * ":", 16), f(tenths), "   error ", Float64(f(tenths) - exact))
end

from session import files, run

# A formula that cancels and its rewrite, then two sets of numbers summed four ways.
files("cancellation.jl", "summation.jl")
run("""
julia cancellation.jl
julia summation.jl
""")

$ julia cancellation.jl
x       naive                  stable                 errors: naive, stable
1e-01   0.49958347219741783    0.49958347219742333    1.1e-14, 1.2e-16
1e-03   0.49999995832550326    0.49999995833333472    1.6e-11, 1.4e-17
1e-05   0.50000004137018539    0.49999999999583328    8.3e-08, 1.1e-16
1e-07   0.4996003610813205     0.49999999999999961    8.0e-04, 5.6e-17
1e-09   0                      0.5                    1.0e+00, 8.3e-20
precision(BigFloat) = 256
$ julia summation.jl
exact:         1.000000000000000999999999999999979097786724034603561841114940846736436341757326
left to right: 1.0
right to left: 1.000000000000001
Kahan:         1.000000000000001
exact - 100000: 5.551115123125783e-12
left to right:  100000.00000133288   error 1.3328771242271387e-6
sum (pairwise): 99999.99999999974   error -2.674855892337291e-10
Kahan:          100000.0   error -5.551115123125783e-12

The naive formula lost digits steadily as x shrank: a relative error of 1.1e-14 at x = 0.1, 8.3e-8 at 10⁻⁵ and 8.0e-4 at 10⁻⁷, and at 10⁻⁹ cos(x) rounded to exactly 1.0 and the result was 0. The rewrite subtracts nothing and stayed within 1.2e-16 everywhere. The same idea is behind functions such as log1p and expm1, which the next lesson uses.

Sums lose digits too. After 1.0, each 1e-16 is less than half the gap above 1.0, so adding them one at a time left 1.0, while adding the small ones first gave 1.000000000000001. Kahan’s compensated summation keeps the lost part in a second variable and adds it back; step through its first two passes:

For a million copies of 0.1, a plain loop was off by 1.3e-6, Julia’s sum by 2.7e-10 and Kahan’s sum by 5.6e-12, the smallest error a Float64 result can have, because the exact sum is not a Float64. sum does better than the loop because, as the source of Julia 1.13.1 shows, it adds an array by pairwise summationpairwise summation: Adding numbers by splitting them into halves, adding each half the same way, and adding the two results, so rounding errors grow far more slowly than in one long loop. Julia's sum of an array works this way.: it splits the array in halves until the pieces have 1,024 numbers or fewer, adds each piece, and then the halves. Its documentation promises no particular order, so results can differ in the last bits from a loop’s.

Practise: floats as they really behave

Every question here is settled by a run on this page. Your weakest idea comes up first, and a wrong choice tells you which slip it is:

Interview questions

Warm-up (fresher to 2–5 years): why is 0.1 + 0.2 == 0.3 false in Julia, and how should you compare floating-point results?

Julia keeps a Float64 as a binary fraction, and none of 0.1, 0.2 and 0.3 has an exact binary form, so each becomes the nearest Float64. The stored 0.1 and 0.2 are both slightly too large, their sum rounds to 0.30000000000000004, and the nearest Float64 to 0.3 is a different value, so == is false. Compare with isapprox, written ≈: it accepts a gap up to max(atol, rtol × max(|x|, |y|)), with a default rtol of √eps, about 1.5e-8, and atol 0. For results that should be zero, give an absolute tolerance, isapprox(x, 0.0; atol = 1e-12), chosen from the scale of the problem, because with the defaults x ≈ 0 is the same as x == 0. Use == only for values that are exact by construction, such as small integers stored as floats.

Two harder questions follow, about a test against zero and about a variance that comes out negative. Note your answer first, then open each first step:

Key takeaways

• Compare floats with isapprox (≈), whose default relative tolerance is √eps and absolute tolerance 0, and give an atol in the problem’s units whenever a result may be zero.
• Floats form a grid whose gap, eps(x), doubles at every power of two: 2.22e-16 at 1.0 and 2.0 at 10¹⁶, so a Float64 keeps about 15 to 16 significant digits and loses whole numbers above 2⁵³.
• == follows IEEE 754 (NaN is never equal, the two zeros are), isequal and dictionaries treat NaN as equal and tell −0.0 from 0.0, and === also tells 1 from 1.0; test for NaN with isnan.
• round sends ties to the even neighbour unless you pass another mode, rounding to decimal digits can be inexact because of the binary grid, and mod(5.5, -Inf) gives 5.5 from Julia 1.13 on.
• Subtracting nearly equal numbers loses digits, so rewrite such formulas, and the order of additions changes a sum: sum uses pairwise summation, and Kahan’s compensated sum was the most accurate in the run.