2014/01/25(土)Error Free Transformation (EFT) is NOT error-free
def twosum(a, b) :
x = a + b
tmp = x - a
y = (a - (x - tmp)) + (b - tmp)
return x, y
def split(a) :
tmp = a * (2**27 + 1)
x = tmp - (tmp - a)
y = a - x
return x, y
def twoproduct(a, b) :
x = a * b
a1, a2 = split(a)
b1, b2 = split(b)
y = a2 * b2 - (((x - a1 * b1) - a2 * b1) - a1 * b2)
return x, y
python風に書くとこんな感じ。x,y = twosum(a, b) は、xがa+bを浮動小数点数で普通に計算した値、yはその誤差分になり、x+y=a+bが厳密に成立するとされています。
x,y = twoproduct(a, b) は、xがa*bを浮動小数点数で普通に計算した値、yはその誤差分になり、x+y=a*bが厳密に成立するとされています。
しかし、オーバーフローやアンダーフローに注意が必要です。
xがオーバーフローしたら誤差も何も無いのでまあよしとしましょう。
アンダーフローは、twosumの場合、IEEE754標準のように非正規化数を用いた漸近アンダーフローを備えた浮動小数点数では原理的に発生しないので問題ありません。
しかし、twoproductの場合は、yが非正規化数の領域に入る場合にはyの有効桁数が失われてしまい、無誤差ではなくなってしまいます。すなわち、|x| < 2^(-969) (=-1074+106-1) の場合は誤差が発生します。
また、twoproductで使われているsplitにも注意が必要で、(2^27+1)を乗じている部分でオーバーフローする可能性があります。ただしこれは実装上の工夫でどうにでも回避可能。
ここまではよく知られている問題点ですが、今までにあまり注意を払われて来なかったと思われる問題点を見つけたのでメモしておきます。
まず、twoproductで
a = 6.929001713869936e+236
b = 2.5944475251952003e+71
のケース。x=a*bはオーバーフローしませんが、a1*b1がオーバーフローしてしまい正しくyが計算されません。2つ目は、twosumで、
a = 3.5630624444874539e+307
b = -1.7976931348623157e+308
のケース。x=a+bはオーバーフローしませんが、中間変数tmpがオーバーフローしてしまい正しくyが計算されません。twosum, twoproductを用いた方向付き丸めのエミュレート法を考えていて、これらを見つけました。
精度保証付き数値計算では一つの例外も許されないので、こういう現象には本当に神経を使います。他にまだ気付いていない問題点が存在しないと言い切れるのでしょうか。
2014/01/02(木)あけおめ ＆ kv-0.4.1
正月で時間があるので、現在のライブラリをkv-0.4.1として公開します。
- 区間演算のpowで奇数乗の場合に区間が大きかったのを修正(Thanks to 山口俊樹君)
- webデモに二重積分と常微分方程式の初期値問題を追加
- 丸めモード変更のエミュレートとガンマ関数の試作
まだ研究室の学生くらいしか使っていないと思われるけど、使われて鍛えられるほど完成度が高くなるのと思われるので、どんどん使って苦情を言ってもらえるとありがたいです。
現状で根本的に欠けているのが線形計算まわり。ただ、ここいらへんは高速化しようとするとLAPACKやBLASと連携する必要があり、いろいろ面倒くさい。double以外もdoubleと同等に使えるようにしようとすると更に面倒に。研究的に関心の薄いジャンルなのでどうしても後回しになってしまうけど、これをやらないとなかなか使ってもらえないだろうなあ。
2013/12/26(木)ちょっとだけ
http://verifiedby.me/
に、
精度保証付き数値計算の必要性
を追加した。たまたま通常の数値計算ソフト(matlab)と精度保証付き数値計算の比較の図を作成する必要に迫られたのでちょっとした文章をつけてまとめてみた。
ほんとは「精度保証付き数値計算とは何ぞや」というもうちょっとちゃんとした文章が欲しいところだがなあ。
2013/12/03(火)kv-0.4 & ライブラリのデモ
もう少し書き溜めてから上げようと思っていたけど、現状の例題の範囲では問題はないものの精度保証的に問題がありそうなバグが見つかったので、直したものを公開することに。
ついでに、まだ試作品だが、web上で精度保証を実行できるデモを作成したので、公開してみた。
KVライブラリのデモ
中間コードを生成するわけではなく、内部でC++コンパイラを呼び出すという無茶仕様なので、計算よりもコンパイルに時間がかかるのが問題。サーバをもう少し速いマシンに変えたい。大勢が使って負荷がどうなるのかとか、セキュリティ的にやばそうとかいろいろ気になるが、まあ暫定ということで。とりあえず全解探索と数値積分。ODEとか二重積分とかも暇があったら追加したい。
2013/11/04(月)丸めモードの変更
丸めモードの変更に関しては気を付けなければいけないポイントが大変多く、講義なんかするときにもすぐに忘れそうになるので、IntelのCPUに関してまとめた資料を作ったので、ここにアップしておく。
roundingmode.pdf
なお、kvライブラリにおいては、この資料の7章にあるようにコンパイラに備わっているfesetroundまたは_controlfpを用いており、最速というわけではない。が、Intel以外のCPUを含む割と多くの環境で安定して動くというメリットもある。
IntelのCPUで64bitモードの場合に限って、KV_FASTROUNDというマクロを定義すると、Intrinsicを用い、プログラム起動時に予めレジスタ値を計算し丸め変更時には書き込むだけという速くなりそうな方法を使うようになる。FPUは使わないと仮定してSSE2の丸めモードしか変更しないし、MXCSRレジスタの丸めモード以外の部分の値を全くケアしないので、そのことによる副作用をよく理解した上で使う必要がある。
なお、プログラミング言語や環境によっては丸めモードの変更をしづらいこともあり、丸めモードの変更を行わずに精度保証する方法の研究も進んでいることを付記しておく。