Kahan と ruby の sum

これは何？

こんな記事

を書くに当たり。

Kahan のアルゴリズムを Wikipedia をもとに実装してみたけど、 ruby の Enumerable#sum の結果と合わなかった。

何が違うのか調べてみた。

Wikipedia の Kahan

記事は こちら 。

ruby で書くとこう。

ruby

def kahan(a)
    sum=c=0.0
    a.each do |v|
        y = v-c
        t = sum+y
        c = (t-sum) - y
        sum = t
    end
    sum
end

実行してみると。

ruby

ary = [1e18, 1.0, -1e18]
p( { kahan:kahan(ary), sum:ary.sum } )
# => {:kahan=>0.0, :sum=>1.0}

合わない。
上記の独自実装が Enumerable#sum に負けている。

本家のソースを見る

ソースはこちら。

論文へのリンクがあるけど、お金払わないと読めないようなのでソースを見る。

Wikipedia との違いは

	Wikipedia	ruby
名前	Kahan	Kahan-Babuska
誤差の扱い	随時総和に算入	誤差を総和とは別に積んでいってあとで総和に算入

という感じ。

ruby で書くとこう。

ruby

def kahan_babuska(a)
    sum = c = 0.0
    a.each do |v|
        t = sum+v
        c += v.abs <= sum.abs ? sum-t+v : v-t+sum
        sum=t
    end
    sum+c
end

実行してみると

ruby

ary = [1e18,1,-1e18]
p( { kahan:kahan(ary), kahan_babuska:kahan_babuska(ary), sum:ary.sum } )
# => {:kahan=>0.0, :kahan_babuska=>1.0, :sum=>1.0}

おお。合っている。

もっと精度を良くする

ruby の sum は大変良いんだけど、もちろん誤差が出ることもある。

ruby

ary = [1e40, -1e20, 1, 1e20, -1e40]
p( { kahan:kahan(ary), kahan_babuska:kahan_babuska(ary), sum:ary.sum } )
#=> {:kahan=>0.0, :kahan_babuska=>0.0, :sum=>0.0}

1e40, -1e20, 1, 1e20, -1e40 は、 kahan_babuska をもってしても 1 にはならない。

で。
先程の kahan_babuska の実装を見ると、 c に対して += を繰り返している。
これは kahan_babuska の出番だ。

というわけで、 kahan_babuska の誤差項に kahan_babuska を適用してみた。

ruby

def kb2(a)
    sum = c0 = c1 = 0.0
    a.each do |v|
        t = sum+v
        v0 = v.abs <= sum.abs ? sum-t+v : v-t+sum
        t0 = c0 + v0
        c1 += v0.abs <= c0.abs ? c0-t0+v0 : v0-t0+c0
        c0 = t0
        sum=t
    end
    sum+c0+c1
end

ary = [1e40, -1e20, 1, 1e20, -1e40]
p( { "+":ary.inject(&:+), sum:ary.sum, kb2:kb2(ary) } )
#=> {:+=>0.0, :sum=>0.0, :kb2=>1.0}

おお。二段 kahan-babuska ならちゃんと 1 になる。

計算量

ループ内の処理を比べてみた。

	加算	絶対値	比較
単純加算	1	0	0
Kahan	4	0	0
Kahan-Babuska	4	2	1
2段 Kahan-Babuska	7	4	2

そんなに誤差が嫌なら浮動小数点数を使うべきではないようにも思うけど、必要な場面に遭遇したら払う気になるコストではあると思う。