Numbers, floating point and rounding

Type 0.1 + 0.2 == 0.3 into R and the answer is FALSE, although 0.1 + 0.2 prints as 0.3. Nothing is broken: R stores numbers as binary fractions, most decimal numbers have no exact binary form, and the same fact explains why round(2.675, 2) is 2.67 and why three shares of 1,000.00, each rounded to 333.33, add up to 999.99.

On runs of R 4.6.1 you will see what R stores for a number and compare doubles safely, predict how round() treats halves, recognise where integers overflow and doubles stop counting by ones, and print numbers in any format without changing the values that calculations use.

Doubles are binary fractions, which is why 0.1 plus 0.2 is not 0.3

A number such as 0.1 is stored as a doubledouble: R's type for real numbers: a binary fraction with 53 bits of precision, about 15 significant decimal digits, in the IEC 60559 (IEEE 754) format. Most decimals, such as 0.1, are stored as the nearest such fraction.: a binary fraction with 53 bits of precision, in the IEC 60559 format that the help page of double says every R platform uses. A fraction whose denominator, in lowest terms, is not a power of 2, such as 1/10, has no exact binary form, just as 1/3 has no exact decimal one, so R stores the double nearest to it. Step through the four numbers of the comparison:

The run prints the same four values to twenty decimal places with sprintf(), and then the comparisons:

from session import run, version

version()
# R reads binary.R as if each line were typed at the console.
run("""
R --vanilla --quiet -f binary.R
""")

# R version 4.6.1 -- "Happy Hop"
$ R --vanilla --quiet -f binary.R
> # The sum prints as 0.3, but it is not 0.3.
> 0.1 + 0.2
[1] 0.3
> 0.1 + 0.2 == 0.3
[1] FALSE
> # Twenty decimal places show the stored values.
> sprintf("%.20f", c(0.1, 0.2, 0.1 + 0.2, 0.3))
[1] "0.10000000000000000555" "0.20000000000000001110" "0.30000000000000004441"
[4] "0.29999999999999998890"
> (0.1 + 0.2) - 0.3
[1] 5.551115e-17
> # They are neighbouring doubles, 2^-54 apart.
> (0.1 + 0.2) - 0.3 == 2^-54
[1] TRUE
> # The same happens inside a sequence.
> seq(0, 1, by = 0.1)[4] == 0.3
[1] FALSE
> # Compare with a tolerance instead.
> all.equal(0.1 + 0.2, 0.3)
[1] TRUE
> isTRUE(all.equal(0.1 + 0.2, 0.3))
[1] TRUE
> abs((0.1 + 0.2) - 0.3) < 1e-9
[1] TRUE
> # A real difference is reported, not hidden.
> all.equal(1, 1.1)
[1] "Mean relative difference: 0.1"
> isTRUE(all.equal(1, 1.1))
[1] FALSE

The sum lands on the double just above 0.3 and the literal 0.3 on the one just below, 2^-54 apart, which the run confirms. The fourth value of seq(0, 1, by = 0.1) carries the same kind of error, so comparing it with 0.3 fails too.

The fix is to compare with a tolerance. all.equal(x, y) accepts differences smaller than about 1.5e-8, measured relative to the size of the numbers, and returns TRUE or a sentence that describes the difference, so inside if () write isTRUE(all.equal(x, y)). When you know the scale of your data, abs(x - y) < tol with your own tol says the same thing explicitly.

round() goes to the even digit, but only for halves it can see

round(c(0.5, 1.5, 2.5)) returns 0 2 2. For a half, the help page of round() says R expects the IEC 60559 rule: go to the even digit, so 0.5 rounds down to 0 and 1.5 up to 2. Halves then round down as often as up, so a long column of rounded values does not drift upwards. But a decimal “half” such as 2.675 is usually not stored as a half at all:

The run puts true halves beside decimal ones, with the stored values printed to twenty places:

from session import run, version

version()
# R reads rounding.R as if each line were typed at the console.
run("""
R --vanilla --quiet -f rounding.R
""")

# R version 4.6.1 -- "Happy Hop"
$ R --vanilla --quiet -f rounding.R
> # A half that is stored exactly goes to the even digit.
> round(c(0.5, 1.5, 2.5, -0.5, -1.5))
[1]  0  2  2  0 -2
> round(0.125, 2)
[1] 0.12
> # Most decimal "halves" are stored a little off.
> round(2.675, 2)
[1] 2.67
> sprintf("%.20f", 2.675)
[1] "2.67499999999999982236"
> round(0.15, 1)
[1] 0.1
> sprintf("%.20f", 0.15)
[1] "0.14999999999999999445"
> # The other ways to round.
> ceiling(2.1)
[1] 3
> floor(-2.1)
[1] -3
> trunc(-2.7)
[1] -2
> signif(123456, 2)
[1] 120000
> round(123456, -3)
[1] 123000

