Kahan summation

Note to language implementers

This task touches on the details of number representation that is often hidden for convenience in many programming languages. You may need to understand the difference between precision in calculations and precision in presenting answers for example; or understand decimal arithmetic as distinct from the usual binary based floating point schemes used in many languages.

The Kahan summation algorithm is a method of summing a series of numbers represented in a limited precision in a way that minimises the loss of precision in the result.

The task is to follow the previously linked Wikipedia articles algorithm and its worked example by:

Do all arithmetic in decimal to a precision of six digits.
Write a function/method/procedure/subroutine/... to perform Kahan summation on an ordered collection of numbers, (such as a list of numbers).
Create the three numbers a, b, c equal to 10000.0, 3.14159, 2.71828 respectively.
Show that the simple left-to-right summation, equivalent to (a + b) + c gives an answer of 10005.8
Show that the Kahan function applied to the sequence of values a, b, c results in the more precise answer of 10005.9

If your language does not have six digit decimal point, but does have some fixed precision floating point number system then state this and try this alternative task using floating point numbers:

Follow the same steps as for the main task description above but for the numbers a, b, c use the floating-point values 1.0, epsilon, -epsilon respectively where epsilon is determined by its final value after the following:

<lang python>epsilon = 1.0 while 1.0 + epsilon != 1.0:

   epsilon = epsilon / 2.0</lang>

The above should ensure that (a + b) + c != 1.0 whilst their Kahan sum does equal 1.0 . Only if this is not the case are you then free to find three values with this property of your own.

If your language has decimals, but not six digits of precision then you may follow the outline given in Kahan_summation#Python:_Arbitrary_precision_Decimal which uses ideas from the floating point case above.

In general

Show all output on this page.
If the floating point calculations are used then answers may depend on the hardware platform or compiler/interpreter tool-chain used for a language.
Slight deviations from the task description should be explained and the the subsequent difference in output explained.
All examples should have constants chosen to clearly show the benefit of Kahan summing!

ALGOL 68

Defines and uses a 6 digit float decimal type to solve the main "1000.0 + pi + e" task. <lang algol68># Kahan summation - Algol 68 doesn't have a float decimal mode, however as #

only 6 digits of precision are required and assuming INT is 32 bits, we can #
emulate this fairly easily #

MODE DEC is a float decimal with 6 digits #
We normalise the mantissa so that the first digit is non-zero if #
the number is not 0 #
so e.g. 1 is represented as 100000 E0 #
10 is represented as 100000 E1 #
etc. #
we only provide enough operations for the task #

MODE DEC = STRUCT( BOOL negative, INT mantissa, INT exponent );

constructs a string representation of a DEC number #

OP TOSTRING = ( DEC number )STRING:

   IF mantissa OF number = 0
   THEN
       # number is zero #
       " 0.00000E +0"
   ELSE
       # non-zero number #
       STRING result = whole( mantissa OF number, 0 );
       ( IF negative OF number
         THEN
             # number is negative - show leading "-" sign #
             "-"
         ELSE
             # positive number - show leading space #
             " "
         FI
       + result[ 1 : 1 ]
       + "."
       + result[ 2 :   ]
       + "E"
       + whole( exponent OF number, 3 )
       )
   FI

negates a DEC number #

OP - = ( DEC a )DEC: ( NOT negative OF a, mantissa OF a, exponent OF a );

adds two DEC numbers #

OP + = ( DEC a, DEC b )DEC:

   IF exponent OF a < exponent OF b
   THEN
       # the left operand has a smaller exponent than the right,             #
       # reverse the operands                                                #
       b + a

   ELSE
       # the left operand's exponent is at least as big as the right's       #

       # removes a digit from a DEC number with rounding, the exponent is    #
       # unchanged                                                           #
       OP REDUCE = ( DEC a )DEC: ( negative OF a
                                 , ( ( mantissa OF a OVER 10 )
                                     + IF mantissa OF a MOD 10 > 4 THEN 1 ELSE 0 FI
                                   )
                                 , exponent OF a
                                 );
       # adds an extra digit to a DEC number, exponent is unchanged          #
       OP EXTEND = ( DEC a )DEC: ( negative OF a, mantissa OF a * 10, exponent OF a );
       # signs the mantissa                                                  #
       OP APPLYSIGN = ( DEC a )DEC:( FALSE, IF negative OF a THEN - mantissa OF a ELSE mantissa OF a FI, exponent OF a );
       # unsigns the mantissa and sets negative accordingly                  #
       OP EXTRACTSIGN = ( DEC a )DEC: ( mantissa OF a < 0, ABS ( mantissa OF a ), exponent OF a );

       # give each number an extra digit                                     #
       DEC result  := EXTEND a;
       DEC operand := EXTEND b;

       # scale the values so that they have the same exponent                #

       WHILE exponent OF operand < exponent OF result
       DO
           operand := REDUCE operand;
           exponent OF operand PLUSAB 1
       OD;

       # sign the mantissas                                                  #
       result  := APPLYSIGN result;
       operand := APPLYSIGN operand;

       # add                                                                 #
       mantissa OF result PLUSAB mantissa OF operand;

       # remove the extra digit, round and unsign the mantissa               #
       result := REDUCE EXTRACTSIGN result;

       # if the mantissa is < 100000 and not 0, increase the mantissa and    #
       # reduce the exponent so that the first digit is not 0                #
       IF mantissa OF result /= 0
       THEN
           # non-zero mantissa - scale if necessary                          #
           WHILE mantissa OF result < 100000
           DO
               mantissa OF result TIMESAB 10;
               exponent OF result MINUSAB  1
           OD
       FI;

       result

   FI # + #

