3

Can someone please help me understand why -1 + 1 <> 0 ?

Can someone please help me understand why I get three different values between the built-in function consum(), my function ct(), and Excel when they are all doing the same thing?

Now, I'm pretty sure the answer is a 'round' issue, but I can't figure out where that part of this issue is coming from. I mean, this all 'seems pretty straight forward.

In R, when I build the sequence 'a' and then run cumsum(a) I don't get the result of 0 like I expect to get. I also get a different answer if I try to calculate the same value using a function. Finally, I get a third answer if I try to calculate the same value using Excel.

This is what I get using cumsum():

> a<- seq(-1, 1, by=.1)
> a
 [1] -1.0 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1  0.0  0.1  0.2  0.3
[15]  0.4  0.5  0.6  0.7  0.8  0.9  1.0
> cumsum(a)
 [1] -1.000000e+00 -1.900000e+00 -2.700000e+00 -3.400000e+00 -4.000000e+00
 [6] -4.500000e+00 -4.900000e+00 -5.200000e+00 -5.400000e+00 -5.500000e+00
[11] -5.500000e+00 -5.400000e+00 -5.200000e+00 -4.900000e+00 -4.500000e+00
[16] -4.000000e+00 -3.400000e+00 -2.700000e+00 -1.900000e+00 -1.000000e+00
[21]  1.110223e-15

I wrote a quick function to test this and expected to get the same answer (or 0), but I get a totally different answer. Here is my function with its results:

ct<- function(x){
        result = 0
        for(i in 1:length(x)){
           cat(i, ": Result = ", result, " + ", x[i], " = ", result + x[i], "\n")
           result = result + x[i]
        }
}

> ct(a)
1 : Result =  0  +  -1  =  -1 
2 : Result =  -1  +  -0.9  =  -1.9 
3 : Result =  -1.9  +  -0.8  =  -2.7 
4 : Result =  -2.7  +  -0.7  =  -3.4 
5 : Result =  -3.4  +  -0.6  =  -4 
6 : Result =  -4  +  -0.5  =  -4.5 
7 : Result =  -4.5  +  -0.4  =  -4.9 
8 : Result =  -4.9  +  -0.3  =  -5.2 
9 : Result =  -5.2  +  -0.2  =  -5.4 
10 : Result =  -5.4  +  -0.1  =  -5.5 
11 : Result =  -5.5  +  0  =  -5.5 
12 : Result =  -5.5  +  0.1  =  -5.4 
13 : Result =  -5.4  +  0.2  =  -5.2 
14 : Result =  -5.2  +  0.3  =  -4.9 
15 : Result =  -4.9  +  0.4  =  -4.5 
16 : Result =  -4.5  +  0.5  =  -4 
17 : Result =  -4  +  0.6  =  -3.4 
18 : Result =  -3.4  +  0.7  =  -2.7 
19 : Result =  -2.7  +  0.8  =  -1.9 
20 : Result =  -1.9  +  0.9  =  -1 
21 : Result =  -1  +  1  =  4.440892e-16

If I change the last line in the for loop to this, then I get the expected answer of 0:

result = round(result + x[I], digits = 2)

In Excel, using the same logic as in my ct() function, I get a final result of -2.886580E-15 (without rounding the values).

GEOCHET
  • 21,119
  • 15
  • 74
  • 98
dave
  • 147
  • 1
  • 1
  • 12

2 Answers2

5

This is the nature of using a fixed-precision representation with values that it cannot represent exactly.

Just like 1/3 can't be represented exactly with a fixed number of decimal places, 0.1 can't be represented exactly with a fixed number of binary places. So just like 3 x (1/3) can't possibly give you 1 with a fixed number of decimal places, adding multiples of 0.1 will never give you exactly 1 in fixed-precision binary.

So, let's look at six-precision decimal representation just to see this more clearly (this is used to indicate values, as opposed to representations):
1 -> 1.000000
1/3 -> .333333
2/3 -> .666667
3 -> 3.000000

This gives:

1/3 + 2/3 -> 0.333333 + 0.666667 -> 1.000000 -> 1 (yay)

1/3 + 1/3 -> 0.333333 + 0.333333 -> 0.666666 (not 2/3, oh well)

3 * 1/3 -> 3.00000 * 0.333333 -> .999999 (not 1, oh well)

How you handle this is up to you, but this should be expected behavior.

To address your last question, why doing the "same thing" two different ways can produce different results, it comes from intermediate rounding. If you've ever done a computation with a calculator, writing down some of the partial intermediate results, you know it can make a difference which intermediate results you write down.

David Schwartz
  • 179,497
  • 17
  • 214
  • 278
  • 3
    Also useful, *R* FAQ #7.31, [*"why doesn't R think these numbers are equal"*](http://cran.r-project.org/doc/FAQ/R-FAQ.html#Why-doesn_0027t-R-think-these-numbers-are-equal_003f). – r2evans May 21 '15 at 19:49
  • 2
    I'll +1 anything that cites [What Every Computer Scientist Should Know About Floating-Point Arithmetic](http://www.validlab.com/goldberg/paper.pdf), even indirectly. – David Schwartz May 21 '15 at 19:53
  • 2
    and.. [Why are these numbers not equal?](http://stackoverflow.com/questions/9508518/why-are-these-numbers-not-equal) – Pierre L May 21 '15 at 19:53
  • @DavidSchwartz, perhaps I'm placating the impatience of many readers. Although truly placating would just paste the text itself instead of a link ... – r2evans May 21 '15 at 19:57
  • Okay, I've read through the links and I "started" plowing though the ACM PDF. These basically give me the answer for why -1 +1 <> 0. I would appreciate some insight as to why my function, `ct(a)`, doesn't give the same answer as `cumsum(a)` since it does the same thing as `cumsum(a)`. (Doesn't it?) – dave May 22 '15 at 14:03
  • @dave See my update to the answer. It happens when which intermediate results you round affects the final result. For example, if the code says `(A+B)+(C+D)`, it may be that `(A+B)` gets rounded, `(C+D)` gets rounded, both get rounded, or neither gets rounded, depending on the availability of floating point registers and what combination of those variables happen to already be in registers when the code begins. – David Schwartz May 22 '15 at 20:11
0

I'd guess it's just rounding issues. If you use the seq.int function to make a vector from -10 to 10 and then do cumsum you get a sum of 0:

> seq.int(-10,10,1)
[1] -10  -9  -8  -7  -6  -5  -4  -3  -2  -1   0   1   2   3   4   5   6   7   8   9  10
> cumsum(seq.int(-10,10,1))
[1] -10 -19 -27 -34 -40 -45 -49 -52 -54 -55 -55 -54 -52 -49 -45 -40 -34 -27 -19 -10   0

If you really want to do a sequence between -1 and 1 then just divide the integer sequence by 10L. 

cumsum(seq.int(-10,10,1)/10L)
[1] -1.0 -1.9 -2.7 -3.4 -4.0 -4.5 -4.9 -5.2 -5.4 -5.5 -5.5 -5.4 -5.2 -4.9 -4.5 -4.0 -3.4 -2.7
[19] -1.9 -1.0  0.0

You are still going to deal with some rounding errors as always but this appears to be below R's threshold of rounding to 0.

Matt Mills
  • 588
  • 1
  • 6
  • 14