16

Can somebody please explain me the following output. I know that it has something to do with floating point precision, but the order of magnitue (difference 1e308) surprises me.

0: high precision

> 1e-324==0
[1] TRUE
> 1e-323==0
[1] FALSE

1: very unprecise

> 1 - 1e-16 == 1
[1] FALSE
> 1 - 1e-17 == 1
[1] TRUE
smci
  • 32,567
  • 20
  • 113
  • 146
  • 4
    Yes, this is the [R FAQ 7.31](http://cran.r-project.org/doc/FAQ/R-FAQ.html#Why-doesn_0027t-R-think-these-numbers-are-equal_003f) –  Jul 20 '14 at 06:14
  • 2
    In particular, try `options(digits = 22); x <- 1 - 1e-16; y <- 1 - 1e-17 == 1; print(y); print(x)` – Robert Krzyzanowski Jul 20 '14 at 06:16
  • @DavidArenburg I don't think that the “duplicate” you offer is the same question at all. The OP does not seem to say they are surprised about some value having to be represented identically but about the density of representable numbers near 0 and near 1. – Pascal Cuoq Jul 20 '14 at 07:31
  • @PascalCuoq, maybe you are right and I was too harsh. i'll leave it to Roman to decide – David Arenburg Jul 20 '14 at 07:36
  • @Pascal this question is about the density of representable numbers near 0 and near 1, not about the fact that some simple decimal numbers are not represented exactly, and it was closed with the wrong duplicate. – Pascal Cuoq Jul 20 '14 at 07:38
  • 2
    @DavidArenburg you mean Roman, the user with a gold badge who just closed the question without reading the comments or looking for the real duplicate if there indeed is one? The question is in good hands. – Pascal Cuoq Jul 20 '14 at 07:40
  • @PascalCuoq, chill, he just reopened it. It was misunderstanding caused mainly by me – David Arenburg Jul 20 '14 at 07:42
  • @Pascal the question asked here is not the one that is answered at http://cran.r-project.org/doc/FAQ/R-FAQ.html#Why-doesn_0027t-R-think-these-numbers-are-equal_003f – Pascal Cuoq Jul 20 '14 at 08:59
  • @PascalCuoq When I answered, there was no duplication flagged. It is what I meant. –  Jul 20 '14 at 09:10
  • 1
    @user3710546: In fairness the numerical examples illustrated in [R FAQ 7.31](https://cran.r-project.org/doc/FAQ/R-FAQ.html#Why-doesn_0027t-R-think-these-numbers-are-equal_003f) are seriously awful; there's no addition example, so it fails to illustrate any of `.Machine$ulp, eps, min.exp, xmin`. `(1 - 1e-17 == 1) == TRUE` is a much clearer and simpler example. – smci Feb 02 '17 at 08:07

2 Answers2

16

R uses IEEE 754 double-precision floating-point numbers.

Floating-point numbers are more dense near zero. This is a result of their being designed to compute accurately (the equivalent of about 16 significant decimal digits, as you have noticed) over a very wide range.

Perhaps you expected a fixed-point system with uniform absolute precision. In practice fixed-point is either wasteful or the ranges of each intermediate computation must be carefully estimated beforehand, with drastic consequences if they are wrong.

Positive floating-point numbers look like this, schematically:

+-+-+-+--+--+--+----+----+----+--------+--------+--------+--
0

The smallest positive normal double-precision number is 2 to the power of the minimal exponent. Near one, the double-precision floating-point numbers are already spread quite wide apart. There is a distance of 2-53 from one to the number below it, and a distance of 2-52 from one to the number above it.

Pascal Cuoq
  • 79,187
  • 7
  • 161
  • 281
  • Nearly 17 significant digits - `1 - 1e-17` would be different from `1` if it supported full 17 digit precision. – Patricia Shanahan Jul 20 '14 at 08:37
  • @PatriciaShanahan I wanted to convey that a number of decimal digit was only a rough measure of the precision of binary floating-point anyway, and ended up using “equivalent” which wasn't strong enough. I have added “about”. – Pascal Cuoq Jul 20 '14 at 09:02
  • 2
    While we're nitpicking, 'nearly 16 digits' would be a more accurate description than 'nearly 17 digits'. (2^53 is roughly 10^15.95). – Mark Dickinson Jul 20 '14 at 09:10
9

As per @PascalCuoq's answer, because the FP precision of an R double on an IEEE 754-compliant platform (e.g. an x86) is not quite 16 decimal digits:

# Machine ULP in decimal digits...
.Machine$double.ulp.digits * log10(2)
-15.65...

# Note the direct relationship between ULP digits and EPS:
.Machine$double.ulp.digits = -52 
2.22 e-16 = .Machine$double.eps == 2^.Machine$double.ulp.digits

Hence 1 - 1e-16 is already very close to ULP, and 1 - 1e-17 is beyond ULP, and gets rounded to FP 1.0

See the R documentation for .Machine : "Numerical Characteristics of the Machine" . In particular look at the difference between EPS and ULP.

(ULP is defined wrt the FP number 1. The bigger your FP number gets, the bigger the value of the last bit, and the more crude a rounding operation is)

As to where the quantity 1e-323 comes from: you're confusing ULP with the minimum-representable FP value, which is far smaller.

The minimum-representable normalized positive FP value has exponent e-308, as per IEEE 754 Double-precision examples ...

# Minimum-representable normalized positive FP value is...
.Machine$double.xmin
2.225..e-308

# ...which would correspond to this base-2 exponent...
log10(.Machine$double.xmin) / log10(2)
-1022
# ...or this base-10 exponent...
.Machine$double.min.exp * log10(2)
-307.65...

But we can get slightly smaller, if we use an UNnormalized FP number, i.e. all the leading mantissa bits are 0. Hence as you empirically discovered, the minimum-representable UNnormalized positive FP value was somewhere between 1e-324 and 1e-323. That's because we have 52 mantissa bits, hence the numerical value of the LSB is 2^51 or 10^15.35 smaller:

# Exponent of Minimum-representable UNnormalized positive FP value 
log10(.Machine$double.xmin) - (.Machine$double.digits * log10(2))
-323.607...

(Why can we not empirically discover it exactly? Because IEEE-754 internally carries a few guard digits before rounding)

(Note also the parameter that says the representation is base-2: .Machine$double.base = 2 )

smci
  • 32,567
  • 20
  • 113
  • 146
  • 4
    To the downvoters and everyone: please actually comment on what if anything is wrong. I put some effort into this one. – smci Feb 04 '17 at 04:18