subtracts two DEC numbers #

OP - = ( DEC a, DEC b )DEC: a + ( - b );

performs Kahan summation on an array of DEC numbers #
The numbers must be ordered highest to lowest #
Algorithm as per the Wikipaedia page, except we start with sum #
set to the first element of the list, not 0 #

OP KAHANSUM = ( []DEC a )DEC:

   BEGIN

       DEC sum := a[LWB a];
       DEC c   := ( FALSE, 0, 0 );

       FOR i FROM LWB a + 1 TO UPB a
       DO
           DEC y := a[i] - c;
           DEC t := sum + y;
           c     := ( t - sum ) - y;
           sum   := t
       OD;

       sum
   END # KAHANSUM #

test the Kahan summation and 6 digit decimal floating point #

main: (

   # construct DEC values a = 10000.0, b = 3.14159, c = 2.71828              #

   DEC a := ( FALSE, 100000, 4 );
   DEC b := ( FALSE, 314159, 0 );
   DEC c := ( FALSE, 271828, 0 );

   # simple sum #
   print( ( ( "Simple: " + TOSTRING a + " + " + TOSTRING b + " + " + TOSTRING c + " = " + TOSTRING ( ( a + b ) + c ) )
          , newline
          )
        );
   # Kahan #
   []DEC list = ( a, b, c );
   print( ( ( "Kahan : " + TOSTRING a + " + " + TOSTRING b + " + " + TOSTRING c + " = " + TOSTRING KAHANSUM list )
          , newline
          )
        );

   SKIP

) </lang>

Output:

Simple: 1.00000E +4 + 3.14159E +0 + 2.71828E +0 = 1.00058E +4 Kahan : 1.00000E +4 + 3.14159E +0 + 2.71828E +0 = 1.00059E +4

EchoLisp

No fixed point addition. EchoLisp floating point numbers are always stored as double precision floating point numbers, following the international IEEE 754 standard. <lang scheme>

using floating point arithmetic the maximum number of decimals is 17. Floating point arithmetic is not always 100% accurate even with simple data

(+ 0.2 0.1 0.3) → 0.6000000000000001 ;; (1)

👀 Will Kahan method do better ?? Kahan procedure

(define ( K+ nums ( sum 0.0) (c 0.0) (y) (t)) (while (!null? nums) (set! y (- (first nums) c)) (set! t (+ sum y)) (set! c (- (- t sum) y)) (set! sum t) (set! nums (rest nums))) sum)

👏 👏 Kahan improves on the above (1)

