2015/10/11(日)数学関数の精度を調べる
現代の一般的なプログラミング環境では、浮動小数点演算はIEEE 754 Std.に基づいて行われるのが一般的です。そしてIEEE 754 Std.においては、演算に完全精度 (フォーマットで表現可能な範囲内で最善の精度、0.5ulp) を要求しているのは
- 加算
- 減算
- 乗算
- 除算
- 平方根
従って、C言語などから使うexp, sinなどの数学関数の返す値は信用することは出来ません。kvライブラリの区間演算でも数学関数は全て自力で実装しています。区間演算ライブラリの一つである PROFIL/BIAS では数学関数をC言語の組み込み関数で計算した後succ/predで区間を広げて計算結果としていますが、これは誤りであり、精度保証付き数値計算を行うためのツールとしては役に立ちません。
Intelの64bit環境では、数学関数は全てソフトウェアで計算されるのが一般的です。従って、計算精度はlibmの実装に依存することになります。たとえ完全精度を出す義務が無かったとしても、長い間のアルゴリズムの改良により完全精度が出るようになっているかも知れません。そこで、
- CPU: Intel Core M 5Y70
- OS: Ubuntu 14.04 LTS 64bit
| 関数名 | 下限 | 上限 | 結果 |
|---|---|---|---|
| exp | -700 | 700 | OK |
| log | 1e-308 | 1e308 | OK |
| sin | -1e10 | 1e10 | OK |
| cos | -1e10 | 1e10 | OK |
| tan | -1e10 | 1e10 | OK |
| asin | -1 | 1 | OK |
| acos | -1 | 1 | OK |
| atan | -10 | 10 | OK |
| sinh | -700 | 700 | NG |
| cosh | -700 | 700 | OK |
| tanh | -700 | 700 | NG |
| asinh | -1000 | 1000 | OK |
| acosh | 1 | 1000 | OK |
| atanh | -1 | 1 | NG |
sinh, tanh, atanhについて、精度不足となった例を示しておきます。
sinh(0.74251904152333736419677734375)=0.8126541801799147535234624228905886...
と計算されたが、真の値はdoubleの区間
[0.8126541801799148645457648854062426,0.8126541801799149755680673479218968]
の間にある。
tanh(0.48066521994769573211669921875)=0.4467762045341044374602290645270841...
と計算されたが、真の値はdoubleの区間
[0.4467762045341044929713802957849111,0.4467762045341045484825315270427382]
の間にある。
atanh(0.2377787786535918712615966796875)=0.2424184572356452571639806592429522...
と計算されたが、真の値はdoubleの区間
[0.2424184572356452016528294279851252,0.2424184572356452294084050436140388]
の間にある。
もう記録は残っていませんが、10年以上前に調べたときはもっと精度に問題があった記憶があるので、libmの品質は上がっているのかもしれません。案外精度良かったんだな、という印象です。
2015/10/05(月)kv-0.4.24
今回の主な変更点は、
- ライセンスを明記し、MITライセンスとした。
- interval.hppのproper_subsetを修正し、proper_subset([0,1],[0,2])が真になるようにした。
- beta.hppを追加。
- eig.hppを追加(高安 亮紀氏の提供による)。
ライセンスは今まで曖昧にしていましたが、某所から問い合わせがあって、MITライセンスと明記することにしました。大雑把に言えば、著作権表示を残しさえすれば自由にしていいよ、責任は取らないよ、というライセンスと理解しています。「煮るなり焼くなり好きにして」という気分だったのですが、やはり明記されていないと使いにくいということはありそうです。
proper_subset(X,Y)は、XがYの真部分集合のとき真、という関数ですが、X=[0,1], Y=[0,2]のようなとき偽を返していたという問題がありました(X=[0.5,1]なら真、すなわち両側に余裕がないと真になっていなかった)。この関数を使っていたのはKrawczyk法の部分だけですが、今まではこのバグによりわずかに性能が低下していた可能性があります。
eig.hppは、非対称実行列の固有値と固有ベクトルを(近似的に)計算するeig関数と、非対称実行列の固有値を精度保証付きで計算するveig関数を含んでいます。高安 亮紀氏に提供していただきました。ありがとうございます。
(2015年10月6日 追記)
kv-0.4.15、シンプルな区間演算ライブラリで公開した簡単な区間演算ライブラリを、0.4.24相当にアップデートしました。interval-simple-0.4.24.tar.gz
特徴を再掲しておきます。
- kvライブラリの区間演算の部分を切り出して、簡単な区間演算ライブラリを作ってみました。
- テンプレートを排除しました。従って内部型はdoubleのみですが、プログラムが読みやすくなっています。
- テンプレートを排除したおかげでboostに依存しません。単体で使えます。
- 必要なファイルはinterval.hpp一つです。
- 内部型のカスタマイズが出来ないこと以外は、kvライブラリに含まれるinterval.hppと機能は同等です。数学関数も多数あります。
- 丸めの変更は単純にfesetroundを使っています。
- 文字列との方向丸め付き相互変換機能もあります。従って、文字列からその文字列の表す数値を含むように区間を初期化できるし、表示も真の値を含むように外側に丸めるようになっています。
- volatileを使った最適化対策はきちんとやってます。
2015/09/10(木)応用数理学会年会
2015/08/06(木)kv-0.4.23
今回は、allsol.hppにある非線形方程式の全解を探索するプログラムに、[0,∞]や[-∞,∞]などの無限大を含む領域を探索区間として指定できるようにしました。もっとも、必ず全ての解を求めて停止するわけでは無く、「求まる問題なら求まる」といったレベルです。ある問題を解いていてたまたま必要になったので実装したのですが、まだまだ改良が必要と思います。こうして文章にすると簡単そうですが、実は結構苦戦しました。無限大は難しい。
その他、vleq.hppの仕様を少し変更したり、区間のmidを精密にしたりしました。ベキ級数演算のpdfのドキュメントにも少し加筆しました。
2015/06/29(月)kv-0.4.22
今回は、defint_power_r, defint_log_r, defint_power_autostep_r, defint_log_autostep_rの4つの端点特異性を持つ積分の関数のbug fixです。これら右側の端点特異性を扱う関数は、時間を反転させて左側versionに任せているのですが、その時間反転部分にバグがありました。これらの関数を使う例題を作っていなかったので、気付きませんでした。
その他、細かいところをいくつか修正しています。