0.125 is 1/8, a binary fraction, so it is stored exactly, it is a true half, and round(0.125, 2) goes to the even digit, 0.12. 2.675 and 0.15 are stored a little below their printed values, so they round down, to 2.67 and 0.1. The help page warns about exactly this: the rounding rule applies to the represented number, not the printed one.

ceiling() rounds up, floor() down, trunc() towards zero, signif() to a number of significant digits, and a negative digits rounds to tens, hundreds or thousands: round(123456, -3) is 123000.

Integers overflow to NA, and doubles stop counting by ones at 2^53

R’s integers are 32-bit, so .Machine$integer.max is 2,147,483,647, and one more is out of range. The help page of integer says that as.integer() turns a value beyond the range into NA, and integer arithmetic does the same; both warn, as the run shows. Doubles have the opposite trade-off: a far larger range, about 2e+308, with 53 bits of precision, so whole numbers are exact only up to 2^53, 9,007,199,254,740,992.

from session import run, version

version()
# R reads limits.R as if each line were typed at the console.
run("""
R --vanilla --quiet -f limits.R
""")

# R version 4.6.1 -- "Happy Hop"
$ R --vanilla --quiet -f limits.R
> # Integers are 32-bit: this is the largest one.
> .Machine$integer.max
[1] 2147483647
> .Machine$integer.max + 1L
[1] NA
Warning message:
In .Machine$integer.max + 1L : NAs produced by integer overflow
> as.integer(2^31)
[1] NA
Warning message:
NAs introduced by coercion to integer range 
> # The same sum in doubles is no problem.
> .Machine$integer.max + 1
[1] 2147483648
> # Doubles hold whole numbers exactly up to 2^53.
> format(2^53, scientific = FALSE)
[1] "9007199254740992"
> 2^53 == 2^53 + 1
[1] TRUE
> # The gap between neighbouring doubles grows with size.
> 2^(floor(log2(c(1, 1e6, 2^53))) - 52)
[1] 2.220446e-16 1.164153e-10 2.000000e+00
> .Machine$double.eps
[1] 2.220446e-16

The second-to-last line shows why. The gap between one double and the next grows with the size of the number: about 2.2e-16 near 1, about 1.2e-10 near a million, and 2 at 2^53, so 2^53 + 1 has nowhere to go and rounds back to 2^53. .Machine$double.eps, 2.220446e-16, is that gap at 1: its help page defines it as the smallest positive number x for which 1 + x != 1.

The last digits can differ between computers, so never test them

sum(rep(0.1, 10)) == 1 is FALSE in the run below, on an Apple-silicon Mac. The help page of sum() says it uses extended-precision accumulators where possible, which is platform-dependent, and capabilities("long.double") says whether this build of R has a long double wider than a double. Writing R Extensions adds that current ARM CPUs have neither extended precision nor long doubles:

from session import run, version

version()
# R reads platform.R as if each line were typed at the console.
run("""
R --vanilla --quiet -f platform.R
""")

# R version 4.6.1 -- "Happy Hop"
$ R --vanilla --quiet -f platform.R
> # Does this build of R have a long double wider than a double?
> capabilities("long.double")
long.double 
      FALSE 
> # sum() adds in long double where there is one; here there is not.
> sum(rep(0.1, 10)) == 1
[1] FALSE
> print(sum(rep(0.1, 10)), digits = 17)
[1] 0.99999999999999989
> # A loop adds in plain doubles on every computer.
> total <- 0
> for (v in rep(0.1, 10)) total <- total + v
> print(total, digits = 17)
[1] 0.99999999999999989
> # A check that gives the same answer on every computer.
> isTRUE(all.equal(sum(rep(0.1, 10)), 1))
[1] TRUE
> # Exact at special values on every platform since R 4.6.0.
> exp(0) == 1
[1] TRUE
> log1p(0) == 0
[1] TRUE

On this Apple-silicon build there is no long double, so sum() adds in plain doubles, exactly as the loop does, and both reach 0.99999999999999989. On a build where capabilities("long.double") is TRUE, sum() adds the same ten values in a wider format, and the same comparison can come out TRUE. Code whose result depends on that last bit behaves differently from one computer to the next; code that compares with all.equal() does not.

format(), formatC() and sprintf() change the print, not the value

x <- 1234567.891 prints as 1234568, because R shows 7 significant digits by default. Printing is a view, never a change: format(), formatC() and sprintf() return text for display, and options(digits = 4) changes how numbers print, while x itself keeps every digit it had.

from session import run, version

version()
# R reads formatting.R as if each line were typed at the console.
run("""
R --vanilla --quiet -f formatting.R
""")