(K+ '( 0.2 0.1 0.3)) → 0.6

using epsilon such as (1.0 + epsilon) = 1.0 compute 1 + epsilon - epsilon

(define eps 1.0)

    (while (!= (+ 1.0 eps) 1.0) (set! eps (// eps 2.0)))

eps → 1.1102230246251565e-16

(define a 1.0 ) (define b eps) (define c (- b))

(writeln 'standard-add= (+ a b c)) (writeln 'Kahan-add= (K+ (list a b c)))

   → standard-add= 0.9999999999999999 
   → Kahan-add= 1

NB. The 17-nth decimal we gain using Kahan method will probably be lost in subsequent calculations. Kahan algorithm will be useful for fixed-point decimal calculations, which EchoLisp does not have.

</lang>

F#

This example is incomplete. Slight deviations from the task description should be explained and the the subsequent difference in output explained. All examples should have constants chosen to clearly show the benefit of Kahan summing! - The J-language example for example explains its differences well. Please ensure that it meets all task requirements and remove this message.

Solving the alternative task with an recursive algorithm and IEEE 754 double precision datatype. <lang fsharp>let ε =

   let rec δ λ = 
       match λ+1.0 with 
       | 1.0 -> λ
       | _ -> δ(λ/2.0)
   δ(1.0)

let Σ numbers =

   let rec Σ numbers σ c = 
       match numbers with
       | [] -> σ
       | number :: rest -> 
           let y = number + c
           let t = σ + y
           let z = ((t - σ) - y)
           Σ rest t z
   Σ numbers 0.0 0.0

[<EntryPoint>] let main argv =

   let (a, b, c) = ( 1.0, ε, -ε )
   let input = [ a; b; c ]
   printfn "ε:     %e" ε
   printfn "Sum:   %e" (input |> List.sum)
   printfn "Kahan: %e" (input |> Σ)
   0</lang>

Output:

ε:     1.110223e-016
Sum:   1.000000e+000
Kahan: 1.000000e+000

Fortran

As the person who provided the worked example for the Wikipaedia article, I am hoist by my own petard! So, consider the IBM1620, a decimal computer, where the user could specify the number of digits to be used in arithmetic: it shall be six decimal digits instead of the default eight. The following source file is Second Fortran (of 1958). Variables were not usually declared (there are no INTEGER or REAL statements), so those whose name starts with I,...,N were integers (known as "fixed-point"), otherwise floating-point. Array bounds were not checked, and a subprogram receiving an array as a parameter needed to know merely that it was an array with the bound specified for one-dimensional arrays unimportant. Standard output was via the console typewriter (a sturdy device), via a TYPE statement. Text literals were available only via the "H"-format code, which consisted of a count before the H, followed by exactly that number of characters, or else... A lot of time and patience attended miscounts. All text is of course in capitals only. Indentation is a latter-day affectation: many decks of cards were prepared with minimal spacing. <lang Fortran>

     REAL FUNCTION SUMC(A,N)

COMPENSATED SUMMATION. C WILL NOT STAY ZERO, DESPITE MATHEMATICS.

      DIMENSION A(12345)
       S = 0.0
       C = 0.0
       DO 1 I = 1,N
         Y = A(I) - C
         T = S + Y
         C = (T - S) - Y
         S = T
   1   CONTINUE
       SUMC = S
     END
     DIMENSION A(3)
     A(1) = 10000.0
     A(2) = 3.14159
     A(3) = 2.71828
     TYPE 1, A(1) + A(2) + A(3)
     TYPE 1, SUMC(A,3)
   1 FORMAT (6HSUM = ,F12.1)
     END

</lang> Alas, I no longer have access to an IBM1620 (or an emulator) whereby to elicit output. Exactly this source is acceptable to a F90/95 compiler, but it on current computers will use binary arithmetic and a precision not meeting the specification. Fortran as a language does not have a "decimal" type (as say in Cobol or Pl/i) but instead the compiler uses the arithmetic convenient for the cpu it is producing code for. This may be in base ten as with the IBM1620 and some others, or base two, or base eight, or base sixteen. Similarly, the system's standard word size may be 16, 18 (PDP 15), 32, 48 (B6700), 64, or even 128 bits. Nor is that the end of the variations. The B6700 worked in base eight which meant that a floating-point number occupied 47 bits. The 48'th bit was unused in arithmetic, including that of integers. All of this means that number size specifications such as REAL*8 are not as straightforward as one might imagine, which is why REAL and DOUBLE are preferred over REAL*4 and REAL*8. And, even on binary computers where only simple sizes such as 32 or 64 bit numbers are offered there are still problems because with floating-point there is the question of assumed-one leading-bit for normalised numbers. The IBM pc and its followers uses this for 32- and 64-bit floating-point, but not for 80-bit floating point.

So much for the first phase.

Go

Go has no floating point decimal type. Its floating point types are float64 and float32, implementations of IEEE 754 binary64 and binary32 formats. Alternative task: <lang go>package main

import "fmt"

type float float64

func kahan(s ...float) float {

   var tot, c float
   for _, x := range s {
       y := x - c
       t := tot + y
       c = (t - tot) - y
       tot = t
   }
   return tot

}

func epsilon() float {

   ε := float(1)
   for 1+ε != 1 {
       ε /= 2
   }
   return ε

}

func main() {

   a := float(1)
   b := epsilon()
   c := -b
   fmt.Println("Left associative:", a+b+c)
   fmt.Println("Kahan summation: ", kahan(a, b, c))
   fmt.Println("Epsilon:         ", b)

}</lang>

Output:

Left associative: 0.9999999999999999
Kahan summation:  1
Epsilon:          1.1102230246251565e-16

With float defined as float32,

Left associative: 0.99999994
Kahan summation:  1
Epsilon:          5.9604645e-08

J

J does not implement 6 digit floating point addition (though the programmer could of course emulate such a thing - if it were specified precisely enough). But J does have IEEE-754 floating point numbers. So...

Required code (for the current draft of this task):

<lang J>epsilon=:3 :0

 epsilon=. 1.0
 while. 1.0 ~: 1.0 + epsilon do.
   epsilon=. epsilon % 2.0
 end.

)

a=: 1.0 b=: epsilon c=: -epsilon

KahanSum=:3 :0

 input=. y
 c=. 0.0
 sum=. 0.0
 for_i. i.#input do.
   y=. (i{input) - c
   t=. sum + y
   c=. (t - sum) - y
   sum=. t
 end.

)</lang>

Required results:

<lang J> (a+b)+c 1

  KahanSum a,b,c

1</lang>

There are several issues here, but the bottom line is that several assumptions behind the design of this task conflict with the design of J.

But let's take them one by one.

First, the value of epsilon:

<lang J> epsilon 5.68434e_14</lang>

This is because J's concept of equality, in the context of floating point numbers, already includes an epsilon. Floating point numbers are inherently imprecise, and this task has assumed that the language has assumed that floating point numbers are precise. So let's implement something closer to this previously unstated assumption:

<lang J>epsilon=:3 :0

 epsilon=. 1.0
 while. 1.0 (~:!.0) 1.0 + epsilon do.
   epsilon=. epsilon % 2.0
 end.

)</lang>

Now we have a smaller value for this explicit value of epsilon:

<lang J> epsilon 1.11022e_16</lang>

With this new value for epsilon, let's try the problem again:

<lang J> (a+b)+c 1

  KahanSum a,b,c

1</lang>

Oh dear, we still are not failing in the manner clearly specified as a task requirement.

Well, again, the problem is that the task has assumed that the language is treating floating point values as precise when their actual accuracy is less than their precision. This time, the problem is in the display of the result. So let's also override that mechanism and ask J to display results to 16 places of precision:

<lang J> ":!.16 (a+b)+c 0.9999999999999999

  ":!.16 KahanSum a,b,c

1</lang>

Voila!

PARI/GP}

No decimals here; the example below uses 64-bit binary arithmetic.

This is a partial solution only; I haven't found appropriate values that fail for normal addition. (Experimentation should produce these, in which case the existing code should work with the inputs changed.) <lang parigp>Kahan(v)=my(s=0.,c=0.,y,t);for(i=1,#v,y=v[i]-c;t=s+y;c=t-s-y;s=t);s; epsilon()=my(e=1.); while(e+1. != 1., e>>=1); e; \p 19 e=epsilon(); Kahan([1.,e,-e])</lang>

Output:

%1 = 1.000000000000000000

PHP

Script

<lang php><?php // Main task

// 1. Do all arithmetic in decimal to a precision of six digits. // PHP does not have a way of restricting overall precision

// 2. Write a function/method/procedure/subroutine/... to perform Kahan // summation on an ordered collection of numbers, (such as a list of numbers). ) function kahansum($input) {

   $sum = $c = 0;
   foreach($input as $item) {
       $y = $item + $c;
       $t = $sum + $y;
       $c = ($t - $sum) - $y;
       $sum = $t;
   }
   return $sum;

}

