On 30-Sep-09 19:32:46, Peter Dalgaard wrote: > Martin Batholdy wrote: >> hum, >> >> can you explain that a little more detailed? >> Perhaps I miss the background knowledge - but it seems just absurd to >> me. >> >> 0.1+0.1+0.1 is 0.3 - there is no rounding involved, is there? >> >> why is >> x <- 0.1 + 0.1 +0.1 >> not equal to >> y <- 0.3 > > Remember that this is in BINARY arithmetic. It's really not any > stranger than the fact that 1/3 + 1/3 != 2/3 in finite accuracy > decimal arithmetic (0.33333 + 0.33333 = 0.66666 != 0.66667).
Let me perhaps try to spell it out. Just as in decimal 1/7 = 0.142857142857142857142857... with the sequence 142857 recurring for ever (i.e. 1/7 is not exactly representable in decimal), so 1/10 = 0.0001100110011001100110011... in binary, with the sequence 0011 recurring for ever, so 1/10 is not exactly representable in binary (though of course it is in decimal). When you write (above) "x <- 0.1 + 0.1 +0.1", you are writing in decimal. R does not do decimal arithmetic (it could be implemented in a special package perhaps, as it is in the Unix program 'bc'; but it isn't there), and instead does its arithmetic internally in binary. Therefore when you enter 0.1 + 0.1 + 0.1 R will convert each "0.1" into an internal binary representation, and add these up. On the other hand, when you enter "Y <- 0.3", R will convert 0.3 into an internal binary representation which is NOT identical to the result of adding up the three internal representations of "0.1". These internal representations are to a fixed number of binary places. Example in R: The binary fraction 0.0001100 is (1/2^4 + 1/2^5) (and these are fine since they are reciprocals of power of 2 and so exactly represented). The next chunk in the expansion of 1/10 is 0.00011001100 which is the above, plus 1/16 of itself. And 0.000110011001100 is the above plus 1/16 plus 1/16^2 of itself, and so on ... Code: M<-1; k<-1/16 for(i in (1:14)){ print((1/2^4 + 1/2^5)*M, 17) M<-M+k; k<-k/16 } # [1] 0.09375 # [1] 0.099609375 # [1] 0.0999755859375 # [1] 0.09999847412109375 # [1] 0.0999999046325684 # [1] 0.0999999940395355 # [1] 0.099999999627471 # [1] 0.0999999999767169 # [1] 0.0999999999985448 # [1] 0.09999999999990905 # [1] 0.0999999999999943 # [1] 0.09999999999999964 # [1] 0.1 # [1] 0.1 (with no further change to 17 decimal places). Now let's do "0.1 + 0.1 + 0.1" explicitly in binary arithmetic (to, say, 27 binary places throughout): 0.000110011001100110011001100 * see below + 0.000110011001100110011001100 * see below ------------------------------- 0.001100110011001100110011000 + 0.000110011001100110011001100 * see below ------------------------------- 0.010011001100110011001100100 =============================== which is incorrect, since the last 2 binary digits in the result should be "10", not "00". The reason is that the next 2 digits in the truncated representation of 1/10, at the three places marked "*" above, are "11", and these would get added in if they had been present. But they are not present, because of truncation, and so do not get added in. To 27 binary places, the correct representation of 3/10 (truncated to 27 binary places)is 0.010011001100110011001100110 Hence (omitting various further details ... ) 3/10 == (1/10+1/10+1/10) # [1] FALSE Because of the cumulative effects of truncation errors like the above, and rounding errors which can also be introduced along the way, the method (1:9)/10 or 0.1*(1:9) is better, since each term only acquires error from the specific arithmetic operation which gives rise to it, and these errors do not accumulate. PS: I'm aware of all the grandmothers on the list who do not need instruction in ovisuction, but this sort of question occurs so frequently -- and is not spelled out in detail in the docs -- that I thought it would be useful to actually exhibit an example of the working for the benefit of people who do wonder just what is going on then things do not seem to match up as expected! Ted. -------------------------------------------------------------------- E-Mail: (Ted Harding) <ted.hard...@manchester.ac.uk> Fax-to-email: +44 (0)870 094 0861 Date: 30-Sep-09 Time: 21:53:25 ------------------------------ XFMail ------------------------------ ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.