[Rd] code for sum function

[Ben Bolker]

This SO question may be of interest:

https://stackoverflow.com/questions/38589705/difference-between-rs-sum-and-armadillos-accu/

  which points out that sum() isn't doing anything fancy *except* using
extended-precision registers when available.  (Using Kahan's algorithm
does come at a computational cost ...)

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
> 
> 
> Bill Dunlap
> TIBCO Software
> wdunlap tibco.com
> 
> 
> On Tue, Feb 19, 2019 at 10:36 AM Paul Gilbert <pgilbert902 using gmail.com> wrote:
> 
>> (I didn't see anyone else answer this, so ...)
>>
>> You can probably find the R code in src/main/ but I'm not sure. You are
>> talking about a very simple calculation, so it seems unlike that the
>> algorithm is the cause of the difference. I have done much more
>> complicated things and usually get machine precision comparisons. There
>> are four possibilities I can think of that could cause (small) differences.
>>
>> 0/ Your code is wrong, but that seems unlikely on such a simple
>> calculations.
>>
>> 1/ You are summing a very large number of numbers, in which case the sum
>> can become very large compared to numbers being added, then things can
>> get a bit funny.
>>
>> 2/ You are using single precision in fortran rather than double. Double
>> is needed for all floating point numbers you use!
>>
>> 3/ You have not zeroed the double precision numbers in fortran. (Some
>> compilers do not do this automatically and you have to specify it.) Then
>> if you accidentally put singles, like a constant 0.0 rather than a
>> constant 0.0D+0, into a double you will have small junk in the lower
>> precision part.
>>
>> (I am assuming you are talking about a sum of reals, not integer or
>> complex.)
>>
>> HTH,
>> Paul Gilbert
>>
>> On 2/14/19 2:08 PM, Rampal Etienne wrote:
>>> Hello,
>>>
>>> I am trying to write FORTRAN code to do the same as some R code I have.
>>> I get (small) differences when using the sum function in R. I know there
>>> are numerical routines to improve precision, but I have not been able to
>>> figure out what algorithm R is using. Does anyone know this? Or where
>>> can I find the code for the sum function?
>>>
>>> Regards,
>>>
>>> Rampal Etienne
>>>
>>> ______________________________________________
>>> R-devel using r-project.org mailing list
>>> https://stat.ethz.ch/mailman/listinfo/r-devel
>>
>> ______________________________________________
>> R-devel using r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-devel
>>
> 
> 	[[alternative HTML version deleted]]
> 
> ______________________________________________
> R-devel using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>