中点はどうやって計算するべきか

2016/12/28(水)中点はどうやって計算するべきか

今回は非常に単純な話で、IEEE754のdoubleで表現された数a,bが与えられたとき、その中点(a+b)/2はどう計算するのが最善でしょうか、という問題です。もちろん誤差は入りますが、出来ればnearestにしたい、すなわち、真の(a+b)/2をIEEE754のルールで丸めたものを計算することを目標にします。IEEE754の内部表現は2進数ですから、2で割る操作は指数部を1減らすだけなので普通は誤差は入りません。しかし、極端な場合を考えると案外難しいのです。

(a+b)*0.5

まずは普通に足して2で割る方法。「2で割る」は「0.5を掛ける」にした方がコンパイラに優しいでしょう。ごく普通の方法で安定度抜群ですが、例えば

a=2¹⁰²³
b=2¹⁰²³

の場合、a+bがオーバーフローしてしまう問題があります。

a0.5+b0.5

オーバーフローを避けて先に2で割る方法です。計算量は少し増えてしまいますが。こちらは、アンダーフロー領域で問題になります。例えば、

a=2^-1074
b=2^-1074

のとき、a*0.5、b*0.5はともにアンダーフローで0になってしまい、結果も0になってしまいます。

a+(b-a)*0.5

これもよく見る式です。IEEE754では関係ありませんが、内部演算が10進数の場合にこうやって計算することが推奨されているのを見かけたことがあります。しかし、

a=2¹⁰²³
b=-2¹⁰²³

でオーバーフローしてしまいます。

a+b0.5-a0.5

上の式でオーバーフローを避けたものです。今度こそ良さそうですが、

a=5*2^-1074
b=7*2^-1074

としてみましょう。このとき、IEEE754の偶数丸めルールの関係で、

a*0.5=2*2^-1074
b*0.5=4*2^-1074

となります。よって計算結果は7*2^-1074となってしまい、真の6*2^-1074になりません。

a+b0.5-a0.5 (ただし上向き丸めで計算)

INTLAB version 9で採用している計算方法です。a*0.5とb*0.5に入る丸めの方向が逆向きになることが無いので、先の例だとうまく計算できます。しかし残念ながら必ずnearestにはならないようで、Mさんがちゃんと計算できない例を発見しました。

a=(2⁵³-1)*2^-1073
b=7*2^-1074

とすると、

(nearest)=2^-1021+2^-1073
(計算値)=2^-1021+2^-1072

と値がずれてしまいます。a*0.5は無誤差ですが、b*0.5と加算で両方上向き丸めになった結果、真値(2^-1021+2^-1073+2^-1075)を丸めたものとずれてしまいました。

結論?

kvライブラリの区間の中点を計算するmid関数では、aとbの絶対値が両方共1より大きいときはa*0.5+b*0.5で、それ以外のときは(a+b)*0.5で計算しています。このようにifで場合分けをしていいなら話は簡単なのですが、matlabの場合はベクトル単位で一括処理しないとパフォーマンスが出ないのでなるべくifを使いたくないという事情があります。if文禁止で完璧な中点の計算アルゴリズムはまだ見つかっていません。誰か考えて！