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