From: Prof Brian Ripley <ripley_at_stats.ox.ac.uk>

Date: Sun 19 Feb 2006 - 10:58:56 GMT

Date: Sun 19 Feb 2006 - 10:58:56 GMT

> var(rep(0.02, 10))

[1] 1.337451e-35

> options(digits=18)

> sum(rep(0.02, 10))

[1] 0.19999999999999998

> sum(rep(0.02, 10)) -0.2

[1] -2.775557561562891e-17

> sum(rep(0.02, 10))/10 -0.02

[1] -3.469446951953614e-18

> z <- sum(rep(0.02, 10))/10 -0.02

> 10*z^2/9

[1] 1.3374513502689138e-35

## sum(int) can overflow, so convert here. if(is.integer(x)) x <- as.numeric(x) ## use one round of iterative refinement res <- sum(x)/n if(is.finite(res)) res + sum(x-res)/n else res

This is a well-known technique in numerical linear algebra, and improves the accuracy whilst doubling the cost. (This is about to be replaced by an internal function to allow the intermediate result to be stored in a long double.)

Note the is.finite(res) there. R works with the extended IEC60059 (aka IEEE754) quantities of Inf, -Inf and NaN (NA being a special NaN). The rearranged computations do not work correctly for those quantities. So although they can be more 'efficient' (in terms of flops), they have to be supplemented by some other calculation to ensure that the specials are handled correctly. Remember once again that we get both speed and accuracy advantages by keeping computations in FPU registers, so complicating the code has considerable cost.

R-devel will shortly use long doubles for critical intermediate results and iterative refinement for calculations of means. This may be slower but it would be an extreme case in which the speed difference was detectable. Higher accuracy has a cost too: there are several packages (dprep and mclust are two) whose results are critically dependent on fine details of computations and will for example infinite-loop if an optimized BLAS is used on AMD64.

-- Brian D. Ripley, ripley@stats.ox.ac.uk Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UK Fax: +44 1865 272595 ______________________________________________ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-develReceived on Sun Feb 19 22:01:39 2006

*
This archive was generated by hypermail 2.1.8
: Mon 20 Feb 2006 - 03:21:41 GMT
*