// 3. Create the three numbers a, b, c equal to 10000.0, 3.14159, 2.71828 // respectively. $input = array(10000.0, 3.14159, 2.71828); list($a, $b, $c) = $input;

// 4. show that the simple left-to-right summation, equivalent to (a + b) + c // gives an answer of 10005.8 ) $sumabc = ($a + $b) + $c; echo "Main task - Left to right summation: "; echo sprintf("%6.1f", $sumabc).PHP_EOL; // 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 echo "Main task - Kahan summation: "; echo sprintf("%6.1f", kahansum($input)).PHP_EOL; // 10005.9

// Let's use the substask $epsilon = 1.0; while ((1.0 + $epsilon) != 1.0) {

   $epsilon /= 2.0;

} echo "Trying the subtask".PHP_EOL."epsilon: "; echo $epsilon.PHP_EOL; // 1.1102230246252E-16

$a = 1.0; $b = $epsilon; $c = -$epsilon;

$input = array($a, $b, $c);

// left-to-right summation $sumabc = ($a + $b) + $c; echo "Sub task - Left to right summation: "; echo sprintf("%.1f", $sumabc).PHP_EOL;

// kahan summation echo "Sub task - Kahan summation: "; echo sprintf("%.1f", kahansum($input)).PHP_EOL;

// but, are they really the same or is an artifact? echo "Are the results the same?".PHP_EOL; var_dump((($a + $b) + $c) === kahansum($input));

// ok, then what is the difference? echo "Difference between the operations: "; $diff = (($a + $b) + $c) - kahansum($input); echo $diff.PHP_EOL; </lang>

Script output

<lang bash> $ php kahansum.php Main task - Left to right summation: 10005.9 Main task - Kahan summation: 10005.9 Trying the subtask epsilon: 1.1102230246252E-16 Sub task - Left to right summation: 1.0 Sub task - Kahan summation: 1.0 Are the results the same? bool(false) Difference between the operations: 1.1102230246252E-16 </lang>

Python

Python: Decimal

<lang python>>>> from decimal import * >>> >>> getcontext().prec = 6 >>> >>> def kahansum(input):

   summ = c = 0
   for num in input:
       y = num - c
       t = summ + y
       c = (t - summ) - y
       summ = t
   return summ

