[R] Inverse of FAQ 7.31.

[Rolf Turner]

Why does R think these numbers ***are*** equal?

In a somewhat bizarre set of circumstances I calculated

     x0 <- 0.03580067
     x1 <- 0.03474075
     y0 <- 0.4918823
     y1 <- 0.4474461
     dx <- x1 - x0
     dy <- y1 - y0
     xx <- (x0 + x1)/2
     yy <- (y0 + y1)/2
     chk <- yy*dx - xx*dy + x0*dy - y0*dx

If you think about it ***very*** carefully ( :-) ) you'll see that 
``chk'' ought to be zero.

Blow me down, R gets 0.  Exactly.  To as many significant digits/decimal 
places
as I can get it to print out.

But .... I wrote a wee function in C to do the *same* calculation and 
dyn.load()-ed
it and called it with .C().  And I got -1.248844e-19.

This is of course zero, to all floating point arithmetic intents and 
purposes.  But if
I name the result returned by my call to .C() ``xxx'' and ask

     xxx >= 0

I get FALSE whereas ``chk >= 0'' returns TRUE (as does ``chk <= 0'', of 
course).
(And inside my C function, the comparison ``xxx >= 0'' yields ``false'' 
as well.)

I was vaguely thinking that raw R arithmetic would be equivalent to C 
arithmetic.
(Isn't R written in C?)

Can someone explain to me how it is that R (magically) gets it exactly 
right, whereas
a call to .C() gives the sort of ``approximately right'' answer that one 
might usually
expect?  I know that R Core is ***good*** but even they can't make C do 
infinite
precision arithmetic. :-)

This is really just idle curiosity --- I realize that this phenomenon is 
one that I'll simply have
to live with.  But if I can get some deeper insight as to why it occurs, 
well, that would
be nice.

     cheers,

         Rolf Turner