Rosetta Code

Kahan summation: Difference between revisions

Page Discussion

Kahan summation (view source)

Revision as of 16:23, 4 June 2021

20,342 bytes added ,  2 years ago
m
→‎{{header|Phix}}: added syntax colouring the hard way
m (→‎{{header|Phix}}: added syntax colouring the hard way)
Line 1,554:
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 aphixonline)-->
<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>
-- NB: only suitable for this very specific example
<span style="color: #000080;font-style:italic;">-- NB: only suitable for this very specific example</span>
string s = trim(sprintf("%6g",a))
<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>
{{a}} = scanf(s,"%f")
<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>
return a
<span style="color: #008080;">return</span> <span style="color: #000000;">a</span>
end function
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
 
function roundsum(sequence s)
<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>
atom res = 0
<span style="color: #004080;">atom</span> <span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span>
for i=1 to length(s) do
<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>
res = round6(res+s[i])
<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>
end for
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
return res
<span style="color: #008080;">return</span> <span style="color: #000000;">res</span>
end function
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
function kahanroundsum(sequence s)
<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>
atom tot = 0, c = 0
<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>
for i=1 to length(s) do
<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>
atom x = s[i],
<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>
y = round6(x-c),
<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>
t = round6(tot+y)
<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>
c = round6(round6(t-tot)-y)
<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>
tot = t
<span style="color: #000000;">tot</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">t</span>
end for
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
return tot
<span style="color: #008080;">return</span> <span style="color: #000000;">tot</span>
end function
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
 
procedure main()
<span style="color: #008080;">procedure</span> <span style="color: #000000;">main</span><span style="color: #0000FF;">()</span>
sequence s = {10000.0, 3.14159, 2.71828}
<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>
string fmt = "%5.1f\n"
<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>
printf(1,fmt,roundsum(s))
<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>
printf(1,fmt,kahanroundsum(s))
<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>
end procedure
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span>
main()</lang>
<span style="color: #000000;">main</span><span style="color: #0000FF;">()</span>
<!--</lang>-->
{{out}}
<pre>
Line 1,594 ⟶ 1,596:
</pre>
Alternative task using floats and the explicitly calculated epsilon:
<!--<lang Phix>function kahansum(sequence sphixonline)-->
<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>
atom tot = 0, c = 0
<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>
for i=1 to length(s) do
<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>
atom x = s[i],
<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>
y = x - c,
<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>
t = tot + y
<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>
c = (t - tot) - y
<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>
tot = t
<span style="color: #000000;">tot</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">t</span>
end for
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
return tot
<span style="color: #008080;">return</span> <span style="color: #000000;">tot</span>
end function
<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>
function epsilon()
<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>
atom e=1
<span style="color: #000000;">b</span> <span style="color: #0000FF;">:=</span> <span style="color: #000000;">epsilon</span><span style="color: #0000FF;">(),</span>
while 1+e!=1 do
<span style="color: #000000;">c</span> <span style="color: #0000FF;">:=</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">b</span>
e/=2
<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>
end while
<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>
return e
<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>
end function
<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>
procedure main2()
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span>
atom a = 1.0,
<span style="color: #000000;">main2</span><span style="color: #0000FF;">()</span>
b := epsilon(),
<!--</lang>-->
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
Line 1,630 ⟶ 1,634:
{"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,639 ⟶ 1,649:
</pre>
Above and beyond:
<!--<lang Phix>function NeumaierSum(sequence sphixonline)-->
<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>
atom res = s[1],
<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>
c = 0.0 -- A running compensation for lost low-order bits.
<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>
for i = 2 to length(s) do
<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>
atom t = res + s[i]
<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>
if abs(res) >= abs(s[i]) then
<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>
c += (res - t) + s[i] -- If res is bigger, low-order digits of s[i] are lost.
<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>
else
<span style="color: #008080;">else</span>
c += (s[i] - t) + res -- Else low-order digits of res are lost
<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>
end if
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
res = t
<span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">t</span>
end for
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
return res + c -- Correction only applied once at the very end
<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>
end function
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
 
procedure main3()
<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>
sequence s = {1e300,1,-1e300}
<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>
printf(1,"Simple sum: %d\n",sum(s))
<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>
printf(1,"Kahan sum : %d\n",kahansum(s))
<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>
printf(1,"Neumaier sum : %d\n",NeumaierSum(s))
<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>
end procedure
<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>
main3()</lang>
<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>
<!--</lang>-->
{{out}}
<pre>
Petelomax
7,795

edits

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