>>> a, b, c = [Decimal(n) for n in '10000.0 3.14159 2.71828'.split()] >>> a, b, c (Decimal('10000.0'), Decimal('3.14159'), Decimal('2.71828')) >>> >>> (a + b) + c Decimal('10005.8') >>> kahansum([a, b, c]) Decimal('10005.9') >>> >>> >>> sum([a, b, c]) Decimal('10005.8') >>> # it seems Python's sum() doesn't make use of this technique. >>> >>> # More info on the current Decimal context: >>> getcontext() Context(prec=6, rounding=ROUND_HALF_EVEN, Emin=-999999, Emax=999999, capitals=1, clamp=0, flags=[Inexact, Rounded], traps=[InvalidOperation, DivisionByZero, Overflow]) >>> >>> >>> ## Lets try the simple summation with more precision for comparison >>> getcontext().prec = 20 >>> (a + b) + c Decimal('10005.85987') >>> </lang>

Python: Floats

This was written as proof that the floating-point sub-task could work and is better off displayed, so...

<lang python>>>> eps = 1.0 >>> while 1.0 + eps != 1.0: eps = eps / 2.0

>>> eps 1.1102230246251565e-16 >>> (1.0 + eps) - eps 0.9999999999999999 >>> kahansum([1, eps, -eps]) 1.0 >>> >>> >>> # Info on this implementation of floats >>> import sys >>> sys.float_info sys.float_info(max=1.7976931348623157e+308, max_exp=1024, max_10_exp=308, min=2.2250738585072014e-308, min_exp=-1021, min_10_exp=-307, dig=15, mant_dig=53, epsilon=2.220446049250313e-16, radix=2, rounds=1) >>> </lang> Note that the kahansum function from the decimal example is duck typed and adjusts to work with the number type used in its argument list.

Note also that the math.fsum function employs an even more precise algorithm for the summation of floating point numbers.

Python: Arbitrary precision Decimal

Some languages have a decimal type, but cannot alter its precision to six digits. This example mimmicks this case byusing the default Python decimal precision and a variant of epsilon finding that divides by ten instead of two. <lang python>>>> from decimal import localcontext, Decimal >>> >>> with localcontext() as ctx: one, ten = Decimal('1.0'), Decimal('10') eps = one while one + eps != one: eps = eps / ten print('eps is:', eps) print('Simple sum is:', (one + eps) - eps) print('Kahan sum is:', kahansum([one, eps, -eps]))

eps is: 1E-28 Simple sum is: 0.9999999999999999999999999999 Kahan sum is: 1.000000000000000000000000000 >>> </lang>

R

R can only limit the number of digits being displayed, not the ones used in calculation. In fact one of the base types is numeric, implemented as standard floating point number.

Therefore, trying to use the main rules will not work as we can see in the code below

an implementation of the Kahan summation algorithm

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

}

The three numbers a, b, c equal to 10000.0, 3.14159, 2.71828 respectively.

a <- 10000.0 b <- 3.14159 c <- 2.71828

just to make sure, let's look at the classes of these three numbers

sapply(c(a, b, c), FUN=Class)

[1] "numeric" "numeric" "numeric"

The simple left-to-right summation: (a + b) + c

(a + b) + c

[1] 10005.86

The Kahan summation applied to the values a, b, c

input <- c(a, b, c) kahansum(input)

[1] 10005.86

</lang>

Apparently there is no difference between the two approaches. So let's try the alternative steps given in the task.

We first calculate epsilon <lang r> epsilon = 1.0 while ((1.0 + epsilon) != 1.0) {

   epsilon = epsilon / 2.0

} epsilon

[1] 1.11022e-16

</lang>

and use it to test the left-to-right summation and the Kahan function

<lang r> a <- 1.0 b <- epsilon c <- -epsilon

left-to-right summation

(a + b) + c

[1] 1

kahan summation

kahansum(c(a, b, c))

[1] 1

</lang>

It might seem that again there is no difference, but let's explore a bit more and see if both results are really the same

<lang r> (a + b) + c == kahansum(c(a, b, c))

FALSE

ok, then what is the difference?

((a + b) + c) - kahansum(c(a, b, c))

[1] -1.110223e-16

</lang>

The absolute value of the difference, is very close to the value we obtained for epsilon

Just to make sure we are not fooling ourselves, let's see if the value of epsilon is stable using different divisors for the generating function

machine epsilon

mepsilon <- function(divisor) {

 guess <- 1.0
 while ((1.0 + guess) != 1.0) {
   guess <- guess / divisor
 }
 guess

}

let's try from 2 to 1000

epsilons <- sapply(2:1000, FUN=mepsilon) summary(epsilons)

Min. 1st Qu. Median Mean 3rd Qu. Max.
2.439e-19 1.905e-18 5.939e-18 1.774e-17 2.238e-17 1.110e-16

