Explain 1 - 3*(4/3 - 1) = 2.2204e-16 in Matlab

Question

I am trying to understand the double precision numbers in Matlab. Why is this 1 - 3*(4/3 - 1) not equal to zero?

Floating point does not represent decimal numbers exactly. See [What every computer scientist should know about floating-point arithmetic](https://ece.uwaterloo.ca/~dwharder/NumericalAnalysis/02Numerics/Double/paper.pdf). — Robert Harvey, Mar 12 '13 at 03:24
My unsubstantiated guess would be that due to the underlying representation of doubles, it gets calculated as 2.22044604925031E-16. — lzcd, Mar 12 '13 at 03:25
Related question: [Is this a Matlab bug? Do you have the same issue?](http://stackoverflow.com/questions/13699596/is-this-a-matlab-bug-do-you-have-the-same-issue/13699708#13699708) — Eitan T, Mar 12 '13 at 08:23
You may be interested in reading my article "The Answer is One (Unless You Use Floating-Point)": http://www.exploringbinary.com/the-answer-is-one-unless-you-use-floating-point/ — Rick Regan, Mar 12 '13 at 12:33

score 6 · Accepted Answer · answered Mar 12 '13 at 12:59

The real number 4/3 is not representable in double precision (or any other binary floating-point format) because it is not a dyadic rational. Thus, when you compute 4/3 in MATLAB, the value you get is rounded to the closest representable double-precision number, which is exactly:

1.3333333333333332593184650249895639717578887939453125

Subtracting 1 from this value is exact (it is a well-known theorem of FP error-analysis that subtracting numbers within a factor of two of each other is exact), so the result of 4/3 - 1 is:

0.3333333333333332593184650249895639717578887939453125

It happens that the result of multiplying this number by three is also exactly representable:

0.9999999999999997779553950749686919152736663818359375

Finally, subtracting from 1.0 is also exact (by the theorem I referenced earlier):

0.0000000000000002220446049250313080847263336181640625

So there is only a single source of rounding error in your computation, owing to the fact that 4/3 cannot be represented as a double, and the final result of the computation is simply that initial error carried forward.

Austin Mullins · Answer 2 · 2013-03-12T03:50:46.167

1

Running this on a scientific calculator - log₂(2.2204e-16) - yields (almost) exactly -52. In other words, Matlab stores 52 bits of precision in a double, with another 5 bits for the exponent (4 + sign) and one for the sign of the significand. This is in line with the IEEE 754 implementation: 53 bits for the significand (52 + sign) and 5 for the exponent. All is well! As always, you should test whether two floating point numbers are close enough, not whether they are exactly equal. An appropriate example in Matlab would look like:

if abs(x - y) < 1e-15
      % Some code to execute when x and y are approximately equal
else
      % Some other code
end

The statement of adherence to IEEE 754 comes from the wikipedia article.

edited Mar 12 '13 at 03:50

answered Mar 12 '13 at 03:36

Austin Mullins

7,307
2
33
48

Whoa, a mod saw my answer. That's a first. – Austin Mullins Mar 12 '13 at 03:49
Also, I edited the answer above by including the "abs" function. Very important! – Austin Mullins Mar 12 '13 at 03:51

score 0 · Answer 3 · answered Mar 12 '13 at 12:53

The trouble starts with the calculation of 4/3. The exact answer is representable in neither a finite number of decimal digits nor a finite number of bits. The result will be stored as 4/3+r where r is a small absolute value signed number representing the difference between the real value of 4/3 and the closest IEEE 754 64-bit binary float to 4/3, the rounding error.

Subtracting 1 results in 1/3+r. Multiplying by 3 gets 1+3r. Subtracting it from 1 gets -3r.

The end result is -3 times the rounding error in the original representation of 4/3.

Explain 1 - 3*(4/3 - 1) = 2.2204e-16 in Matlab

3 Answers3