2015/11/26(木)kv-0.4.27
今回の変更点は、
- 公比が区間な等比級数の和の公式を修正し、それに伴いhypergeom.hppを修正。
- optimize.hppで、求まった区間のうち隣接しているものを結合する機能を追加。
- allsol.hppで、TRIMのアルゴリズムを少し変更。
- constants<affine<T>>::pi()などが機能していなかった問題を修正。
- webデモの問題記述でpi, e, ln2を使えるようにした。
- webデモにaffineとintervalによる関数の値域評価を追加。
2015/11/03(火)kv-0.4.26
今回は、ベキ級数演算で欠けていたasin, acos, atan, asinh, acosh, atanhを追加したのが大きな変更です。これで、これらの関数を含む数値積分や常微分方程式を解けるようになりました。超幾何関数やBessel関数も少し進化しましたが、これらはまだまだ試作段階です。また、Karush-Kuhn-Tucker条件を方程式に直して最適化問題の解を精度保証する試みのための機能をいくつか追加しました。
githubとかで公開できるようになるといいかなあと思いつつ、年寄りは新しいメディアに適応するのにエネルギーが要るのう。
2015/10/15(木)MATLAB 2015b での精度保証
- BLASを他のものに差し替える。速度は遅くなってしまうことが多い。
- matlabを2012aで留めておく。
- マルチコア使用を禁止し、シングルコアで動作させる。当然遅くなる。
で、最近出たmatlab 2015bですが、Intel MKLが再び丸め変更が効くようになったという噂があり、3年ぶりにちゃんと精度保証の出来るmatlabになったかと期待されていました。しかし、Rump先生のテストによると結果は残念なものだったそうです。確かにIntel MKLはちゃんと動くようになったそうですが、今度は行列a, bに対する a .* b のような演算で丸め変更が効かなくなってしまったそうです。このような簡単な計算で外部のBLASを呼ぶとは考えにくく、恐らくこれはmatlab内部で完結した処理と考えられます。matlab自身でマルチスレッドでデータを分割して計算するような実装になっていて、そのとき他スレッドに丸めの向きの情報を伝えていないと予想されます。
たまたまうまく行っていた「他人のふんどしで精度保証する」方法論はそろそろ限界なんじゃないでしょうか。自分のコントロールできる範囲内の事柄で精度が保証できている証拠を構成するようにアルゴリズム全体を設計するべきだと思います。
2015/10/15(木)kv-0.4.25
今回の主な変更は、
- 内部型の異なる区間の間の相互変換を行うinterval-conv.hpp追加
- eig.hppで重複固有値があるとうまく計算できなかった問題を修正
- ガウスの超幾何関数を追加
最近はほんの遊びで追加したmpfrラッパーが使われる頻度が高くなってきており、もう少しきちんと設計しなおした方がいいかもなあ、なんて考えています。それも含めて、混合精度演算をスマートに行える全体の構成をきちんと考える必要を感じています。
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の品質は上がっているのかもしれません。案外精度良かったんだな、という印象です。