</lang>

It would seem that it makes a differences what divisor we are using, let's look at the distribution using a text plotting library

<lang r> library(txtplot)

txtboxplot(epsilons)

0 1e-17 2e-17 3e-17 4e-17 5e-17
|---+-----------+----------+-----------+-----------+-----------+--------|
+----+------------------+
--| | |-------------------------------------
+----+------------------+

txtdensity(epsilons)

6e+16 +--+---------+----------+----------+---------+----------+----------+
| ** |
| *** |
5e+16 + * ** +
| * |
4e+16 + * +
| * |
| * |
3e+16 + ** +
| * |
2e+16 + ** +
| ** |
| *** |
1e+16 + **** +
| ************* |
0 + ************************************ +
+--+---------+----------+----------+---------+----------+----------+
0 2e-17 4e-17 6e-17 8e-17 1e-16

</lang>

Definitely the epsilon generating function is not stable and gives a positively skewed (right-tailed) distribution.

We could consider using the median value of the series of epsilons we have estimated, but because the precision for the base numeric type (class) in R is ~16 decimals, using that value will be in practice indistinguishable from using zero.

<lang r> epsilon <- median(epsilons) epsilon

[1] 5.939024e-18

a <- 1.0 b <- epsilon c <- -epsilon

left-to-right summation

(a + b) + c

[1] 1

kahan summation

kahansum(c(a, b, c))

[1] 1

are they the same?

(a + b) + c == kahansum(c(a, b, c))

TRUE

</lang>

Racket

Racket doesn't have arbitrary fixed precision numbers, but we can use single precision float point numbers that have approximately 8 decimal places. To reproduce the original task we have to replace the 10000.0 with 1000000.0.

We then compare this result to the double precision result of the same numbers. The double precision numbers have almost 16 decimal places, so with 1000000.0 the usual summation nd the Kahan summation give the same result.

Finally we try the alternative task version for single precision float point numbers. <lang Racket>#lang racket

(define (sum/kahan . args)

 (define-values (sum c)
   (for/fold ([sum 0] [c 0]) ([num args])
     (define y (- num c))
     (define t (+ sum y))
     (values t (- (- t sum) y))))
 sum)

(displayln "Single presition flonum") (+ 1000000.0f0 3.14159f0 2.71828f0) (sum/kahan 1000000.0f0 3.14159f0 2.71828f0)

(displayln "Double presition flonum") (+ 1000000.0 3.14159 2.71828) (sum/kahan 1000000.0 3.14159 2.71828)</lang>

Output:

Single presition flonum
1000005.8f0
1000005.9f0
Double presition flonum
1000005.85987
1000005.85987

Alternative task version <lang Racket>(define epsilon.f0

 (let loop ([epsilon 1.f0])
   (if (= (+ 1.f0 epsilon) 1.f0)
       epsilon
       (loop (/ epsilon 2.f0)))))

(displayln "Alternative task, single precision flonum") epsilon.f0

(+ 1.f0 epsilon.f0 (- epsilon.f0)) (sum/kahan 1.f0 epsilon.f0 (- epsilon.f0))</lang>

Output:

Alternative task, single precision flonum
5.9604645f-008
0.99999994f0
1.0f0

REXX

Programming note: It wasn't clear what precision to use for the 2^nd part of this task, so an arbitrary precision
of 30 (decimal digits) was chosen. The default precision for REXX is 9 decimal digits.

vanilla version

<lang rexx>/*REXX program demonstrates simple addition versus using Kahan summation*/ numeric digits 6 /*use numeric (decimal) digits=6.*/ call show 10000.0, 3.14169, 2.71828 /*invoke SHOW to sum and display.*/

numeric digits 30 epsilon=1.0

                     do  while 1.0+epsilon \= 1.0
                     epsilon=epsilon / 2.0
                     end   /*while*/

call show 1.0, epsilon, -epsilon /*invoke SHOW to sum and display.*/ exit /*stick a fork in it, we're done.*/ /*──────────────────────────────────KAHAN subroutine────────────────────*/ kahan: procedure; sum=0; c=0

                             do j=1  for arg()       /*do for each arg.*/
                             y=arg(j)-c              /*sub C from arg. */
                             t=sum+y                 /*use a temp sum. */
                             c=(t-sum)-y             /*same as: t-sum-y*/
                             sum=t                   /*define the sum. */
                             end   /*j*/

return sum /*──────────────────────────────────SHOW subroutine─────────────────────*/ show: procedure; parse arg a,b,c /*obtain the args.*/ say 'decimal digits =' digits() /*show # dec digs.*/ say ' a = ' left(, a>=0) a /*show A justified*/ say ' b = ' left(, b>=0) b /* " B " */ say ' c = ' left(, c>=0) c /* " C " */ say 'simple summation of a,b,c = ' (a+b)+c /*same as a+b+c */ say 'Kahan summation of a,b,c = ' kahan(a,b,c) /*sum via Kahan. */ say; say copies('▒',70); say /*display a fence.*/ return</lang> Output note: This is a work in progress and attempts are being made to determine a suitable/workable epsilon.

