[Rd] code for sum function

[Berend Hasselman]

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