[Rd] code for sum function
Berend Hasselman
bhh @end|ng |rom x@4@||@n|
Wed Feb 20 11:11:50 CET 2019
>
> On 2019-02-19 2:08 p.m., William Dunlap via R-devel wrote:
>> The algorithm does make a differece. You can use Kahan's summation
>> algorithm (https://en.wikipedia.org/wiki/Kahan_summation_algorithm) to
>> reduce the error compared to the naive summation algorithm. E.g., in R
>> code:
>>
>> naiveSum <-
>> function(x) {
>> s <- 0.0
>> for(xi in x) s <- s + xi
>> s
>> }
>> kahanSum <- function (x)
>> {
>> s <- 0.0
>> c <- 0.0 # running compensation for lost low-order bits
>> for(xi in x) {
>> y <- xi - c
>> t <- s + y # low-order bits of y may be lost here
>> c <- (t - s) - y
>> s <- t
>> }
>> s
>> }
>>
>>> rSum <- vapply(c(1:20,10^(2:7)), function(n) sum(rep(1/7,n)), 0)
>>> rNaiveSum <- vapply(c(1:20,10^(2:7)), function(n) naiveSum(rep(1/7,n)), 0)
>>> rKahanSum <- vapply(c(1:20,10^(2:7)), function(n) kahanSum(rep(1/7,n)), 0)
>>>
>>> table(rSum == rNaiveSum)
>>
>> FALSE TRUE
>> 21 5
>>> table(rSum == rKahanSum)
>>
>> FALSE TRUE
>> 3 23
If you use the vector c(1,10^100,1,-10^100) as input then
sum, naiveSum or kahanSum will all give an incorrect answer.
All return 0 instead of 2.
From the wikipedia page we can try the pseudocode given of the modification by Neumaier.
My R version (with a small correction to avoid cancellation?) is
neumaierSum <- function (x)
{
s <- 0.0
z <- 0.0 # running compensation for lost low-order bits
for(xi in x) {
t <- s + xi
if( abs(s) >= abs(xi) ){
b <- (s-t)+xi # intermediate step needed in R otherwise cancellation
z <- z+b # If sum is bigger, low-order digits of xi are lost.
} else {
b <- (xi-t)+s # intermediate step needed in R otherwise cancellation
z <- z+b # else low-order digits of sum are lost
}
s <- t
}
s+z # correction only applied once in the very end
}
testx <- c(1,10^100,1,-10^100)
neumaierSum(testx)
gives 2 as answer.
Berend Hasselman
More information about the R-devel
mailing list