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

[Thomas Lumley]

On Fri, 10 Sep 2010, Duncan Murdoch wrote:

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

As a little more detail, 'long double' is in the C99 standard and seems to be fairly widely implemented, so code using it is likely to compile.   The Standard, as usual, doesn't define exactly what type it is, and permits it to be a synonym for 'double', so you may not get any extra precision.

On Intel chips it is likely to be the 80-bit type, but the Sparc architecture doesn't have any larger hardware type.  Radford Neal has recently reported much slower results on Solaris with long double, consistent with Wikipedia's statement that long double is sometimes a software-implemented 128-bit type on these systems.


>> 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.
>> 
>> ###
>>  
>
> ______________________________________________
> R-devel at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>

Thomas Lumley
Professor of Biostatistics
University of Washington, Seattle