From murdoch.duncan at gmail.com  Fri Sep 10 13:00:25 2010
From: murdoch.duncan at gmail.com (Duncan Murdoch)
Date: Fri, 10 Sep 2010 07:00:25 -0400
Subject: [Rd] Non identical numerical results from R code vs C/C++ code?
In-Reply-To: <4C8A0C86.8030307@cbio.uct.ac.za>
References: <4C8A0C86.8030307@cbio.uct.ac.za>
Message-ID: <4C8A0FC9.9060400@gmail.com>

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
>



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