最終更新: 2022/1/3

4倍精度演算ライブラリ

柏木雅英

1. はじめに

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

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

加減乗除 (近似計算)
平方根 (近似計算)
等号、不等号
abs, floor, frexp
exp, log, sin, cos, tan, sinh, cosh, tanh, asin, acos, atan, asinh, acosh, atanh, expm1, log1p, powなどの数学関数
10進文字列との相互変換

の機能を持つ。これだけでは大したことはできず、rdd.hppの中にある

丸めの向きの指定可能な加減乗除
丸めの向きの指定可能な平方根
方向付きの10進文字列との相互変換

を使って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 2018/6/9) の4節までにまとめた。 (簡単には簡易メモ。)

5. 4倍精度区間演算

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

6. コンパイルオプション

コンパイルオプションで高速化できる場合があり、詳細は 20. 丸めモードの変え方とコンパイルオプションまとめを見て欲しい。

7. Intel CPUでの32bit環境での問題

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して、これを実行している。コンパイラによってはこれがうまく動作しないなど考えられるので、何か問題があれば報告して欲しい。

(version 0.4.54で仕様を変更) ddの演算が正しく動くためにははIEEE754に完全に従っている必要があり、例えばIntel CPUの32bit環境はIEEE754に従わないことがあるので、本ライブラリはコンパイルできないようになっている。具体的には、#include <cfloat>した後マクロFLT_EVAL_METHODを調べ、それが0でないならばエラーでコンパイルできないようにした。 Intel CPUの32bit環境でどうしてもkvを使いたい場合は、コンパイルオプションとして -msse2 -mfpmath=sse を付けるなどして SSE2を使わせるようにすると、上記マクロが0になってコンパイルできるようになる。このマクロが定義されないコンパイラ、あるいはSSE2を持たないCPUは、古すぎるためサポートしないことにした。

8. frexpの問題点

frexpは、

    double x, y;
    int i;
    y = frexp(x, &i);

のように呼び出すと、xを入力として 0.5 ≤ |y| < 1、x=y×2ⁱを満たすような yとiを計算してくれる関数である。 (IEEE754の指数部とは1ずれるが) 指数部を取り出してくれる関数としてよく用いられる。 ddに対しても、

    kv::dd x, y;
    int i;
    y = frexp(x, &i);

のように同様な関数を実装してある。

しかし、例えばx=1+2^-1074のように ddの上位と下位の大きさが大きく離れている場合を考えると、 frexpの出力はy=0.5+2^-1075、i=1となるべきであるが、 2^-1075がdoubleで表現できないためそういう出力は不可能であるという問題点がある。すなわち、doubleではfrexpは誤差の発生しない関数であるが、 ddでは誤差が発生する場合があることに注意が必要である。なお、指数部iの方は常に信用できる。

9. おわりに

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

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

IEEE 754のdoubleと同様な無限大が実装されている。内部表現は、(±∞, 0)である。