最終更新: 2016/11/20

4倍精度演算ライブラリ

柏木 雅英

1. はじめに

4倍精度演算ライブラリは、例えば High-Precision Software Directory で公開されているDavid H. Baileyによるソフトウェア群と同様、 KnuthのtwosumアルゴリズムとDekkerのtwoproductアルゴリズムを用いて 倍精度(double)計算で擬似的に4倍精度計算を実現するものである。

本ライブラリ(dd.hpp) は、

の機能を持つ。これだけでは大したことはできず、rdd.hppの中にある を使ってinterval.hppと組み合わせて、4倍精度区間演算 を実現できるところに特徴がある。

2. ファイル構成

dd.hpp
(下請け: conv-dd.hpp, fpu53.hpp, convert.hpp, constants.hpp)
rdd.hpp
(下請け: hwround.hpp, rdd-hwround.hpp, rdd-nohwround.hpp)

3. 使い方

test-dd.hppを見てください。 区間の両端に使って4倍精度区間演算を行う例は、test-idd.hppを見てください。

4. 4倍精度演算

アルゴリズムは、 double-double演算、double-double区間演算まとめ (Version 2016/11/20) の4節までにまとめた。 (簡単には簡易メモ。)

5. 4倍精度区間演算

4倍精度区間演算は、丸めの向きの指定可能な4倍精度加減乗除、平方根を使って 実装されている。 それらのアルゴリズムは、 double-double演算、double-double区間演算まとめ (Version 2016/11/20) の5節にまとめた。

(2014/4/9に機能追加) -DKV_NOHWROUNDを付けてコンパイルすると、 ハード的な丸め方向の変更を行わず、方向付き丸めをエミュレートで実現 するようになる。

6. Intel CPUの過剰精度問題

IntelのCPUは、昔からあるFPUと、最近のCPUにのみ備わっているSIMD演算用のSSE2の 2種類の浮動小数点演算器を持っている。

FPUの内部精度は 全長80bit, 仮数部64bitで、IEEE754の全長64bit, 仮数部53bitと 異なることはよく知られている。 これと本ライブラリで使っている擬似多倍長法とは極めて相性が悪く、様々な 異常現象を引き起こす。例えばtwosumではいったん足して再度引くことによって エラーの計測を行っているが、なまじ中間変数の精度が高いために、正しくエラーが 測定できないというようなことが起こる。

SSE2の方はIEEE754に完全に従っている。

IntelのCPUでも、64bit環境だと基本的にSSE2のみを用いもはやFPUは使わないので、 このような異常現象は起きない。が、わざと-mfpmath=387とかすればその限りではない。 また、32bitCPUでも、-msse2 -mfpmath=sseとか付ければ安全になる。

小さな親切大きなお世話でFPUには困ったものだが、幸いなことにFPU Control Registerに FPUの計算精度を設定するフラグがあり、これを64bitから53bitに変更すれば この問題は避けられる。本ライブラリでは、fpu53.hをincludeして、 これを実行している。 コンパイラによってはこれがうまく動作しないなど考えられるので、何か問題があれば 報告して欲しい。

7. おわりに

dd.hpp単体で使ったときの数学関数 (近似計算) は面倒なのでまだ実装してないが、 やろうと思えばすぐなのでそのうち実装するつもり。精度保証には関係ないからなー。 interval<dd>に対しては数学関数は使えます! dd.hpp単体で使ったときの数学関数 (近似計算) を実装しました(2014年11月12日)。 interval<dd>に対しての精度保証された数学関数は元々使えます。

conv-dd.hppに含まれているstringtodd, ddtostringは、dd型と10進数文字列の 相互変換を行うもので、実装は面倒だった。 単体でも使えるので、ddの入出力に困ってる人はぜひ使ってみて欲しい。