Rosetta Code

Kahan summation: Difference between revisions

Page Discussion

Kahan summation (view source)

Revision as of 23:46, 6 March 2024

50,071 bytes added ,  2 months ago
Grouping BASIC dialects
m (→‎{{header|Perl 6}}: reduce the echo)
(Grouping BASIC dialects)
 
(47 intermediate revisions by 24 users not shown)
Line 4:
: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.
</small>
 
 
The [[wp:Kahan summation algorithm|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 &nbsp; '''Kahan summation algorithm''' &nbsp; is also known as the &nbsp; '''summation with carry algorithm'''.)
 
 
'''The task''' is to follow the previously linked Wikipedia articles algorithm and its [[wp:Kahan_summation_algorithm#Worked_example|worked example]] by:
Line 17 ⟶ 22:
'''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 <code>a, b, c</code> use the floating-point values <code>1.0, epsilon, -epsilon</code> respectively where epsilon is determined by its final value after the following:
<langsyntaxhighlight lang="python">epsilon = 1.0
while 1.0 + epsilon != 1.0:
epsilon = epsilon / 2.0</langsyntaxhighlight>
The above should ensure that <code>(a + b) + c != 1.0</code> 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.
 
Line 31 ⟶ 36:
* All examples should have constants chosen to clearly show the benefit of Kahan summing!
<br><br>
 
=={{header|11l}}==
{{trans|Python}}
 
<syntaxhighlight lang="11l">F kahansum(input)
V summ = 0.0s
V c = 0.0s
L(num) input
V y = num - c
V t = summ + y
c = (t - summ) - y
summ = t
R summ
 
V eps = 1.0s
L 1.0s + eps != 1.0s
eps = eps / 2.0s
 
print(‘Epsilon = ’eps)
print(‘(a + b) + c = #.7’.format((1.0s + eps) - eps))
print(‘Kahan sum = #.7’.format(kahansum([1.0s, eps, -eps])))</syntaxhighlight>
 
{{out}}
<pre>
Epsilon = 5.960464e-8
(a + b) + c = 0.9999999
Kahan sum = 1.0000000
</pre>
 
=={{header|Ada}}==
<syntaxhighlight lang="ada">with Ada.Text_Io; use Ada.Text_Io;
 
procedure Kahan is
 
type Kahan_Summer is
record
Sum : Float := 0.0;
C : Float := 0.0;
end record;
 
procedure Add (Acc : in out Kahan_Summer; Right : Float) is
Y : constant Float := Right - Acc.C;
T : constant Float := Acc.Sum + Y;
begin
Acc.C := (T - Acc.Sum) - Y;
acc.Sum := T;
end Add;
 
function Sum (Acc : Kahan_Summer) return Float
is (Acc.Sum);
 
function Epsilon return Float is
E : Float := 1.000;
begin
while 1.0 + E /= 1.0 loop
E := E / 2.0;
end loop;
return E;
end Epsilon;
 
A : constant Float := 1.000;
B : constant Float := Epsilon;
C : constant Float := -B;
D : constant Float := (A + B) + C;
K : Kahan_Summer;
 
package Float_Io is new Ada.Text_Io.Float_Io (Float);
use Float_Io;
begin
Add (K, A);
Add (K, B);
Add (K, C);
 
Default_Exp := 0;
Default_Aft := 12;
Put ("Epsilon : "); Put (B); New_Line;
Put ("(A + B) - C : "); Put (D); New_Line;
Put ("Kahan sum : "); Put (Sum (K)); New_Line;
end Kahan;</syntaxhighlight>
{{out}}
<pre>
Epsilon : 0.000000059605
(A + B) - C : 0.999999940395
Kahan sum : 1.000000000000
</pre>
 
=={{header|ALGOL 68}}==
Defines and uses a 6 digit float decimal type to solve the main "1000.0 + pi + e" task.
<langsyntaxhighlight 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 #
Line 191 ⟶ 281:
SKIP
)
</syntaxhighlight>
</lang>
{{out}}<pre>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</pre>
 
=={{header|Arturo}}==
 
<syntaxhighlight lang="arturo">kahanSum: function [inp][
c: 0
result: 0
loop inp 'val [
y: val - c
t: result + y
c: (t-result) - y
result: t
]
return result
]
 
eps: 1.0
while [1 <> 1 + eps]->
eps: eps / 2
 
a: 1.0
b: eps
c: neg eps
 
print ["Sum of 1.0, epsilon and -epsilon for epsilon:" eps]
 
print ["Is result equal to 1.0?"]
print ["- simple addition:" 1 = (a + b) + c]
print ["- using Kahan sum:" one? kahanSum @[a b c]]</syntaxhighlight>
 
{{out}}
 
<pre>Sum of 1.0, epsilon and -epsilon for epsilon: 1.110223024625157e-16
Is result equal to 1.0?
- simple addition: false
- using Kahan sum: true</pre>
 
=={{header|Asymptote}}==
{{trans|FreeBASIC}}
<syntaxhighlight lang="Asymptote">real a, b, c;
 
real KahanSum(real a, real b, real c) {
real sum = 0.0, y, t;
c = 0.0;
for (real i = 1; i <= a; ++i) {
y = i - c;
t = sum + y;
c = (t - sum) - y;
sum = t;
}
return sum;
}
 
real epsilon() {
real eps = 1;
while (1 + eps != 1) {
eps /= 2;
}
return eps;
}
 
a = 1.0;
b = epsilon();
c = -b;
 
real s = (a + b) + c;
real k = KahanSum(a, b, c);
real d = k - s;
write("Epsilon = " + string(b));
write("(a + b) + c = " + string(s));
write("Kahan sum = " + string(k));
write("Delta = " + string(d));</syntaxhighlight>
{{out}}
<pre>Epsilon = 1.11022302462516e-16
(a + b) + c = 1
Kahan sum = 1
Delta = 1.11022302462516e-16</pre>
 
=={{header|AWK}}==
<syntaxhighlight lang="awk">
# syntax: GAWK -f KAHAN_SUMMATION.AWK
# converted from C
BEGIN {
epsilon = 1
while (1 + epsilon != 1) {
epsilon /= 2
}
arr[1] = a = 1.0
arr[2] = b = epsilon
arr[3] = c = -b
printf("Epsilon = %18.16g\n",b)
printf("(a+b)+c = %18.16f\n",(a+b)+c)
printf("Kahan sum = %18.16f\n",kahan_sum(arr))
exit(0)
}
function kahan_sum(nums, c,i,sum,t,y) {
for (i=1; i<=length(nums); i++) {
y = nums[i] - c
t = sum + y
c = (t - sum) - y
sum = t
}
return(sum)
}
</syntaxhighlight>
{{out}}
<pre>
Epsilon = 1.110223024625157e-016
(a+b)+c = 0.9999999999999999
Kahan sum = 1.0000000000000000
</pre>
 
=={{header|BASIC}}==
==={{header|FreeBASIC}}===
<syntaxhighlight lang="vbnet">Dim Shared As Double a, b, c
 
Function KahanSum (a As Double, b As Double, c As Double) As Double
Dim As Double sum = 0.0, i, y, t
c = 0.0
For i = 1 To a
y = i - c
t = sum + y
c = (t - sum) - y
sum = t
Next i
Return sum
End Function
 
Function epsilon() As Double
Dim As Double eps = 1
While (1 + eps <> 1)
eps /= 2
Wend
Return eps
End Function
 
a = 1.0
b = epsilon()
c = -b
 
Dim As Double s = (a + b) + c
Dim As Double k = KahanSum(a, b, c)
Dim As Double d = k - s
Print "Epsilon ="; b
Print "(a + b) + c ="; s
Print "Kahan sum ="; k
Print "Delta ="; d
Sleep</syntaxhighlight>
{{out}}
<pre>
Epsilon = 1.110223024625157e-016
(a + b) + c = 0.9999999999999999
Kahan sum = 1
Delta = 1.110223024625157e-016
</pre>
 
==={{header|FutureBasic}}===
FB has proper decimal numbers supporting mantissas and exponents. But conversion to and from floating point numbers (or strings) makes it easier and more readable for this task to be completed with doubles as are many other examples here.
<syntaxhighlight lang="futurebasic">
_elements = 3
 
local fn Epsilon as double
double eps = 1.0
while ( 1.0 + eps != 1.0 )
eps = eps / 2.0
wend
end fn = eps
 
 
local fn KahanSum( nums(_elements) as double, count as long ) as double
double sum = 0.0
double c = 0.0
double t, y
long i
for i = 0 to count - 1
y = nums(i) - c
t = sum + y
c = (t - sum) - y
sum = t
next
end fn = sum
 
 
local fn DoKahan
double a = 1.0
double b = fn Epsilon
double c = -b
double fa[_elements]
fa(0) = a : fa(1) = b : fa(2) = c
printf @"Epsilon = %.9e", b
printF @"(a + b) + c = %.9e", (a + b) + c
printf @"Kahan sum = %.9e", fn KahanSum( fa(0), 3 )
printf @"Delta = %.9e", fn KahanSum( fa(0), 3 ) - ((a + b) + c)
end fn
 
fn DoKahan
 
HandleEvents</syntaxhighlight>
{{output}}
<pre>
Epsilon = 1.110223025e-16
(a + b) + c = 1.000000000e+00
Kahan sum = 1.000000000e+00
Delta = 1.110223025e-16
</pre>
 
==={{header|Gambas}}===
{{trans|FreeBASIC}}
<syntaxhighlight lang="vbnet">Public a As Float
Public b As Float
Public c As Float
 
Function KahanSum(a As Float, b As Float, c As Float) As Float
 
Dim sum As Float = 0.0
Dim i As Float, y As Float, t As Float
 
c = 0.0
For i = 1 To a
y = i - c
t = sum + y
c = (t - sum) - y
sum = t
Next
Return sum
 
End Function
 
Function epsilon() As Float
 
Dim eps As Float = 1
While (1 + eps <> 1)
eps /= 2
Wend
Return eps
 
End Function
 
Public Sub Main()
a = 1.0
b = epsilon()
c = -b
Dim s As Float = (a + b) + c
Dim k As Float = KahanSum(a, b, c)
Dim d As Float = k - s
Print "Epsilon = "; b
Print "(a + b) + c = "; s
Print "Kahan sum = "; k
Print "Delta = "; d
End</syntaxhighlight>
{{out}}
<pre>Epsilon = 1,11022302462516E-16
(a + b) + c = 1
Kahan sum = 1
Delta = 1,11022302462516E-16</pre>
 
==={{header|True BASIC}}===
<syntaxhighlight lang="qbasic">
FUNCTION epsilon
LET eps = 1
DO while (1+eps <> 1)
LET eps = eps/2
LOOP
LET epsilon = eps
END FUNCTION
 
FUNCTION kahansum(a, b, c)
LET sum = 0
LET c = 0
FOR i = 1 to a
LET y = i-c
LET t = sum+y
LET c = (t-sum)-y
LET sum = t
NEXT i
LET kahansum = sum
END FUNCTION
 
LET a = 1
LET b = epsilon
LET c = -b
LET s = (a+b)+c
LET k = kahansum(a, b, c)
LET d = k-s
PRINT "Epsilon ="; b
PRINT "(a + b) + c ="; s
PRINT "Kahan sum ="; k
PRINT "Delta ="; d
END</syntaxhighlight>
{{out}}
<pre>Epsilon = 1.110223e-16
(a + b) + c = 1.
Kahan sum = 1
Delta = 1.110223e-16</pre>
 
==={{header|Visual Basic .NET}}===
{{trans|C#}}
<syntaxhighlight lang="vbnet">Module Module1
 
Function KahanSum(ParamArray fa As Single()) As Single
Dim sum = 0.0F
Dim c = 0.0F
For Each f In fa
Dim y = f - c
Dim t = sum + y
c = (t - sum) - y
sum = t
Next
Return sum
End Function
 
Function Epsilon() As Single
Dim eps = 1.0F
While 1.0F + eps <> 1.0F
eps /= 2.0F
End While
Return eps
End Function
 
Sub Main()
Dim a = 1.0F
Dim b = Epsilon()
Dim c = -b
Console.WriteLine("Epsilon = {0}", b)
Console.WriteLine("(a + b) + c = {0}", (a + b) + c)
Console.WriteLine("Kahan sum = {0}", KahanSum(a, b, c))
End Sub
 
End Module</syntaxhighlight>
{{out}}
<pre>Epsilon = 1.110223E-16
(a + b) + c = 1
Kahan sum = 1</pre>
 
==={{header|Yabasic}}===
{{trans|FreeBASIC}}
<syntaxhighlight lang="vbnet">a = 1.0
b = epsilon()
c = -b
 
s = (a + b) + c
k = KahanSum(a, b, c)
d = k - s
print "Epsilon = ", b
print "(a + b) + c = ", s
print "Kahan sum = ", k
print "Delta = ", d
end
 
sub KahanSum (a, b, c)
sum = 0.0
c = 0.0
for i = 1 to a
y = i - c
t = sum + y
c = (t - sum) - y
sum = t
next i
return sum
end sub
 
sub epsilon()
eps = 1
while (1 + eps <> 1)
eps = eps / 2.0
wend
return eps
end sub</syntaxhighlight>
{{out}}
<pre>
Epsilon = 1.11022e-16
(a + b) + c = 1
Kahan sum = 1
Delta = 1.11022e-16
</pre>
 
=={{header|C}}==
{{trans|C++}}
<syntaxhighlight lang="c">#include <stdio.h>
#include <stdlib.h>
 
float epsilon() {
float eps = 1.0f;
while (1.0f + eps != 1.0f) eps /= 2.0f;
return eps;
}
 
float kahanSum(float *nums, int count) {
float sum = 0.0f;
float c = 0.0f;
float t, y;
int i;
for (i = 0; i < count; ++i) {
y = nums[i] - c;
t = sum + y;
c = (t - sum) - y;
sum = t;
}
return sum;
}
 
int main() {
float a = 1.0f;
float b = epsilon();
float c = -b;
float fa[3];
 
fa[0] = a;
fa[1] = b;
fa[2] = c;
 
printf("Epsilon = %0.12f\n", b);
printf("(a + b) + c = %0.12f\n", (a + b) + c);
printf("Kahan sum = %0.12f\n", kahanSum(fa, 3));
 
return 0;
}</syntaxhighlight>
{{out}}
<pre>Epsilon = 0.000000059605
(a + b) + c = 0.999999940395
Kahan sum = 1.000000000000</pre>
 
=={{header|C sharp|C#}}==
<syntaxhighlight lang="csharp">using System;
 
namespace KahanSummation {
class Program {
static float KahanSum(params float[] fa) {
float sum = 0.0f;
float c = 0.0f;
foreach (float f in fa) {
float y = f - c;
float t = sum + y;
c = (t - sum) - y;
sum = t;
}
 
return sum;
}
 
static float Epsilon() {
float eps = 1.0f;
while (1.0f + eps != 1.0f) eps /= 2.0f;
return eps;
}
 
static void Main(string[] args) {
float a = 1.0f;
float b = Epsilon();
float c = -b;
Console.WriteLine("Epsilon = {0}", b);
Console.WriteLine("(a + b) + c = {0}", (a + b) + c);
Console.WriteLine("Kahan sum = {0}", KahanSum(a, b, c));
}
}
}</syntaxhighlight>
{{out}}
<pre>Epsilon = 1.110223E-16
(a + b) + c = 1
Kahan sum = 1</pre>
 
=={{header|C++}}==
<syntaxhighlight lang="cpp">#include <iostream>
 
float epsilon() {
float eps = 1.0f;
while (1.0f + eps != 1.0f) eps /= 2.0f;
return eps;
}
 
float kahanSum(std::initializer_list<float> nums) {
float sum = 0.0f;
float c = 0.0f;
for (auto num : nums) {
float y = num - c;
float t = sum + y;
c = (t - sum) - y;
sum = t;
}
return sum;
}
 
int main() {
using namespace std;
 
float a = 1.f;
float b = epsilon();
float c = -b;
 
cout << "Epsilon = " << b << endl;
cout << "(a + b) + c = " << (a + b) + c << endl;
cout << "Kahan sum = " << kahanSum({ a, b, c }) << endl;
 
return 0;
}</syntaxhighlight>
{{out}}
<pre>Epsilon = 5.96046e-08
(a + b) + c = 1
Kahan sum = 1</pre>
 
=={{header|Crystal}}==
{{works with|crystal|0.31.1}}
{{trans|C++}}
 
<syntaxhighlight lang="crystal">def epsilon
eps = 1.0_f32
while 1.0_f32 + eps != 1.0_f32
eps /= 2.0_f32
end
 
eps
end
 
def kahan(nums)
sum = 0.0_f32
c = 0.0_f32
nums.each do |num|
y = num - c
t = sum + y
c = (t - sum) - y
sum = t
end
 
sum
end
 
a = 1.0_f32
b = epsilon
c = -b
 
puts "Epsilon = #{b}"
puts "Sum = #{a + b + c}"
puts "Kahan sum = #{kahan([a, b, c])}"</syntaxhighlight>
 
{{out}}
<pre>Epsilon = 5.9604645e-8
Sum = 0.99999994
Kahan sum = 1.0
</pre>
 
=={{header|D}}==
<syntaxhighlight lang="d">import std.stdio;
 
float kahanSum(float[] fa...) {
float sum = 0.0;
float c = 0.0;
foreach (f; fa) {
float y = f - c;
float t = sum + y;
c = (t - sum) - y;
sum = t;
}
return sum;
}
 
void main() {
float a = 1.0;
float b = float.epsilon;
float c = -b;
writefln("Epsilon = %0.8e", b);
writefln("(a + b) + c = %0.8f", ((a + b) + c));
writefln("Kahan sum = %0.8f", kahanSum(a, b, c));
}</syntaxhighlight>
{{out}}
<pre>Epsilon = 1.19209290e-07
(a + b) + c = 1.00000000
Kahan sum = 1.00000000</pre>
 
=={{header|Dart}}==
{{trans|C++}}
<syntaxhighlight lang="dart">double epsilon() {
double eps = 1.0;
while (1.0 + eps != 1.0) {
eps /= 2.0;
}
return eps;
}
 
double kahanSum(List<double> nums) {
double sum = 0.0;
double c = 0.0;
for (var num in nums) {
double y = num - c;
double t = sum + y;
c = (t - sum) - y;
sum = t;
}
return sum;
}
 
void main() {
double a = 1.0;
double b = epsilon();
double c = -b;
 
print("Epsilon = $b");
print("(a + b) + c = ${(a + b) + c}");
print("Kahan sum = ${kahanSum([a, b, c])}");
print("Delta = ${kahanSum([a, b, c]) - (a + b) + c}");
}</syntaxhighlight>
{{out}}
<pre>Epsilon = 1.1102230246251565e-16
(a + b) + c = 0.9999999999999999
Kahan sum = 1
Delta = -1.1102230246251565e-16</pre>
 
=={{header|EchoLisp}}==
No fixed point addition. EchoLisp floating point numbers are always stored as double precision floating point numbers, following the international IEEE 754 standard.
<langsyntaxhighlight 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 :
Line 236 ⟶ 936:
 
;;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.</langsyntaxhighlight>
 
=={{header|F sharp|F#}}==
{{incomplete|F sharp|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.}}
Solving the alternative task with an recursive algorithm and IEEE 754 double precision datatype.
<langsyntaxhighlight lang="fsharp">let ε =
let rec δ λ =
match λ+1.0 with
Line 266 ⟶ 966:
printfn "Sum: %e" (input |> List.sum)
printfn "Kahan: %e" (input |> Σ)
0</langsyntaxhighlight>
{{out}}
<pre>ε: 1.110223e-016
Sum: 1.000000e+000
Kahan: 1.000000e+000</pre>
 
=={{header|Factor}}==
Factor's floats use IEEE 754 double-precision format. Note the <code>kahan-sum</code> word already exists in Factor's standard library (in the <code>math.extras</code> vocabulary). See the implementation [https://docs.factorcode.org/content/word-kahan-sum%2Cmath.extras.html here].
{{works with|Factor|0.99 2019-10-06}}
<syntaxhighlight lang="factor">USING: io kernel literals math math.extras prettyprint sequences ;
 
CONSTANT: epsilon $[ 1.0 [ dup 1 + 1 number= ] [ 2 /f ] until ]
 
${ 1.0 epsilon dup neg }
[ "Left associative: " write sum . ]
[ "Kahan summation: " write kahan-sum . ] bi
"Epsilon: " write epsilon .</syntaxhighlight>
{{out}}
<pre>
Left associative: 0.9999999999999999
Kahan summation: 1.0
Epsilon: 1.110223024625157e-16
</pre>
 
=={{header|Fortran}}==
Line 277 ⟶ 995:
 
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.
<langsyntaxhighlight Fortranlang="fortran"> FUNCTION SUMC(A,N)
COMPENSATED SUMMATION. C WILL NOT STAY ZERO, DESPITE MATHEMATICS.
DIMENSION A(12345)
Line 298 ⟶ 1,016:
1 FORMAT (6HSUM = ,F12.1)
END
</syntaxhighlight>
</lang>
Alas, I no longer have access to an IBM1620 (or an emulator) whereby to elicit output. Despite being a language style over half a century old, 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 - or even base three on some Russian computers. Similarly, the system's standard word size may be 16, 18 (PDP 15), 32, 36, 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 4'''7''' bits. The 48'th bit was unused in arithmetic, including that of integers. Since a CHARACTER type was unavailable, text was placed up to six eight-bit (EBCDIC) codes to a word, and characters differing in the value of the topmost bit were indistinguishable. A special operator .IS. was introduced to supplement .EQ.
 
Line 305 ⟶ 1,023:
===Second phase: decimal precision via binary floating point===
Decimal fractions are nearly always recurring sequences in binary (i.e. except for powers of two such as three eights, etc) so at once there are approximations and these may combine unhelpfully. A value such as 3.14159 uses up nearly all of the precision of a 32-bit binary floating-point number, and is not represented exactly. This applies in turn to the various intermediate results as the summation proceeds, and one might be unlucky with the demonstration. The TRUNC6 approach alas doesn't work, as both summations come out as 10005·8, but with ROUND6, the desired results are presented...
<langsyntaxhighlight Fortranlang="fortran"> REAL FUNCTION TRUNC6(X) !Truncate X to six-digit decimal precision.
REAL X !The number.
REAL A !Its base ten log.
Line 370 ⟶ 1,088:
1 FORMAT (A1,"Sum = ",F12.1)
END
</syntaxhighlight>
</lang>
Output is
<pre>SSum = 10005.8
Line 397 ⟶ 1,115:
 
It is a small relief to see that 4*ATAN(1) gives a correctly-rounded value - at least with this system. Specifying a decimal value, even if with additional digits, may not yield the best result in binary. To find out exactly what binary values are being used, one needs a scheme as follows:
<langsyntaxhighlight Fortranlang="fortran"> MODULE PROBE
CONTAINS
CHARACTER*76 FUNCTION FP8DIGITS(X,BASE,W) !Full expansion of the value of X in BASE.
Line 544 ⟶ 1,262:
1 FORMAT (A1,"Sum = ",F12.1)
666 FORMAT (A8,":",A)
END</langsyntaxhighlight>
Because of the use of a function returning a CHARACTER value, and an array with a lower bound of zero, F90 features are needed, and to reduce the complaints over non-declaration of types of invoked functions, the MODULE protocol was convenient. Startlingly, if function SUMC was invoked with the wrong value of N (after removing a fourth element, Euler's constant, 0·5772156649015..., but forgetting to change the 4 back to 3) there was no complaint from the array bound-checking, and a fourth element was accessed! (It was zero) Changing the declaration from the old-style of A(N) to A(:) revived the bound checking, but at the loss of a documentation hint.
 
Line 607 ⟶ 1,325:
 
F90 offers interesting functions related to the floating-point number scheme, such as FRACTION(x) which returns the mantissa that for binary floating point will be in the range [0·5,1) so FRACTION(3·0) = 0·75, and EXPONENT(x), which will be for powers of two in a binary scheme. RADIX(x) will reveal the base and DIGITS(x) the number of digits of precision, while EPSILON(x) give the smallest value that can be distinguished from one. In these functions the parameter's value is irrelevant, it is present only to select between single and double precision, or perhaps quadruple precision. However, revealing their values in base ten is not entirely helpful, so using FP8DIGITS as before, proceed as follows:
<langsyntaxhighlight Fortranlang="fortran"> USE PROBE
REAL A(3)
REAL*4 X4E,X4L,X4S,X4C
Line 689 ⟶ 1,407:
665 FORMAT (A1,"Sum = ",F12.1)
666 FORMAT (A8,":",A)
END</langsyntaxhighlight>
Output, again with some minor editing. Because free-format output presents numbers in "scientific format" with a value in the units digit while "E"-type format does so with a zero units digit (or no units digit at all if cramped), the 1P prefix is used to cause a shift of one place and the exponent is adjusted accordingly. The same shift can be applied to "F"-type format but in that case there is no exponent to correct.
<pre>
Line 780 ⟶ 1,498:
=={{header|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:
<langsyntaxhighlight lang="go">package main
 
import "fmt"
Line 812 ⟶ 1,530:
fmt.Println("Kahan summation: ", kahan(a, b, c))
fmt.Println("Epsilon: ", b)
}</langsyntaxhighlight>
{{out}}
<pre>Left associative: 0.9999999999999999
Line 828 ⟶ 1,546:
 
Required code (for the current draft of this task):
<langsyntaxhighlight Jlang="j">epsilon=:3 :0''
epsilon=. 1.0
while. 1.0 ~: 1.0 + epsilon do.
Line 849 ⟶ 1,567:
sum=. t
end.
)</langsyntaxhighlight>
 
Required results:
 
<langsyntaxhighlight Jlang="j"> (a+b)+c
1
KahanSum a,b,c
1</langsyntaxhighlight>
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.
 
Line 862 ⟶ 1,580:
 
First, the value of epsilon:
<langsyntaxhighlight Jlang="j"> epsilon
5.68434e_14</langsyntaxhighlight>
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:
<langsyntaxhighlight Jlang="j">epsilon=:3 :0''
epsilon=. 1.0
while. 1.0 (~:!.0) 1.0 + epsilon do.
epsilon=. epsilon % 2.0
end.
)</langsyntaxhighlight>
Now we have a smaller value for this explicit value of epsilon:
<langsyntaxhighlight Jlang="j"> epsilon
1.11022e_16</langsyntaxhighlight>
With this new value for epsilon, let's try the problem again:
 
<langsyntaxhighlight Jlang="j"> (a+b)+c
1
KahanSum a,b,c
1</langsyntaxhighlight>
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:
 
<langsyntaxhighlight Jlang="j"> ":!.16 (a+b)+c
0.9999999999999999
":!.16 KahanSum a,b,c
1</langsyntaxhighlight>
Voila!
See also http://keiapl.info/anec/#fuzz
 
=={{header|Java}}==
{{trans|Kotlin}}
As is noted in the Kotlin implementation, the JVM does not provide the desired decimal type, so the alternate implementation is provided instead.
<langsyntaxhighlight Javalang="java">public class KahanSummation {
private static float kahanSum(float... fa) {
float sum = 0.0f;
Line 920 ⟶ 1,639:
System.out.println("Kahan sum = " + kahanSum(a, b, c));
}
}</langsyntaxhighlight>
{{out}}
<pre>Epsilon = 5.9604645E-8
(a + b) + c = 0.99999994
Kahan sum = 1.0</pre>
 
=={{header|Julia}}==
Julia can use its BigFloat data type to avoid floating point epsilon errors, as shown below.
<syntaxhighlight lang="julia">epsilon() = begin eps = 1.0; while 1.0 + eps != 1.0 eps /= 2.0 end; eps end
 
function kahansum(arr)
tot = temp = 0.0
for x in arr
y = x - temp
t = tot + y
temp = (t - tot) - y
tot = t
end
return tot
end
 
const a = 1.0
const ep = epsilon()
const b = -ep
const v = [a, ep, b]
 
println("Epsilon is $ep")
println("(a + ep) + -ep = ", (a + ep) + b)
println("Kahan sum is ", kahansum(v))
println("BigFloat sum is ", (BigFloat(a) + ep) + b)
</syntaxhighlight>{{out}}
<pre>
Epsilon is 1.1102230246251565e-16
(a + ep) + -ep = 0.9999999999999999
Kahan sum is 1.0
BigFloat sum is 1.0
</pre>
 
=={{header|Kotlin}}==
Kotlin does not have a 6 digit decimal type. The closest we have to it is the 4-byte floating point type, Float, which has a precision of 6 to 9 significant decimal digits - about 7 on average. So performing the alternative task:
<langsyntaxhighlight lang="scala">// version 1.1.2
 
fun kahanSum(vararg fa: Float): Float {
Line 955 ⟶ 1,706:
println("(a + b) + c = ${(a + b) + c}")
println("Kahan sum = ${kahanSum(a, b, c)}")
}</langsyntaxhighlight>
{{out}}
<pre>Epsilon = 5.9604645E-8
(a + b) + c = 0.99999994
Kahan sum = 1.0</pre>
 
=={{header|Lua}}==
{{trans|C}}
<syntaxhighlight lang="lua">function epsilon()
local eps = 1.0
while 1.0 + eps ~= 1.0 do
eps = eps / 2.0
end
return eps
end
 
function kahanSum(nums)
local sum = 0.0
local c = 0.0
for _,num in ipairs(nums) do
local y = num - c
local t = sum + y
c = (t - sum) - y
sum = t
end
return sum
end
 
-- main
local a = 1.0
local b = epsilon()
local c = -b
local fa = {a,b,c}
 
print(string.format("Epsilon = %0.24f", b))
print(string.format("(a + b) + c = %0.24f", (a+b)+c))
print(string.format("Kahan sum = %0.24f", kahanSum(fa)))</syntaxhighlight>
{{out}}
<pre>Epsilon = 0.000000000000000111022302
(a + b) + c = 0.999999999999999890000000
Kahan sum = 1.000000000000000000000000</pre>
 
=={{header|Modula-2}}==
<syntaxhighlight lang="modula2">MODULE KahanSummation;
{{Output?}}
<lang modula2>MODULE KahanSummation;
FROM RealStr IMPORT RealToStr;
FROM Terminal IMPORT WriteString,WriteLn,ReadChar;
Line 1,024 ⟶ 1,811:
 
ReadChar;
END KahanSummation.</langsyntaxhighlight>
{{out}}
<pre>Epsilon = 1.1102230E-16
(a + b) + c = 1.00000000
Kahan sum = 1.00000000</pre>
 
=={{header|Nim}}==
In version 1.4.x, Nim doesn’t provide decimal point arithmetic. So we use floating points with the alternative task.
 
The standard module <code>std/sums</code> contains the procedure <code>sumKbn</code> which uses the Kahan-Babuška-Neumaier algorithm, an improvement of Kahan algorithm. So, we compute the sum of 1.0, epsilon and -epsilon using a translation of Kahan algorithm as described in Wikipedia article and also using <code>sumKbn</code>. As expected, in this simple case there is no difference between these algorithms, but in other cases, results may differ.
 
<syntaxhighlight lang="nim">import std/sums
 
# "std/sums" proposes the "sumKbn" function which uses the
# Kahan-Babuška-Neumaier algorithm, an improvement of Kahan algorithm.
 
func kahanSum[T](input: openArray[T]): T =
var c = T(0)
for val in input:
let y = val - c
let t = result + y
c = (t - result) - y
result = t
 
template isOne[T](n: T): string =
if n == 1: "yes" else: "no"
 
var epsilon = 1.0
while 1 + epsilon != 1:
epsilon = epsilon / 2
 
let a = 1.0
let b = epsilon
let c = -epsilon
 
echo "Computing sum of 1.0, epsilon and -epsilon for epsilon = ", epsilon, '.'
echo "Is result equal to 1.0?"
echo "- simple addition: ", (a + b + c).isOne
echo "- using Kahan sum: ", kahanSum([a, b, c]).isOne
echo "- using stdlib: ", sumKbn([a, b, c]).isOne</syntaxhighlight>
 
{{out}}
<pre>Computing sum of 1.0, epsilon and -epsilon for epsilon = 1.110223024625157e-16.
Is result equal to 1.0?
- simple addition: no
- using Kahan sum: yes
- using stdlib: yes</pre>
 
=={{header|Objeck}}==
{{trans|C#}}
<syntaxhighlight lang="objeck">class KahanSum {
function : Main(args : String[]) ~ Nil {
a := 1.0f;
b := Epsilon();
c := -b;
 
abc := (a + b) + c;
fa := Float->New[3]; fa[0] := a; fa[1] := b; fa[2] := c;
sum := KahanSum(fa);
 
"Epsilon = {$b}"->PrintLine();
"(a + b) + c = {$abc}"->PrintLine();
"Kahan sum = {$sum}"->PrintLine();
}
 
function : KahanSum(fa : Float[]) ~ Float {
sum := 0.0f;
c := 0.0f;
each(i : fa) {
f := fa[i];
y := f - c;
t := sum + y;
c := (t - sum) - y;
sum := t;
};
 
return sum;
}
 
function : Epsilon() ~ Float {
eps := 1.0f;
while(1.0f + eps <> 1.0f) { eps /= 2.0f; };
return eps;
}
}</syntaxhighlight>
 
{{output}}
<pre>
Epsilon = 1.11022e-16
(a + b) + c = 1
Kahan sum = 1
</pre>
 
=={{header|PARI/GP}}==
Line 1,030 ⟶ 1,908:
 
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.)
<langsyntaxhighlight 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])</langsyntaxhighlight>
{{out}}
<pre>%1 = 1.000000000000000000</pre>
 
=={{header|Perl 6}}==
{{trans|PythonRaku}}
<syntaxhighlight lang="perl">use strict;
Perl&nbsp;6 does not offer a fixed precision decimal. It ''does'' have IEEE 754 floating point numbers so let's try implementing the arbitrary precision option as shown in Python. Need to explicitly specify scientific notation numbers to force floating point Nums.
use warnings;
use feature 'say';
 
sub kahan {
<lang perl6>constant ε = (1e0, */2e0 … *+1e0==1e0)[*-1];
my(@nums) = @_;
 
sub kahan (*@nums) {
my $summ = my $c = 0e0;
for @nums ->my $num (@nums) {
my $y = $num - $c;
my $t = $summ + $y;
Line 1,055 ⟶ 1,934:
}
 
my $eps = 1;
say 'Epsilon: ', ε;
do { $eps /= 2 } until 1e0 == 1e0 + $eps;
 
say 'Simple sumEpsilon: ', ((1e0 + ε)' - ε).fmt: "%.16f"$eps;
say 'Simple sum: ' . sprintf "%.16f", ((1e0 + $eps) - $eps);
say 'Kahan sum: ' . sprintf "%.16f", kahan(1e0, $eps, -$eps);</syntaxhighlight>
{{out}}
<pre>Epsilon: 1.11022302462516e-16
Simple sum: 0.9999999999999999
Kahan sum: 1.0000000000000000</pre>
 
=={{header|Phix}}==
say 'Kahan sum: ', kahan(1e0, ε, -ε).fmt: "%.16f";</lang>
Phix does not have fixed precision decimals, just the usual IEEE 754 floating point numbers.<br>
Using enforced rounding. From the manual: "For the %f format, precision specifies how many digits follow <br>
the decimal point character, whereas for %g, the precision specifies how many significant digits to print",<br>
hence we can use printf(%g), followed immediately by scanf(), to emulate limited precision decimals.
<!--<syntaxhighlight lang="phix">(phixonline)-->
<span style="color: #008080;">function</span> <span style="color: #000000;">round6</span><span style="color: #0000FF;">(</span><span style="color: #004080;">atom</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">)</span>
<span style="color: #000080;font-style:italic;">-- NB: only suitable for this very specific example</span>
<span style="color: #004080;">string</span> <span style="color: #000000;">s</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">trim</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">sprintf</span><span style="color: #0000FF;">(</span><span style="color: #008000;">"%6g"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">a</span><span style="color: #0000FF;">))</span>
<span style="color: #0000FF;">{{</span><span style="color: #000000;">a</span><span style="color: #0000FF;">}}</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">scanf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">s</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"%f"</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">a</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">roundsum</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">s</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">s</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
<span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">round6</span><span style="color: #0000FF;">(</span><span style="color: #000000;">res</span><span style="color: #0000FF;">+</span><span style="color: #000000;">s</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">])</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">res</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">kahanroundsum</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">s</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">tot</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">c</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">s</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">x</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">s</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">],</span>
<span style="color: #000000;">y</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">round6</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;">-</span><span style="color: #000000;">c</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">t</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">round6</span><span style="color: #0000FF;">(</span><span style="color: #000000;">tot</span><span style="color: #0000FF;">+</span><span style="color: #000000;">y</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">c</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">round6</span><span style="color: #0000FF;">(</span><span style="color: #000000;">round6</span><span style="color: #0000FF;">(</span><span style="color: #000000;">t</span><span style="color: #0000FF;">-</span><span style="color: #000000;">tot</span><span style="color: #0000FF;">)-</span><span style="color: #000000;">y</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">tot</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">t</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">tot</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">procedure</span> <span style="color: #000000;">main</span><span style="color: #0000FF;">()</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">s</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">10000.0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">3.14159</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">2.71828</span><span style="color: #0000FF;">}</span>
<span style="color: #004080;">string</span> <span style="color: #000000;">fmt</span> <span style="color: #0000FF;">=</span> <span style="color: #008000;">"%5.1f\n"</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">fmt</span><span style="color: #0000FF;">,</span><span style="color: #000000;">roundsum</span><span style="color: #0000FF;">(</span><span style="color: #000000;">s</span><span style="color: #0000FF;">))</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">fmt</span><span style="color: #0000FF;">,</span><span style="color: #000000;">kahanroundsum</span><span style="color: #0000FF;">(</span><span style="color: #000000;">s</span><span style="color: #0000FF;">))</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span>
<span style="color: #000000;">main</span><span style="color: #0000FF;">()</span>
<!--</syntaxhighlight>-->
{{out}}
<pre>
<pre>Epsilon: 1.1102230246251565e-16
10005.8
10005.9
</pre>
Alternative task using floats and the explicitly calculated epsilon:
<!--<syntaxhighlight lang="phix">(phixonline)-->
<span style="color: #008080;">function</span> <span style="color: #000000;">kahansum</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">s</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">tot</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">c</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">s</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">x</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">s</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">],</span>
<span style="color: #000000;">y</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">x</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">c</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">t</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">tot</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">y</span>
<span style="color: #000000;">c</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">t</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">tot</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">y</span>
<span style="color: #000000;">tot</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">t</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">tot</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">epsilon</span><span style="color: #0000FF;">()</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">e</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span>
<span style="color: #008080;">while</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">+</span><span style="color: #000000;">e</span><span style="color: #0000FF;">!=</span><span style="color: #000000;">1</span> <span style="color: #008080;">do</span>
<span style="color: #000000;">e</span><span style="color: #0000FF;">/=</span><span style="color: #000000;">2</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">while</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">e</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">procedure</span> <span style="color: #000000;">main2</span><span style="color: #0000FF;">()</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">a</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1.0</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">b</span> <span style="color: #0000FF;">:=</span> <span style="color: #000000;">epsilon</span><span style="color: #0000FF;">(),</span>
<span style="color: #000000;">c</span> <span style="color: #0000FF;">:=</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">b</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">s</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,</span><span style="color: #000000;">b</span><span style="color: #0000FF;">,</span><span style="color: #000000;">c</span><span style="color: #0000FF;">}</span>
<span style="color: #004080;">string</span> <span style="color: #000000;">fmt</span> <span style="color: #0000FF;">=</span> <span style="color: #008080;">iff</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">machine_bits</span><span style="color: #0000FF;">()=</span><span style="color: #000000;">32</span><span style="color: #0000FF;">?</span><span style="color: #008000;">"%.16f\n"</span><span style="color: #0000FF;">:</span><span style="color: #008000;">"%.19f\n"</span><span style="color: #0000FF;">)</span>
<span style="color: #0000FF;">?{</span><span style="color: #008000;">"Epsilon:"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">b</span><span style="color: #0000FF;">}</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"Simple sum: "</span><span style="color: #0000FF;">&</span><span style="color: #000000;">fmt</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">sum</span><span style="color: #0000FF;">(</span><span style="color: #000000;">s</span><span style="color: #0000FF;">))</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"Kahan sum : "</span><span style="color: #0000FF;">&</span><span style="color: #000000;">fmt</span><span style="color: #0000FF;">,</span><span style="color: #000000;">kahansum</span><span style="color: #0000FF;">(</span><span style="color: #000000;">s</span><span style="color: #0000FF;">))</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span>
<span style="color: #000000;">main2</span><span style="color: #0000FF;">()</span>
<!--</syntaxhighlight>-->
{{out}}
32-bit
<pre>
{"Epsilon:",1.110223024e-16}
Simple sum: 0.9999999999999998
Kahan sum : 1.0000000000000000
</pre>
browser (also 32-bit but obviously some kind of subtly different rounding)
<pre>
{"Epsilon:",1.110223025e-16}
Simple sum: 0.9999999999999999
Kahan sum: : 1.0000000000000000
</pre>
64-bit
<pre>
{"Epsilon:",5.421010862e-20}
Simple sum: 0.9999999999999999999
Kahan sum : 1.0000000000000000000
</pre>
Above and beyond:
<!--<syntaxhighlight lang="phix">(phixonline)-->
<span style="color: #008080;">function</span> <span style="color: #000000;">NeumaierSum</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">s</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">s</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">],</span>
<span style="color: #000000;">c</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0.0</span> <span style="color: #000080;font-style:italic;">-- A running compensation for lost low-order bits.</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">2</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">s</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">t</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">res</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">s</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span>
<span style="color: #008080;">if</span> <span style="color: #7060A8;">abs</span><span style="color: #0000FF;">(</span><span style="color: #000000;">res</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">>=</span> <span style="color: #7060A8;">abs</span><span style="color: #0000FF;">(</span><span style="color: #000000;">s</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">])</span> <span style="color: #008080;">then</span>
<span style="color: #000000;">c</span> <span style="color: #0000FF;">+=</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">res</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">t</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">s</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> <span style="color: #000080;font-style:italic;">-- If res is bigger, low-order digits of s[i] are lost.</span>
<span style="color: #008080;">else</span>
<span style="color: #000000;">c</span> <span style="color: #0000FF;">+=</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">s</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">t</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">res</span> <span style="color: #000080;font-style:italic;">-- Else low-order digits of res are lost</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">t</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">res</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">c</span> <span style="color: #000080;font-style:italic;">-- Correction only applied once at the very end</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">kahansum</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">s</span><span style="color: #0000FF;">)</span> <span style="color: #000080;font-style:italic;">-- (copy of above)</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">tot</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">c</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">s</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">x</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">s</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">],</span>
<span style="color: #000000;">y</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">x</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">c</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">t</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">tot</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">y</span>
<span style="color: #000000;">c</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">t</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">tot</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">y</span>
<span style="color: #000000;">tot</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">t</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">tot</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">procedure</span> <span style="color: #000000;">main3</span><span style="color: #0000FF;">()</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">s</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">1e300</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">1e300</span><span style="color: #0000FF;">}</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"Simple sum: %d\n"</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">sum</span><span style="color: #0000FF;">(</span><span style="color: #000000;">s</span><span style="color: #0000FF;">))</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"Kahan sum : %d\n"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">kahansum</span><span style="color: #0000FF;">(</span><span style="color: #000000;">s</span><span style="color: #0000FF;">))</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"Neumaier sum : %d\n"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">NeumaierSum</span><span style="color: #0000FF;">(</span><span style="color: #000000;">s</span><span style="color: #0000FF;">))</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span>
<span style="color: #000000;">main3</span><span style="color: #0000FF;">()</span>
<!--</syntaxhighlight>-->
{{out}}
<pre>
Simple sum: 0
Kahan sum : 0
Neumaier sum : 1
</pre>
 
=={{header|PHP}}==
'''Script'''
<langsyntaxhighlight lang="php"><?php
// Main task
 
Line 1,137 ⟶ 2,159:
$diff = (($a + $b) + $c) - kahansum($input);
echo $diff.PHP_EOL;
</syntaxhighlight>
</lang>
'''Script output'''
<langsyntaxhighlight lang="bash">$ php kahansum.php
Main task - Left to right summation: 10005.9
Main task - Kahan summation: 10005.9
Line 1,149 ⟶ 2,171:
bool(false)
Difference between the operations: 1.1102230246252E-16
</syntaxhighlight>
</lang>
 
=={{header|Python}}==
===Python: Decimal===
<langsyntaxhighlight lang="python">>>> from decimal import *
>>>
>>> getcontext().prec = 6
Line 1,189 ⟶ 2,211:
>>> (a + b) + c
Decimal('10005.85987')
>>> </langsyntaxhighlight>
 
===Python: Floats===
This was written as proof that the floating-point sub-task could work and is better off displayed, so...
<langsyntaxhighlight lang="python">>>> eps = 1.0
>>> while 1.0 + eps != 1.0:
eps = eps / 2.0
Line 1,210 ⟶ 2,232:
>>> 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)
>>> </langsyntaxhighlight>
Note that the kahansum function from the decimal example is [[wp:Duck typing|duck typed]] and adjusts to work with the number type used in its argument list.
 
Line 1,217 ⟶ 2,239:
===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.
<langsyntaxhighlight lang="python">>>> from decimal import localcontext, Decimal
>>>
>>> with localcontext() as ctx:
Line 1,232 ⟶ 2,254:
Simple sum is: 0.9999999999999999999999999999
Kahan sum is: 1.000000000000000000000000000
>>> </langsyntaxhighlight>
 
=={{header|R}}==
Line 1,239 ⟶ 2,261:
Therefore, trying to use the main rules will not work as we can see in the code below
 
<langsyntaxhighlight lang="r"># an implementation of the Kahan summation algorithm
kahansum <- function(x) {
ks <- 0
Line 1,269 ⟶ 2,291:
kahansum(input)
# [1] 10005.86
</syntaxhighlight>
</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'''
<langsyntaxhighlight lang="r">epsilon = 1.0
while ((1.0 + epsilon) != 1.0) {
epsilon = epsilon / 2.0
Line 1,279 ⟶ 2,301:
epsilon
# [1] 1.11022e-16
</syntaxhighlight>
</lang>
and use it to test the left-to-right summation and the Kahan function
<langsyntaxhighlight lang="r">a <- 1.0
b <- epsilon
c <- -epsilon
Line 1,292 ⟶ 2,314:
kahansum(c(a, b, c))
# [1] 1
</syntaxhighlight>
</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
<langsyntaxhighlight lang="r">(a + b) + c == kahansum(c(a, b, c))
# FALSE
 
Line 1,300 ⟶ 2,322:
((a + b) + c) - kahansum(c(a, b, c))
# [1] -1.110223e-16
</syntaxhighlight>
</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
<langsyntaxhighlight lang="r"># machine epsilon
mepsilon <- function(divisor) {
guess <- 1.0
Line 1,318 ⟶ 2,340:
# 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
</syntaxhighlight>
</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
<langsyntaxhighlight lang="r">library(txtplot)
 
txtboxplot(epsilons)
Line 1,348 ⟶ 2,370:
# +--+---------+----------+----------+---------+----------+----------+
# 0 2e-17 4e-17 6e-17 8e-17 1e-16
</syntaxhighlight>
</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.
<langsyntaxhighlight lang="r">epsilon <- median(epsilons)
epsilon
# [1] 5.939024e-18
Line 1,370 ⟶ 2,392:
(a + b) + c == kahansum(c(a, b, c))
# TRUE
</syntaxhighlight>
</lang>
 
=={{header|Racket}}==
Line 1,378 ⟶ 2,400:
 
Finally we try the alternative task version for single precision float point numbers.
<langsyntaxhighlight Racketlang="racket">#lang racket
 
(define (sum/kahan . args)
Line 1,394 ⟶ 2,416:
(displayln "Double presition flonum")
(+ 1000000.0 3.14159 2.71828)
(sum/kahan 1000000.0 3.14159 2.71828)</langsyntaxhighlight>
{{out}}
<pre>Single presition flonum
Line 1,404 ⟶ 2,426:
 
Alternative task version
<langsyntaxhighlight Racketlang="racket">(define epsilon.f0
(let loop ([epsilon 1.f0])
(if (= (+ 1.f0 epsilon) 1.f0)
Line 1,414 ⟶ 2,436:
 
(+ 1.f0 epsilon.f0 (- epsilon.f0))
(sum/kahan 1.f0 epsilon.f0 (- epsilon.f0))</langsyntaxhighlight>
{{out}}
<pre>Alternative task, single precision flonum
Line 1,420 ⟶ 2,442:
0.99999994f0
1.0f0</pre>
 
=={{header|Raku}}==
(formerly Perl 6)
{{trans|Python}}
Raku does not offer a fixed precision decimal. It ''does'' have IEEE 754 floating point numbers so let's try implementing the floating point option as shown in Python. Need to explicitly specify scientific notation numbers to force floating point Nums.
 
<syntaxhighlight lang="raku" line>constant ε = (1e0, */2e0 … *+1e0==1e0)[*-1];
 
sub kahan (*@nums) {
my $summ = my $c = 0e0;
for @nums -> $num {
my $y = $num - $c;
my $t = $summ + $y;
$c = ($t - $summ) - $y;
$summ = $t;
}
$summ
}
 
say 'Epsilon: ', ε;
 
say 'Simple sum: ', ((1e0 + ε) - ε).fmt: "%.16f";
 
say 'Kahan sum: ', kahan(1e0, ε, -ε).fmt: "%.16f";</syntaxhighlight>
{{out}}
<pre>Epsilon: 1.1102230246251565e-16
Simple sum: 0.9999999999999999
Kahan sum: 1.0000000000000000
</pre>
 
=={{header|REXX}}==
Line 1,425 ⟶ 2,476:
<br>of &nbsp; '''30''' &nbsp; (decimal digits) was chosen. &nbsp; The default precision for REXX is &nbsp; '''9''' &nbsp; decimal digits.
===vanilla version===
<langsyntaxhighlight lang="rexx">/*REXX program demonstrates simple addition versus using Kahan summation algorithm. */
numeric digits 6 /*use six decimal digits for precision.*/
call show 10000.0, 3.14169, 2.71828 /*invoke SHOW to sum & display numbers.*/
numeric digits 30 /*from now on, use 30 decimal digits.*/
epsilon= 1
do while 1+epsilon \= 1 /*keep looping 'til we can't add unity.*/
epsilon= epsilon / 2 /*halve the value of epsilon variable.*/
end /*while*/
say /*display a blank line before the fence*/
Line 1,438 ⟶ 2,489:
exit /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
kahan: procedure; $= 0; c= 0; do j=1 for arg() /*perform for each argument. */
y=arg(j) - c y= arg(j) - c /*subtract C from argument.*/
t=$ + y t= $ + y /*use a temporary sum (T). */
c=t - $ - y c= t - $ - y /*compute the value of C. */
$=t $= t /*redefine the sum ($). */
end /*j*/
return $
/*──────────────────────────────────────────────────────────────────────────────────────*/
show: procedure; parse arg a,b,c /*obtain the arguments. */
say 'decimal digits =' digits() /*show number of decimal digs*/
say ' a = ' left(''"", a>=0) a a /*display A justified. */
say ' b = ' left(''"", b>=0) b b /* " B " */
say ' c = ' left(''"", c>=0) c c /* " C " */
say 'simple summation of a,b,c = ' a+ b+c c /*compute simple summation. */
say 'Kahan summation of a,b,c = ' kahan(a, b, c) /*sum via Kahan summation. */
return</langsyntaxhighlight>
{{out|output|text=&nbsp; using PC/REXX and also Personal REXX:}}
<pre>decimal digits = 6
Line 1,487 ⟶ 2,538:
Kahan summation of a,b,c = 1.00000000000000000000000000000
</pre>
{{out|output|text=&nbsp; using Regina version 3.3 and earlier:}}
<pre>
decimal digits = 6
Line 1,504 ⟶ 2,555:
simple summation of a,b,c = 0.999999999999999999999999999997
Kahan summation of a,b,c = 1.00000000000000000000000000000</pre>
{{out|output|text=&nbsp; using Regina version 3.4 and later &nbsp; (latest at this time is version 3.9.1):}}
<pre>decimal digits = 6
a = 10000.0
Line 1,522 ⟶ 2,573:
 
===tweaked version===
The following tweaked REXX version causes Regina (version 3.4 and later) to work properly.
<langsyntaxhighlight lang="rexx">/*REXX program demonstrates simple addition versus using Kahan summation algorithm. */
numeric digits 6 /*use six decimal digits for precision.*/
call show 10000.0, 3.14169, 2.71828 /*invoke SHOW to sum & display numbers.*/
numeric digits 30 /*from now on, use 30 decimal digits.*/
epsilon= 1
do while 1+epsilon \= 1 /*keep looping 'til we can't add unity.*/
epsilon= epsilon / 2 /*halve the value of epsilon variable.*/
end /*while*/
say /*display a blank line before the fence*/
Line 1,538 ⟶ 2,589:
exit /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
kahan: procedure; $= 0; c= 0; do j=1 for arg() /*perform for each argument. */
y=arg(j) - c y= arg(j) - c /*subtract C from argument.*/
t=$ + y t= $ + y /*use a temporary sum (T). */
c=t - $ - y c= t - $ - y /*compute the value of C. */
$=t $= t /*redefine the sum ($). */
end /*j*/
return $
/*──────────────────────────────────────────────────────────────────────────────────────*/
show: procedure; parse arg a,b,c /*obtain the arguments. */
say 'decimal digits =' digits() /*show number of decimal digs*/
say ' a = ' left(''"", a>=0) a a /*display A justified. */
say ' b = ' left(''"", b>=0) b b /* " B " */
say ' c = ' left(''"", c>=0) c c /* " C " */
say 'simple summation of a,b,c = ' a+ b+c c /*compute simple summation. */
say 'Kahan summation of a,b,c = ' kahan(a, b, c) /*sum via Kahan summation. */
return</langsyntaxhighlight>
{{out|output|text=&nbsp; using Regina version 3.4 and also for later versions:}}
<pre>decimal digits = 6
a = 10000.0
Line 1,569 ⟶ 2,620:
c = -3.15544362088404722164691426133E-30
simple summation of a,b,c = 1.0000000000000000000000000000001
Kahan summation of a,b,c = 1.0000000000000000000000000000000</pre>
</pre>
'''output''' for ooRexx:
<pre>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
 
=={{header|RPL}}==
ªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªª
{| class="wikitable" ≪
! RPL code
! Comment
|-
|
≪ → input
≪ 0 0
1 input SIZE '''FOR''' j
input j GET SWAP -
DUP2 +
DUP 4 ROLL - ROT -
'''NEXT''' DROP
≫ ≫ '<span style="color:blue">∑KAHAN</span>' STO
≪ .1
'''WHILE''' 1 OVER + 1 ≠ '''REPEAT''' 10 / '''END'''
→ eps
≪ 1 eps + eps NEG +
1 eps eps NEG →LIST <span style="color:blue">∑KAHAN</span>
≫ ≫ '<span style="color:blue">TASK</span>' STO
|
<span style="color:blue">∑KAHAN</span> ''( { a b .. n } → a+b+..+n ) ''
var sum = 0.0, c = 0.0
for i = 1 to input.length do
var y = input[i] - c
var t = sum + y
c = (t - sum) - y
sum = t
next i ; return sum
Calculate epsilon
Display 1 + epsilon - epsilon
Display KAHAN( { 1 epsilon -epsilon })
|}
{{out}}
<pre>
2: 0.999999999999
1: 1
</pre>
 
=={{header|Ruby}}==
decimal digits = 32
Arrays provide a "sum" method, which uses a Kahan-Babuska balancing compensated summation algorithm, according to the C-source.
a = 1.0
<syntaxhighlight lang="ruby">epsilon = 1.0
b = 0.00000000000000000000000000000315544362088404722164691426133
epsilon /= 2 until 1.0 + epsilon == 1.0
c = -0.00000000000000000000000000000315544362088404722164691426133
 
simple summation of a,b,c = 1.0000000000000000000000000000001
a = 1.0
Kahan summation of a,b,c = 1.0000000000000000000000000000000
b = epsilon
c = -b
 
puts "epsilon : #{epsilon}"
puts "(a+b)+c : #{(a+b)+c}"
puts "[a,b,c].sum: #{[a,b,c].sum}"
</syntaxhighlight>
{{out}}
<pre>epsilon : 1.1102230246251565e-16
(a+b)+c : 0.9999999999999999
[a,b,c].sum: 1.0
</pre>
 
=={{header|Rust}}==
<syntaxhighlight lang="rust">
extern crate num;
extern crate permutohedron;
 
use num::Float;
use permutohedron::Heap;
use std::f32;
 
fn find_max(lst: &[f32]) -> Option<f32> {
if lst.is_empty() {
return None;
}
let max = lst.iter().fold(f32::NEG_INFINITY, |a, &b| Float::max(a, b));
Some(max)
}
 
fn with_bits(val: f32, digits: usize) -> f32 {
let num = format!("{:.*}", digits, val);
num.parse::<f32>().unwrap()
}
 
fn kahan_sum(lst: &[f32]) -> Option<f32> {
let mut sum = 0.0f32;
let mut c = 0.0f32;
for i in lst {
let y = *i - c;
let t = sum + y;
c = (t - sum) - y;
sum = t;
}
Some(with_bits(sum, 1))
}
 
fn all_sums(vec: &mut [f32]) -> Vec<f32> {
let mut res = Vec::new();
let mut perms = Heap::new(vec);
loop {
let v = perms.next();
match v {
Some(v) => {
let mut sum = 0.0f32;
for e in &v {
sum += with_bits(*e, 1);
}
res.push(with_bits(sum, 1));
}
None => break,
}
}
res
}
 
#[allow(clippy::approx_constant)]
fn main() {
let v = vec![10_000f32, 3.14159, 2.71828];
let sums = all_sums(&mut v.clone());
let res = kahan_sum(&v).unwrap();
let max = find_max(&sums[..]).unwrap();
println!("max: {} res: {}", max, res);
}
 
</syntaxhighlight>
{{out}}
<pre>
max: 10005.8 res: 10005.9
</pre>
 
ªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªªª</pre>
=={{header|Scala}}==
===IEEE 754 Single precision 32-bit (JavaScript defaults to Double precision.)===
<langsyntaxhighlight Scalalang="scala">import scala.annotation.tailrec
 
object KahanSummation extends App {
Line 1,618 ⟶ 2,786:
println("Kahan sum = " + kahanSum(a, ε, -ε))
}
}</langsyntaxhighlight>
{{Out}}See it running in your browser by [https://scalafiddle.io/sf/lWI0SJC/0 ScalaFiddle (JavaScript, non JVM)] or by [https://scastie.scala-lang.org/CI7dNAKQSlupivKAnYcGJg Scastie (remote JVM)].
Note: JVM float is IEEE-754 32 bit floating point while JavaScript always default to the IEEE 754 standard 64-bit.
 
=={{header|Tcl}}==
===Tcl: Floats===
First, using native floating point we see the same epsilon value as other languages using float64:
 
<langsyntaxhighlight Tcllang="tcl"># make {+ - * /} etc available as commands, for easier expressions
namespace path ::tcl::mathop
 
Line 1,653 ⟶ 2,822:
puts "\tEpsilon is: [set e [epsilon]]"
puts "\tAssociative sum: [expr {1.0 + $e - $e}]"
puts "\tKahan sum: [kahansum 1.0 $e -$e]"</langsyntaxhighlight>
<pre>Epsilon is: 1.1102230246251565e-16
Associative sum: 0.9999999999999999
Line 1,663 ⟶ 2,832:
 
The last stanza exercises the decimal package's different rounding modes, to see what happens there:
<langsyntaxhighlight Tcllang="tcl">package require math::decimal
namespace path ::math::decimal
 
Line 1,697 ⟶ 2,866:
puts "\tAssociative sum $a + $b + $c: [asum $a $b $c]"
puts "\tKahan sum $a + $b + $c: [kahansum $a $b $c]"
}</langsyntaxhighlight>
The results are a little surprising:
{{out}}
Line 1,724 ⟶ 2,893:
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.
 
=={{header|V (Vlang)}}==
{{trans|go}}
Vlang has no floating point decimal type. Its floating point types are f64 and f32, implementations of IEEE 754 binary64 and binary32 formats. Alternative task:
<syntaxhighlight lang="v (vlang)">fn kahan(s ...f64) f64 {
mut tot, mut c := 0.0, 0.0
for x in s {
y := x - c
t := tot + y
c = (t - tot) - y
tot = t
}
return tot
}
fn epsilon() f64 {
mut e := 1.0
for 1+e != 1 {
e /= 2
}
return e
}
fn main() {
a := 1.0
b := epsilon()
c := -b
println("Left associative: ${a+b+c}")
println("Kahan summation: ${kahan(a, b, c)}")
println("Epsilon: $b")
}</syntaxhighlight>
{{out}}
<pre>Left associative: 0.9999999999999999
Kahan summation: 1
Epsilon: 1.1102230246251565e-16
</pre>
With float defined as float32,
<pre>
Left associative: 0.99999994
Kahan summation: 1
Epsilon: 5.9604645e-08</pre>
 
=={{header|Wren}}==
Wren's only 'native' number type (Num) is double-precision floating point and so the alternative task is performed. Although it appears that there is no difference between the left associative sum and Kahan summation, there is in fact a difference of epsilon but the accuracy of System.print (14 significant figures) is insufficient to indicate this directly.
<syntaxhighlight lang="wren">var kahanSum = Fn.new { |a|
var sum = 0
var c = 0
for (f in a) {
var y = f - c
var t = sum + y
c = (t - sum) - y
sum = t
}
return sum
}
 
var epsilon = Fn.new {
var eps = 1
while (1 + eps != 1) eps = eps/2
return eps
}
 
var a = 1
var b = epsilon.call()
var c = -b
 
var s = (a + b) + c
var k = kahanSum.call([a, b, c])
var d = k - s
System.print("Epsilon = %(b)")
System.print("(a + b) + c = %(s)")
System.print("Kahan sum = %(k)")
System.print("Delta = %(d)")</syntaxhighlight>
 
{{out}}
<pre>
Epsilon = 1.1102230246252e-16
(a + b) + c = 1
Kahan sum = 1
Delta = 1.1102230246252e-16
</pre>
 
=={{header|XPL0}}==
{{trans|C}}
XPL0's only 'native' floating point type (real) is double-precision, so
the alternative task is performed. The normal output shows no difference
because of rounding, so the hex output is used to show it.
<syntaxhighlight lang "XPL0">include xpllib; \for Print
 
func real KahanSum(Nums, Count);
real Nums; int Count;
real Sum, C, T, Y;
int I;
[Sum:= 0.0;
C:= 0.0;
for I:= 0 to Count-1 do
[Y:= Nums(I) - C;
T:= Sum + Y;
C:= (T - Sum) - Y;
Sum:= T;
];
return Sum;
];
 
func real Epsilon;
real Eps;
[Eps:= 1.0;
while 1.0 + Eps # 1.0 do
Eps:= Eps/2.0;
return Eps;
];
 
real A, B, C, FA(3), D, K;
[A:= 1.0; B:= Epsilon; C:= -B;
FA(0):= A; FA(1):= B; FA(2):= C;
Print("Epsilon = %0.16f\n", B);
D:= (A + B) + C;
Print("(A + B) + C = %0.16f\n", D);
K:= KahanSum(FA, 3);
Print("Kahan sum = %0.16f\n", K);
Print("(A + B) + C = %x %x\n", D);
Print("Kahan sum = %x %x\n", K);
]</syntaxhighlight>
{{out}}
<pre>
Epsilon = 1.1102230246251600E-016
(A + B) + C = 1.0000000000000000E+000
Kahan sum = 1.0000000000000000E+000
(A + B) + C = FFFFFFFF 3FEFFFFF
Kahan sum = 00000000 3FF00000
</pre>
 
=={{header|zkl}}==
zkl floats are C doubles.
{{trans|go}}
<langsyntaxhighlight lang="zkl">fcn kahanSum(numbers){
sum:=c:=0.0;
foreach x in (vm.arglist){
Line 1,741 ⟶ 3,041:
while(1.0 + e!=1.0){ e/=2 }
e
}</langsyntaxhighlight>
<langsyntaxhighlight 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");</langsyntaxhighlight>
{{out}}
<pre>0.99999999999999988898 Left associative. Delta from 1: 1.11022e-16
Jjuanhdez
2,122

edits

Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/252236...361311"