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@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Thomas Lumley
Professor of Biostatistics
University of Washington, Seattle
______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
