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

[Duncan Murdoch]

On 10/09/2010 7:07 AM, Renaud Gaujoux wrote:
> Thank you Duncan for your reply.
> 
> Currently I am using 'double' for the computations.
> What type should I use for extended real in my intermediate computations?

I think it depends on compiler details.  On some compilers "long double" 
will get it, but I don't think there's a standard type that works on all 
compilers.  (In fact, the 80 bit type on Intel chips isn't necessarily 
supported on other hardware.)  R defines LDOUBLE in its header files and 
it is probably best to use that if you want to duplicate R results.

> The result will still be 'double' anyway right?

Yes, you do need to return type double.

Duncan Murdoch

> 
> 
> 
> On 10/09/2010 13:00, Duncan Murdoch wrote:
>> 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
>>>
> 
>  
> 
> ###
> UNIVERSITY OF CAPE TOWN 
> 
> This e-mail is subject to the UCT ICT policies and e-mail disclaimer published on our website at http://www.uct.ac.za/about/policies/emaildisclaimer/ or obtainable from +27 21 650 4500. This e-mail is intended only for the person(s) to whom it is addressed. If the e-mail has reached you in error, please notify the author. If you are not the intended recipient of the e-mail you may not use, disclose, copy, redirect or print the content. If this e-mail is not related to the business of UCT it is sent by the sender in the sender's individual capacity.
> 
> ###
>  
>