3

I encounted the following precision issue. I define c = a + b, but the result for 'c - a - b' does not equal zero. I am not sure how to resolve this.

'''

In[1]:= 
a = N[Exp[3]];
b = N[Exp[4]];
c = a + b;
c - a - b

Out[4]= 7.10543*10^-15

'''

HB_LUKE
  • 33
  • 3
  • 1
    What happens when you do Chop[%] on the final output? does it become zero? see help on Chop or you can use infinite precision. a = N[Exp[3], Infinity]; b = N[Exp[4], Infinity]; and you will get exact zero without doing Chop – Nasser Sep 27 '21 at 07:47
  • 1
    But note that N[Exp[3], Infinity] just keeps it as Exp[3] and will not make it real number. Once you convert things to approximate numbers, you'll get approximate numbers as result of computation. – Nasser Sep 27 '21 at 07:56
  • 1
    notice also that N[...] use machine precision (unless you give second argument to N). on same PC, I typed the same thing on Matlab, and got same output: a=exp(3); b=exp(4); c=a+b; c-a-b gives 7.1054e-15 This is the nature of using real numbers. Mathematica graphics – Nasser Sep 27 '21 at 08:05
  • 2
    You can avoid some of the problems by using c - (a + b) instead of c - a - b. Now you will get 0. as output. This is because the order makes difference in terms of real numbers computations (large number - small number, vs. difference of two numbers close to each others). see What Every Computer Scientist Should Know About Floating-Point Arithmetic – Nasser Sep 27 '21 at 08:11
  • 1
    Everything is indeed correct. You can see the same phenomenon even on a scientific calculator. As @Nasser said, this is the nature of floating point computations. – yarchik Sep 27 '21 at 09:51
  • Thank you, Nasser and Yarchik! Your comments are very helpful. – HB_LUKE Sep 27 '21 at 20:16

1 Answers1

1

This phenomenon occurs for simple cases too:

0.3 + 1.1 - 1.4

2.22045*10^-16`

These numbers do not have finite length in binary and so they're truncated when expressed as an inexact number. We can see this with RealDigits:

RealDigits[1.1, 2]

{{1, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0,
  1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 
  0, 0, 1, 1, 0, 0, 1, 1, 0, 1, 0}, 1}

The identity $1.1_{10} = 1.0001100110011_2...$ gets truncated to $53$ digits to the right of the decimal (on a $64$ bit machine). The same happens for $0.3$ and $1.4$ — and in fact only dyadic rationals can be represented exactly in floating point arithmetic.

We can see the error by expressing the truncated reals as rational numbers:

SetPrecision[0.3 + 1.1, ∞] - SetPrecision[1.4, ∞]

1/4503599627370496

N[%]

2.22045*10^-16
Greg Hurst
  • 35,921
  • 1
  • 90
  • 136