[R] Odp: precision issue?

[(Ted Harding)]

On 04-Mar-10 13:35:42, Duncan Murdoch wrote:
> On 04/03/2010 7:35 AM, (Ted Harding) wrote:
>> On 04-Mar-10 10:50:56, Petr PIKAL wrote:
>> > Hi
>> > 
>> > r-help-bounces at r-project.org napsal dne 04.03.2010 10:36:43:
>> >> Hi R Gurus,
>> >> 
>> >> I am trying to figure out what is going on here.
>> >> 
>> >> > a <- 68.08
>> >> > b <- a-1.55
>> >> > a-b
>> >> [1] 1.55
>> >> > a-b == 1.55
>> >> [1] FALSE
>> >> > round(a-b,2) == 1.55
>> >> [1] TRUE
>> >> > round(a-b,15) == 1.55
>> >> [1] FALSE
>> >> 
>> >> Why should (a - b) == 1.55 fail when in fact b has been defined
>> >> to be a - 1.55?  Is this a precision issue? How do i correct this?
>> > 
>> > In real world those definitions of b are the same but not in
>> > computer 
>> > world. See FAQ 7.31
>> > 
>> > Use either rounding or all.equal.
>> > 
>> >> all.equal(a-b, 1.55)
>> > [1] TRUE
>> > 
>> > To all, this is quite common question and it is documented in FAQs. 
>> > programs, therefore maybe a confusion from novices. 
>> > 
>> > I wonder if there could be some type of global option which will
>> > get rid of these users mistakes or misunderstandings by setting
>> > some threshold option for equality testing by use "==".
>> > 
>> > Regards
>> > Petr
>> > 
>> >> Alex
>>
>> Interesting suggestion, but in my view it would probably give
>> rise to more problems than it would avoid!
>>
>> The fundamental issue is that many inexperienced users are not
>> aware that once 68.08 has got inside the computer (as used by
>> R and other programs which do fixed-length binary arithmetic)
>> it is no longer 68.08 (though 1.55 is still 1.55).
>>   
> 
> I think most of what you write above and below is true, but your 
> parenthetical remark is not.  1.55 can't be represented exactly in the 
> double precision floating point format used in R.  It doesn't have a 
> terminating binary expansion, so it will be rounded to a binary 
> fraction.  I believe you can see the binary expansion using this code:
> 
> x <- 1.55
> for (i in 1:70) {
>   cat(floor(x)) 
>   if (i == 1) cat(".")
>    x <- 2*(x - floor(x))
> }
> 
> which gives
> 
> 1.100011001100110011001100110011001100110011001100110100000000000000000
> 
> Notice how it becomes a repeating expansion, with 0011 repeated 12 
> times, but then it finishes with 0100000000000000000, because we've run
> out of bits.
> So 1.55 is actually stored as a number which is ever so slightly
> bigger.
> 
> Duncan Murdoch

Of course you are quite right! My parenthetical remark was
an oversight and a blunder -- not so much a typo as a braino,
due to crossing wires between dividing by 2 (1.5 = 1 + 1/2)
and dividing by 10 ( 1.55 != 1 + 1/2 + (1/2)/2) )!

However, your detailed explanation will certainly be useful
to some. And your method of producing the binary expansion
is very neat indeed.

>> [the rest of my stuff snipped]

Ted.

--------------------------------------------------------------------
E-Mail: (Ted Harding) <Ted.Harding at manchester.ac.uk>
Fax-to-email: +44 (0)870 094 0861
Date: 04-Mar-10                                       Time: 20:11:39
------------------------------ XFMail ------------------------------