# R version 4.6.1 -- "Happy Hop"
$ R --vanilla --quiet -f formatting.R
> x <- 1234567.891
> # Printing shows 7 significant digits by default.
> x
[1] 1234568
> # format(), formatC() and sprintf() return text.
> format(x, nsmall = 2)
[1] "1234567.89"
> formatC(x, format = "f", digits = 2, big.mark = ",")
[1] "1,234,567.89"
> sprintf("%.2f", x)
[1] "1234567.89"
> # options(digits =) changes printing, nothing else.
> options(digits = 4)
> 1 / 3
[1] 0.3333
> options(digits = 7)
> # The stored value never changed.
> x == 1234567.891
[1] TRUE
> # An amount of 1,000.00 in its smallest unit, split three ways.
> total <- 100000L
> share <- rep(total %/% 3L, 3)
> extra <- seq_len(total %% 3L)
> share[extra] <- share[extra] + 1L
> share
[1] 33334 33333 33333
> sum(share) == total
[1] TRUE
> sprintf("%.2f", share / 100)
[1] "333.34" "333.33" "333.33"
> # Rounding each share instead loses one unit.
> naive <- rep(round(1000 / 3, 2), 3)
> naive
[1] 333.33 333.33 333.33
> sum(naive)
[1] 999.99

• format(x, nsmall = 2) keeps at least 2 decimals, "1234567.89".
• formatC(x, format = "f", digits = 2, big.mark = ",") adds thousands separators, "1,234,567.89".
• sprintf("%.2f", x) uses C’s formats, the same as in many other languages, "1234567.89".
• options(digits = 4) changes every later print, here 1 / 3 as 0.3333, until you set it back.

The second half of the run shows how to keep money exact. An amount of 1,000.00 held as 100000 of its smallest unit, an integer, splits three ways with integer division, %/%, and the remainder, %%, goes to the first share: 333.34, 333.33 and 333.33, which add up exactly. Rounding each share to two decimals gives 333.33 three times and loses one unit. Keep amounts in whole minor units, compute with them, and format them only for display.

Practise: predict, then check

The last set of lines mixes binary fractions, halves, integer limits and printing. Commit to an answer for each.

from session import run, version

version()
# R reads practice.R as if each line were typed at the console.
run("""
R --vanilla --quiet -f practice.R
""")

# R version 4.6.1 -- "Happy Hop"
$ R --vanilla --quiet -f practice.R
> x <- 0.3 - 0.1
> x == 0.2
[1] FALSE
> isTRUE(all.equal(x, 0.2))
[1] TRUE
> round(c(1.5, 2.5, 3.5))
[1] 2 2 4
> round(1.005, 2)
[1] 1
> .Machine$integer.max + 1
[1] 2147483648
> format(2/3, digits = 3)
[1] "0.667"

0.3 - 0.1 misses 0.2 for the same reason 0.1 + 0.2 misses 0.3, and all.equal() accepts it; the halves went to the even digits 2, 2 and 4; 1.005 is stored just below 1.005, so it rounds to 1; adding a double 1 to the largest integer gives a double, with no overflow; and format() rounded the printed text, not the value.

Last, the ideas of this page meet in one adaptive round: the weakest idea is asked first, each answer comes with its reason, and a question you miss returns once more before the summary:

Interview questions

Warm-up: why is 0.1 + 0.2 == 0.3 FALSE in R, and how should you compare such numbers?

R keeps a number such as 0.1 as a double, a binary fraction with 53 bits of precision. None of 0.1, 0.2 and 0.3 can be written as such a fraction, so R keeps the closest double to each, and the three small errors do not cancel: 0.1 + 0.2 ends up as 0.30000000000000004441 and 0.3 as 0.29999999999999998890, neighbouring doubles about 5.55e-17 apart. == compares the stored values, so it says FALSE, while printing to 7 digits shows both as 0.3. Compare with a tolerance instead: isTRUE(all.equal(a, b)), which accepts relative differences below about 1.5e-8, or abs(a - b) < tol with a tol that suits the scale of the data. The same applies in every language that uses IEEE 754 doubles; it is a property of binary floating point, not of R.

Key takeaways

• A double is a binary fraction with 53 bits of precision, so most decimals are stored inexactly and 0.1 + 0.2 == 0.3 is FALSE; compare with isTRUE(all.equal(x, y)) or abs(x - y) < tol.
• round() sends a true half to the even digit, round(2.5) is 2, while a decimal like 2.675 rounds by its stored value, so round(2.675, 2) is 2.67.
• Integers stop at 2,147,483,647 and overflow to NA with a warning, while doubles count whole numbers exactly only up to 2^53.
• The last bits can differ between computers, as sum() with and without a long double shows, so tests must never depend on them.
• format(), formatC(), sprintf() and options(digits) change only how numbers print; keep money in whole minor units and format it at the end.