Implementation in R of the Kahan summation algorithm

kahansum.R

# Kahan summation
# ref: http://rosettacode.org/wiki/Kahan_summation

# Try using the task rules
# 1. Do all arithmetic in decimal to a precision of six digits.
options(digits=6)

# 2. Write a function/method/procedure/subroutine/... to perform Kahan summation
# on an ordered collection of numbers, (such as a list of numbers).
kahansum <- function(x) {
    ks <- 0
    c <- 0
    for(i in 1:length(x)) {
        y <- x[i] - c
        kt <- ks + y
        c = (kt - ks) - y
        ks = kt
    }
    ks
}

# 3. Create the three numbers a, b, c equal to 10000.0, 3.14159, 2.71828 respectively.
a <- 10000.0
b <- 3.14159
c <- 2.71828

# 4. Show that the simple left-to-right summation, equivalent to (a + b) + c gives an answer of 10005.8
a + b + c
# [1] 10005.9

# 5. Show that the Kahan function applied to the sequence of values a, b, c results in the more precise answer of 10005.9
input <- c(a, b, c)
kahansum(input) 
# [1] 10005.9

# as there is no difference let's try the subtask
epsilon = 1.0
while ((1.0 + epsilon) != 1.0) {
    epsilon = epsilon / 2.0
}
epsilon
# [1] 1.11022e-16

a <- 1.0000000000000011
b <- epsilon
c <- -epsilon

# left-to-right summation
(a + b) + c
# [1] 1

# kahan summation
kahansum(c(a, b, c))
# [1] 1

# but, are they really the same or is an artifact?
(a + b) + c == kahansum(c(a, b, c))
# FALSE

# ok, then what is the difference?
((a + b) + c) - kahansum(c(a, b, c))
# [1] 2.22045e-16