Advanced Search
All Posts : Applications, Examples and Libraries

add, floats, and Kahan sum

vv 12453 Maple
April 27 2018
3

add, floats, and Kahan sum

 

I found an intresting fact about the Maple command add for floating point values.
It seems that add in this case uses a summation algorithm in order to reduce the numerical error.
It is probably the Kahan summation algorithm (see wiki), but I wonder why this fact is not documented.

Here is a simple Maple procedure describing and implementing the algorithm.

 

 

> 

restart;
> 

Digits:=15;

15

 (1)
> 

KahanSum := proc(f::procedure, ab::range)  
local S,c,y,t, i;      # https://en.wikipedia.org/wiki/Kahan_summation_algorithm
S := 0.0;              # S = result (final sum: add(f(n), n=a..b))
c := 0.0;              # c = compensation for lost low-order bits.
for i from lhs(ab) to rhs(ab) do
    y := f(i) - c;     
    t := S + y;              
    c := (t - S) - y;        
    S := t;                  
od;                         
return S
end proc:

 

Now, a numerical example.

 

 

> 

f:= n ->  evalf(1/(n+1/n^3+1) - 1/(n+1+1/(n+1)^3+1));

proc (n) options operator, arrow; evalf(1/(n+1/n^3+1)-1/(n+2+1/(n+1)^3)) end proc

 (2)
> 

n := 50000;
K := KahanSum(f, 1..n);

50000

  

.333313334133301

 (3)
> 

A := add(f(k),k=1..n);

.333313334133302

 (4)
> 

s:=0.0:  for i to n do s:=s+f(i) od:
's' = s;

s = .333313334133413

 (5)
> 

exact:=( 1/3 - 1/(n+1+1/(n+1)^3+1) );

6250249999999900000/18751875067501050009

 (6)
> 

evalf( [errK = K-exact, errA = A-exact, err_for=s-exact] );

[errK = 0., errA = 0.1e-14, err_for = 0.112e-12]

 (7)
> 

evalf[20]( [errK = K-exact, errA = A-exact, err_for=s-exact] );

[errK = -0.33461e-15, errA = 0.66539e-15, err_for = 0.11166539e-12]

 (8)
> 

 


Download KahanSum.mw
evalf add

Please Wait...