[R] Inverse of FAQ 7.31.

[peter dalgaard]

On Aug 2, 2011, at 08:02 , Rolf Turner wrote:

> 
> 
> 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.

I think the long and the short of it is that R lost a couple of bits of precision that C retained. This sort of thing happens if R stores things into 64 bit floating point objects while C keeps them in 80 bit CPU registers. In general, floating point calculations do not obey the laws of math, for example the associative law (i.e., (a+b)-c ?= a+(b-c), especially if b and c are large and nearly equal), so any reordering of expressions by the compiler may give a slightly different result.

-- 
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Email: pd.mes at cbs.dk  Priv: PDalgd at gmail.com
"Døden skal tape!" --- Nordahl Grieg