Kahan summation: Difference between revisions
Grouping BASIC dialects
(Grouping BASIC dialects) |
|||
(34 intermediate revisions by 22 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 '''Kahan summation algorithm''' is also known as the '''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:
<
while 1.0 + epsilon != 1.0:
epsilon = epsilon / 2.0</
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.
<
# only 6 digits of precision are required and assuming INT is 32 bits, we can #
# emulate this fairly easily #
Line 191 ⟶ 281:
SKIP
)
</syntaxhighlight>
{{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++}}
<
#include <stdlib.h>
Line 235 ⟶ 705:
return 0;
}</
{{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++}}==
<
float epsilon() {
Line 274 ⟶ 783:
return 0;
}</
{{out}}
<pre>Epsilon = 5.96046e-08
Line 280 ⟶ 789:
Kahan sum = 1</pre>
=={{header|
{{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
end
eps
end
def kahan(nums)
sum = 0.0_f32
c = 0.0_f32
nums.each do |num|
y = num -
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
Kahan sum
</pre>
=={{header|D}}==
<
float kahanSum(float[] fa...) {
Line 341 ⟶ 851:
writefln("(a + b) + c = %0.8f", ((a + b) + c));
writefln("Kahan sum = %0.8f", kahanSum(a, b, c));
}</
{{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.
<
;; the maximum number of decimals is 17. Floating point arithmetic is not always 100% accurate
;; even with simple data :
Line 388 ⟶ 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.</
=={{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.
<
let rec δ λ =
match λ+1.0 with
Line 418 ⟶ 966:
printfn "Sum: %e" (input |> List.sum)
printfn "Kahan: %e" (input |> Σ)
0</
{{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 429 ⟶ 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.
<
COMPENSATED SUMMATION. C WILL NOT STAY ZERO, DESPITE MATHEMATICS.
DIMENSION A(12345)
Line 450 ⟶ 1,016:
1 FORMAT (6HSUM = ,F12.1)
END
</syntaxhighlight>
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 457 ⟶ 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...
<
REAL X !The number.
REAL A !Its base ten log.
Line 522 ⟶ 1,088:
1 FORMAT (A1,"Sum = ",F12.1)
END
</syntaxhighlight>
Output is
<pre>SSum = 10005.8
Line 549 ⟶ 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:
<
CONTAINS
CHARACTER*76 FUNCTION FP8DIGITS(X,BASE,W) !Full expansion of the value of X in BASE.
Line 696 ⟶ 1,262:
1 FORMAT (A1,"Sum = ",F12.1)
666 FORMAT (A8,":",A)
END</
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 759 ⟶ 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:
<
REAL A(3)
REAL*4 X4E,X4L,X4S,X4C
Line 841 ⟶ 1,407:
665 FORMAT (A1,"Sum = ",F12.1)
666 FORMAT (A8,":",A)
END</
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 932 ⟶ 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:
<
import "fmt"
Line 964 ⟶ 1,530:
fmt.Println("Kahan summation: ", kahan(a, b, c))
fmt.Println("Epsilon: ", b)
}</
{{out}}
<pre>Left associative: 0.9999999999999999
Line 980 ⟶ 1,546:
Required code (for the current draft of this task):
<
epsilon=. 1.0
while. 1.0 ~: 1.0 + epsilon do.
Line 1,001 ⟶ 1,567:
sum=. t
end.
)</
Required results:
<
1
KahanSum a,b,c
1</
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 1,014 ⟶ 1,580:
First, the value of epsilon:
<
5.68434e_14</
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:
<
epsilon=. 1.0
while. 1.0 (~:!.0) 1.0 + epsilon do.
epsilon=. epsilon % 2.0
end.
)</
Now we have a smaller value for this explicit value of epsilon:
<
1.11022e_16</
With this new value for epsilon, let's try the problem again:
<
1
KahanSum a,b,c
1</
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:
<
0.9999999999999999
":!.16 KahanSum a,b,c
1</
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.
<
private static float kahanSum(float... fa) {
float sum = 0.0f;
Line 1,072 ⟶ 1,639:
System.out.println("Kahan sum = " + kahanSum(a, b, c));
}
}</
{{out}}
<pre>Epsilon = 5.9604645E-8
Line 1,080 ⟶ 1,647:
=={{header|Julia}}==
Julia can use its BigFloat data type to avoid floating point epsilon errors, as shown below.
<
function kahansum(arr)
Line 1,102 ⟶ 1,669:
println("Kahan sum is ", kahansum(v))
println("BigFloat sum is ", (BigFloat(a) + ep) + b)
</
<pre>
Epsilon is 1.1102230246251565e-16
Line 1,112 ⟶ 1,679:
=={{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:
<
fun kahanSum(vararg fa: Float): Float {
Line 1,139 ⟶ 1,706:
println("(a + b) + c = ${(a + b) + c}")
println("Kahan sum = ${kahanSum(a, b, c)}")
}</
{{out}}
<pre>Epsilon = 5.9604645E-8
Line 1,147 ⟶ 1,714:
=={{header|Lua}}==
{{trans|C}}
<
local eps = 1.0
while 1.0 + eps ~= 1.0 do
Line 1,175 ⟶ 1,742:
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)))</
{{out}}
<pre>Epsilon = 0.000000000000000111022302
Line 1,182 ⟶ 1,749:
=={{header|Modula-2}}==
<
FROM RealStr IMPORT RealToStr;
FROM Terminal IMPORT WriteString,WriteLn,ReadChar;
Line 1,244 ⟶ 1,811:
ReadChar;
END KahanSummation.</
{{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,254 ⟶ 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.)
<
epsilon()=my(e=1.); while(e+1. != 1., e>>=1); e;
\p 19
e=epsilon();
Kahan([1.,e,-e])</
{{out}}
<pre>%1 = 1.000000000000000000</pre>
=={{header|Perl
{{trans|Raku}}
<syntaxhighlight lang="perl">use strict;
use warnings;
use feature 'say';
sub kahan {
my(@nums) = @_;
my $summ = my $c = 0e0;
for
my $y = $num - $c;
my $t = $summ + $y;
Line 1,280 ⟶ 1,934:
}
my $eps = 1;
do { $eps /= 2 } until 1e0 == 1e0 + $eps;
say '
say 'Simple sum: ' . sprintf "%.16f", ((1e0 + $eps) - $eps);
say 'Kahan sum: ' . sprintf "%.16f", kahan(1e0,
{{out}}
<pre>Epsilon: 1.
Simple sum: 0.9999999999999999
Kahan sum: 1.0000000000000000</pre>
=={{header|Phix}}==
Line 1,296 ⟶ 1,950:
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>
Line 1,336 ⟶ 1,992:
</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
Line 1,372 ⟶ 2,030:
{"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>
Line 1,381 ⟶ 2,045:
</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>
Line 1,412 ⟶ 2,090:
=={{header|PHP}}==
'''Script'''
<
// Main task
Line 1,481 ⟶ 2,159:
$diff = (($a + $b) + $c) - kahansum($input);
echo $diff.PHP_EOL;
</syntaxhighlight>
'''Script output'''
<
Main task - Left to right summation: 10005.9
Main task - Kahan summation: 10005.9
Line 1,493 ⟶ 2,171:
bool(false)
Difference between the operations: 1.1102230246252E-16
</syntaxhighlight>
=={{header|Python}}==
===Python: Decimal===
<
>>>
>>> getcontext().prec = 6
Line 1,533 ⟶ 2,211:
>>> (a + b) + c
Decimal('10005.85987')
>>> </
===Python: Floats===
This was written as proof that the floating-point sub-task could work and is better off displayed, so...
<
>>> while 1.0 + eps != 1.0:
eps = eps / 2.0
Line 1,554 ⟶ 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)
>>> </
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,561 ⟶ 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.
<
>>>
>>> with localcontext() as ctx:
Line 1,576 ⟶ 2,254:
Simple sum is: 0.9999999999999999999999999999
Kahan sum is: 1.000000000000000000000000000
>>> </
=={{header|R}}==
Line 1,583 ⟶ 2,261:
Therefore, trying to use the main rules will not work as we can see in the code below
<
kahansum <- function(x) {
ks <- 0
Line 1,613 ⟶ 2,291:
kahansum(input)
# [1] 10005.86
</syntaxhighlight>
Apparently there is no difference between the two approaches. So let's try the alternative steps given in the task.
We first calculate '''epsilon'''
<
while ((1.0 + epsilon) != 1.0) {
epsilon = epsilon / 2.0
Line 1,623 ⟶ 2,301:
epsilon
# [1] 1.11022e-16
</syntaxhighlight>
and use it to test the left-to-right summation and the Kahan function
<
b <- epsilon
c <- -epsilon
Line 1,636 ⟶ 2,314:
kahansum(c(a, b, c))
# [1] 1
</syntaxhighlight>
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
<
# FALSE
Line 1,644 ⟶ 2,322:
((a + b) + c) - kahansum(c(a, b, c))
# [1] -1.110223e-16
</syntaxhighlight>
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
<
mepsilon <- function(divisor) {
guess <- 1.0
Line 1,662 ⟶ 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>
It would seem that it makes a differences what divisor we are using, let's look at the distribution using a text plotting library
<
txtboxplot(epsilons)
Line 1,692 ⟶ 2,370:
# +--+---------+----------+----------+---------+----------+----------+
# 0 2e-17 4e-17 6e-17 8e-17 1e-16
</syntaxhighlight>
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.
<
epsilon
# [1] 5.939024e-18
Line 1,714 ⟶ 2,392:
(a + b) + c == kahansum(c(a, b, c))
# TRUE
</syntaxhighlight>
=={{header|Racket}}==
Line 1,722 ⟶ 2,400:
Finally we try the alternative task version for single precision float point numbers.
<
(define (sum/kahan . args)
Line 1,738 ⟶ 2,416:
(displayln "Double presition flonum")
(+ 1000000.0 3.14159 2.71828)
(sum/kahan 1000000.0 3.14159 2.71828)</
{{out}}
<pre>Single presition flonum
Line 1,748 ⟶ 2,426:
Alternative task version
<
(let loop ([epsilon 1.f0])
(if (= (+ 1.f0 epsilon) 1.f0)
Line 1,758 ⟶ 2,436:
(+ 1.f0 epsilon.f0 (- epsilon.f0))
(sum/kahan 1.f0 epsilon.f0 (- epsilon.f0))</
{{out}}
<pre>Alternative task, single precision flonum
Line 1,764 ⟶ 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,769 ⟶ 2,476:
<br>of '''30''' (decimal digits) was chosen. The default precision for REXX is '''9''' decimal digits.
===vanilla version===
<
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
end /*while*/
say /*display a blank line before the fence*/
Line 1,782 ⟶ 2,489:
exit /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
kahan: procedure; $= 0; c= 0; do j=1 for arg()
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(
say ' b = ' left(
say ' c = ' left(
say 'simple summation of a,b,c = ' a+ b+
say 'Kahan summation of a,b,c = ' kahan(a, b, c)
return</
{{out|output|text= using PC/REXX and also Personal REXX:}}
<pre>decimal digits = 6
Line 1,867 ⟶ 2,574:
===tweaked version===
The following tweaked REXX version causes Regina (version 3.4 and later) to work properly.
<
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
end /*while*/
say /*display a blank line before the fence*/
Line 1,882 ⟶ 2,589:
exit /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
kahan: procedure; $= 0; c= 0; do j=1 for arg()
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(
say ' b = ' left(
say ' c = ' left(
say 'simple summation of a,b,c = ' a+ b+
say 'Kahan summation of a,b,c = ' kahan(a, b, c)
return</
{{out|output|text= using Regina version 3.4 and also for later versions:}}
<pre>decimal digits = 6
Line 1,913 ⟶ 2,620:
c = -3.15544362088404722164691426133E-30
simple summation of a,b,c = 1.0000000000000000000000000000001
Kahan summation of a,b,c = 1.0000000000000000000000000000000
</pre>
=={{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}}==
Arrays provide a "sum" method, which uses a Kahan-Babuska balancing compensated summation algorithm, according to the C-source.
<syntaxhighlight lang="ruby">epsilon = 1.0
epsilon /= 2 until 1.0 + epsilon == 1.0
a = 1.0
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>
=={{header|Scala}}==
===IEEE 754 Single precision 32-bit (JavaScript defaults to Double precision.)===
<
object KahanSummation extends App {
Line 1,964 ⟶ 2,786:
println("Kahan sum = " + kahanSum(a, ε, -ε))
}
}</
{{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:
<
namespace path ::tcl::mathop
Line 1,999 ⟶ 2,822:
puts "\tEpsilon is: [set e [epsilon]]"
puts "\tAssociative sum: [expr {1.0 + $e - $e}]"
puts "\tKahan sum: [kahansum 1.0 $e -$e]"</
<pre>Epsilon is: 1.1102230246251565e-16
Associative sum: 0.9999999999999999
Line 2,009 ⟶ 2,832:
The last stanza exercises the decimal package's different rounding modes, to see what happens there:
<
namespace path ::math::decimal
Line 2,043 ⟶ 2,866:
puts "\tAssociative sum $a + $b + $c: [asum $a $b $c]"
puts "\tKahan sum $a + $b + $c: [kahansum $a $b $c]"
}</
The results are a little surprising:
{{out}}
Line 2,070 ⟶ 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}}
<
sum:=c:=0.0;
foreach x in (vm.arglist){
Line 2,087 ⟶ 3,041:
while(1.0 + e!=1.0){ e/=2 }
e
}</
<
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");</
{{out}}
<pre>0.99999999999999988898 Left associative. Delta from 1: 1.11022e-16
|