output using PC/REXX and also Personal REXX:

decimal digits = 6
   a =  10000.0
   b =  3.14169
   c =  2.71828
simple summation of a,b,c =  10005.8
Kahan  summation of a,b,c =  10005.9

▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒

decimal digits = 30
   a =   1.0
   b =   0.00000000000000000000000000000315544362088404722164691426133
   c =  -0.00000000000000000000000000000315544362088404722164691426133
simple summation of a,b,c =  1.00000000000000000000000000000
Kahan  summation of a,b,c =  1.00000000000000000000000000000

▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒

output using R4 and also ROO:

decimal digits = 6
   a =   10000.0
   b =   3.14169
   c =   2.71828
simple summation of a,b,c =  10005.8
Kahan  summation of a,b,c = 10005.9

▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒

decimal digits = 30
   a =   1.0
   b =   0.00000000000000000000000000000315544362088404722164691426133
   c =  -0.00000000000000000000000000000315544362088404722164691426133
simple summation of a,b,c =  0.999999999999999999999999999997
Kahan  summation of a,b,c =  1.00000000000000000000000000000

▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒

output using Regina 3.3 and earlier:

decimal digits = 6
   a =   10000.0
   b =   3.14169
   c =   2.71828
simple summation of a,b,c =  10005.8
Kahan  summation of a,b,c =  10005.9

▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒

decimal digits = 30
   a =   1.0
   b =   3.15544362088404722164691426133E-30
   c =  -3.15544362088404722164691426133E-30
simple summation of a,b,c =  0.999999999999999999999999999997
Kahan  summation of a,b,c =  1.00000000000000000000000000000

▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒

output using Regina 3.4 and later (latest at this time is 3.9.0):

decimal digits = 6
   a =   10000.0
   b =   3.14169
   c =   2.71828
simple summation of a,b,c =  10005.8
Kahan  summation of a,b,c =  10005.9

▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒

decimal digits = 30
   a =   1.0
   b =   3.15544362088404722164691426133E-30
   c =  -3.15544362088404722164691426133E-30
simple summation of a,b,c =  1.00000000000000000000000000000
Kahan  summation of a,b,c =  1.00000000000000000000000000000

▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒

tweaked version

The following tweaked REXX version causes Regina (3.4 and later) to work properly/ <lang rexx>/*REXX program demonstrates simple addition versus using Kahan summation*/ numeric digits 6 /*use numeric (decimal) digits=6.*/ call show 10000.0, 3.14169, 2.71828 /*invoke SHOW to sum and display.*/

numeric digits 30 epsilon=1.0

                     do  while  1.0+epsilon \= 1.0
                     epsilon=epsilon / 2.0
                     end   /*while*/
                                      /* [↓]  for Regina 3.4 and later.*/

numeric digits digits()+2 /*bump the precision by 2 digits.*/ call show 1.0, epsilon, -epsilon /*invoke SHOW to sum and display.*/ exit /*stick a fork in it, we're done.*/ /*──────────────────────────────────KAHAN subroutine────────────────────*/ kahan: procedure; sum=0; c=0

                             do j=1  for arg()       /*do for each arg.*/
                             y=arg(j)-c              /*sub C from arg. */
                             t=sum+y                 /*use a temp sum. */
                             c=(t-sum)-y             /*same as: t-sum-y*/
                             sum=t                   /*define the sum. */
                             end   /*j*/

return sum /*──────────────────────────────────SHOW subroutine─────────────────────*/ show: procedure; parse arg a,b,c /*obtain the args.*/ say 'decimal digits =' digits() /*show # dec digs.*/ say ' a = ' left(, a>=0) a /*show A justified*/ say ' b = ' left(, b>=0) b /* " B " */ say ' c = ' left(, c>=0) c /* " C " */ say 'simple summation of a,b,c = ' (a+b)+c /*same as a+b+c */ say 'Kahan summation of a,b,c = ' kahan(a,b,c) /*sum via Kahan. */ say; say copies('▒',70); say /*display a fence.*/ return</lang> output for Regina 3.4 and later versions:

decimal digits = 6
   a =    10000.0
   b =    3.14169
   c =    2.71828
simple summation of a,b,c =  10005.8
Kahan  summation of a,b,c =  10005.9

▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒

decimal digits = 32
   a =    1.0
   b =    3.15544362088404722164691426133E-30
   c =   -3.15544362088404722164691426133E-30
simple summation of a,b,c =  1.0000000000000000000000000000001
Kahan  summation of a,b,c =  1.0000000000000000000000000000000

▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒

