[Rd] Non identical numerical results from R code vs C/C++ code?

[Duncan Murdoch]

On 10/09/2010 6:46 AM, Renaud Gaujoux wrote:
> Hi,
> 
> suppose you have two versions of the same algorithm: one in pure R, the 
> other one in C/C++ called via .Call().
> Assuming there is no bug in the implementations (i.e. they both do the 
> same thing), is there any well known reason why the C/C++ implementation 
> could return numerical results non identical to the one obtained from 
> the pure R code? (e.g. could it be rounding errors? please explain.)
> Has anybody had a similar experience?

R often uses extended reals (80 bit floating point values on Intel 
chips) for intermediate values.  C compilers may or may not do that.
> 
> By not identical, I mean very small differences (< 2.4 e-14), but enough 
> to have identical() returning FALSE. Maybe I should not bother, but I 
> want to be sure where the differences come from, at least by mere curiosity.
> 
> Briefly the R code perform multiple matrix product; the C code is an 
> optimization of those specific products via custom for loops, where 
> entries are not computed in the same order, etc... which improves both 
> memory usage and speed. The result is theoretically the same.

Changing the order of operations will often affect rounding.  For 
example, suppose epsilon is the smallest number such that 1 + epsilon is 
not equal to 1.  Then 1 + (epsilon/2) + (epsilon/2) will evaluate to 
either 1 or 1 + epsilon, depending on the order of computing the additions.

Duncan Murdoch

> 
> Thank you,
> Renaud
>