From renaud at mancala.cbio.uct.ac.za  Fri Sep 10 13:07:34 2010
From: renaud at mancala.cbio.uct.ac.za (Renaud Gaujoux)
Date: Fri, 10 Sep 2010 13:07:34 +0200
Subject: [Rd] Non identical numerical results from R code vs C/C++ code?
In-Reply-To: <4C8A0FC9.9060400@gmail.com>
References: <4C8A0C86.8030307@cbio.uct.ac.za> <4C8A0FC9.9060400@gmail.com>
Message-ID: <4C8A1176.9060906@cbio.uct.ac.za>

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?
The result will still be 'double' anyway right?



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.

###



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