Tcl

Tcl: Floats

First, using native floating point we see the same epsilon value as other languages using float64:

<lang Tcl># make {+ - * /} etc available as commands, for easier expressions namespace path ::tcl::mathop

find epsilon with native floating point:

proc epsilon {} {

   set e 1.0
   while {1 + $e != 1} {
       set e   [/ $e 2]
   }
   return $e

}

kahan sum with native floats:

proc kahansum {args} {

   set sum 0.0
   set c   0.0
   foreach i $args {
       set y   [- $i $c]
       set t   [+ $sum $y]
       set c   [- [- $t $sum] $y]
       set sum $t
   }
   return $sum

}

puts "Native floating point:" puts "\tEpsilon is: [set e [epsilon]]" puts "\tAssociative sum: [expr {1.0 + $e - $e}]" puts "\tKahan sum: [kahansum 1.0 $e -$e]"</lang>

Epsilon is: 1.1102230246251565e-16
Associative sum: 0.9999999999999999
Kahan sum: 1.0

Tcl: Decimals

Library: Tcllib (Package: math::decimal)

For the decimal part of the exercise we can use a the Tcllib library math::decimal. Note how similar the implementation of Kahan sum is: the only changes are fromstr and tostr.

The last stanza exercises the decimal package's different rounding modes, to see what happens there:

<lang Tcl>package require math::decimal namespace path ::math::decimal

proc kahansum {args} {

   set sum [fromstr 0.0]
   set c   [fromstr 0.0]
   foreach i $args {
       set i [fromstr $i]
       set y   [- $i $c]
       set t   [+ $sum $y]
       set c   [- [- $t $sum] $y]
       set sum $t
   }
   return [tostr $sum]

}

proc asum {args} {

   set sum [fromstr 0.0]
   foreach a $args {
       set sum [+ $sum [fromstr $a]]
   }
   return [tostr $sum]

}

setVariable precision 6 set a 10000.0 set b 3.14159 set c 2.71828

foreach rounding {half_even half_up half_down down up floor ceiling} {

   setVariable rounding $rounding
   puts "Rounding mode: $rounding"
   puts "\tAssociative sum $a + $b + $c: [asum $a $b $c]"
   puts "\tKahan       sum $a + $b + $c: [kahansum $a $b $c]"

}</lang>

The results are a little surprising:

Output:

Rounding mode: half_even
        Associative sum 10000.0 + 3.14159 + 2.71828: 10005.8
        Kahan       sum 10000.0 + 3.14159 + 2.71828: 10005.9
Rounding mode: half_up
        Associative sum 10000.0 + 3.14159 + 2.71828: 10005.8
        Kahan       sum 10000.0 + 3.14159 + 2.71828: 10005.9
Rounding mode: half_down
        Associative sum 10000.0 + 3.14159 + 2.71828: 10005.8
        Kahan       sum 10000.0 + 3.14159 + 2.71828: 10005.9
Rounding mode: down
        Associative sum 10000.0 + 3.14159 + 2.71828: 10005.8
        Kahan       sum 10000.0 + 3.14159 + 2.71828: 10005.8
Rounding mode: up
        Associative sum 10000.0 + 3.14159 + 2.71828: 10006.0
        Kahan       sum 10000.0 + 3.14159 + 2.71828: 10005.9
Rounding mode: floor
        Associative sum 10000.0 + 3.14159 + 2.71828: 10005.8
        Kahan       sum 10000.0 + 3.14159 + 2.71828: 10005.8
Rounding mode: ceiling
        Associative sum 10000.0 + 3.14159 + 2.71828: 10006.0
        Kahan       sum 10000.0 + 3.14159 + 2.71828: 10005.9

In no rounding mode are both answers correct. With "down" and "floor" rounding, the Kahan sum is too low (10005.8), but any other rounding makes it correct (10005.9). The Associative largest-to-smallest sum is never correct: "up" and "ceiling" rounding make it too high, while the rest make it low.

zkl

zkl floats are C doubles.

Translation of: go

<lang zkl>fcn kahanSum(numbers){

  sum:=c:=0.0;
  foreach x in (vm.arglist){
     y,t:=x - c, sum + y;
     c=(t - sum) - y;
     sum=t;
  }
  sum

} fcn epsilon{

   e:=1.0;
   while(1.0 + e!=1.0){ e/=2 }
   e

}</lang> <lang zkl>a,b,c,sum:=1.0,epsilon(),-b,a + b + c; sum :"%.20f".fmt(_).println("\tLeft associative. Delta from 1: ",1.0 - sum); kahanSum(a,b,c) :"%.20f".fmt(_).println("\tKahan summation"); b.println("\t\tEpsilon");</lang>

Output:

0.99999999999999988898	Left associative. Delta from 1: 1.11022e-16
1.00000000000000000000	Kahan summation
1.11022e-16		Epsilon