Kahan summation: Difference between revisions
added Phix
Thundergnat (talk | contribs) m (→{{header|Perl 6}}: Hmm. not really sure what inaccuracies are in the example...) |
(added Phix) |
||
Line 1,065:
Simple sum: 0.9999999999999999
Kahan sum: 1.0000000000000000
</pre>
=={{header|Phix}}==
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.
<lang Phix>function round6(atom a)
-- NB: only suitable for this very specific example
string s = trim(sprintf("%6g",a))
{{a}} = scanf(s,"%f")
return a
end function
function roundsum(sequence s)
atom res = 0
for i=1 to length(s) do
res = round6(res+s[i])
end for
return res
end function
function kahanroundsum(sequence s)
atom tot = 0, c = 0
for i=1 to length(s) do
atom x = s[i],
y = round6(x-c),
t = round6(tot+y)
c = round6(round6(t-tot)-y)
tot = t
end for
return tot
end function
procedure main()
sequence s = {10000.0, 3.14159, 2.71828}
string fmt = "%5.1f\n"
printf(1,fmt,roundsum(s))
printf(1,fmt,kahanroundsum(s))
end procedure
main()</lang>
{{out}}
<pre>
10005.8
10005.9
</pre>
Alternative task using floats and the explicitly calculated epsilon:
<lang Phix>function kahansum(sequence s)
atom tot = 0, c = 0
for i=1 to length(s) do
atom x = s[i],
y = x - c,
t = tot + y
c = (t - tot) - y
tot = t
end for
return tot
end function
function epsilon()
atom e=1
while 1+e!=1 do
e/=2
end while
return e
end function
procedure main2()
atom a = 1.0,
b := epsilon(),
c := -b
sequence s = {a,b,c}
string fmt = iff(machine_bits()=32?"%.16f\n":"%.19f\n")
?{"Epsilon:",b}
printf(1,"Simple sum: "&fmt,sum(s))
printf(1,"Kahan sum : "&fmt,kahansum(s))
end procedure
main2()</lang>
{{out}}
32-bit
<pre>
{"Epsilon:",1.110223024e-16}
Simple sum: 0.9999999999999998
Kahan sum : 1.0000000000000000
</pre>
64-bit
<pre>
{"Epsilon:",5.421010862e-20}
Simple sum: 0.9999999999999999999
Kahan sum : 1.0000000000000000000
</pre>
Above and beyond:
<lang Phix>function NeumaierSum(sequence s)
atom res = s[1],
c = 0.0 -- A running compensation for lost low-order bits.
for i = 2 to length(s) do
atom t = res + s[i]
if abs(res) >= abs(s[i]) then
c += (res - t) + s[i] -- If res is bigger, low-order digits of s[i] are lost.
else
c += (s[i] - t) + res -- Else low-order digits of res are lost
end if
res = t
end for
return res + c -- Correction only applied once at the very end
end function
procedure main3()
sequence s = {1e300,1,-1e300}
printf(1,"Simple sum: %d\n",sum(s))
printf(1,"Kahan sum : %d\n",kahansum(s))
printf(1,"Neumaier sum : %d\n",NeumaierSum(s))
end procedure
main3()</lang>
{{out}}
<pre>
Simple sum: 0
Kahan sum : 0
Neumaier sum : 1
</pre>
| |||