2014/06/11(水)NAS2014
今回から、発表に関する資料はなるべく公開したいと考えていて、
ここに関連資料をアップロードします。
- nas2014.pdf 予稿の原稿です。mul_up/down(対策後)の誤植を修正しています。
- test-nohwround.cpp 丸め変更を行った場合と結果が一致することを確認するのに使ったプログラムです。
- bench.cpp 個々の演算のベンチマークに使ったファイルです。
- fivesol.zip 5つの解を持つ回路方程式の全解探索のベンチマークに使ったプログラムです。コンパイルオプションは、-O3 -DNDEBUG -DKV_NOHWROUND (最後のオプションを付けると丸めモードを変更しない本発表の方法)
- nas2014-slide.pdf 発表スライド
2014/06/06(金)kv-0.4.5
主な変更は、数値積分や常微分方程式のstep size controlに関する部分です。
また、digamma (gamma関数の対数微分) も追加しました。
digammaを作っていて内部で使っている数値積分の甘さが気になったのでちゃんと書きたくなった、
というのがstep size controlをいじった動機です。
digammaは元々gammaの改善が目的だったけれど、gammaはまだ特に変わっていません。
また、数値積分の使い方について詳細に記述しました。
他の機能についても詳細に記述したいところだけど、いつになることやら。
lobachevsky, gamma, digammaとそれぞれ異なる特異積分のテクニックを含んでおり、そのへんを汎用的な手法として整理するのもそろそろやりたいところだなあ。
2014/06/02(月)gamma関数とdigamma関数
gamma関数を以前作ったけど、微分が0になる位置を正確に計算できなかったので、
特に負の領域で満足の出来ないものだった。(大きい区間幅の入力に弱い)
某R先生もgamma関数を作ってるけど、こちらは微分が0になる位置を数表で持っていて、
それは数式処理ソフトで作っているようで、少しずるい感じ。
その数式処理ソフトは精度保証されているのか? 単なる高精度計算である可能性は?
精度保証を行うソフトウェアは、その計算結果が精度保証出来ていると思えるための理由を
そのソフトウェア自身の中に持っていて欲しい気がする。
そのへんを完全にするために、digamma関数を何となく作り始めて、何となく動いた気がした。
-digamma(1) = γ (0.5772...=オイラー定数)らしいので、オイラー定数の精度保証もついでに出来たことになる?
2014/06/01(日)kv-0.4.4
定積分のライブラリに初期化し忘れがあり、例えばODEを解かせたあとにそのまま
定積分を行うと結果がおかしくなることがありました。
精度保証ライブラリとしては、計算結果の信頼性に影響を及ぼす可能性のあるバグはみっともないので、
こういう類のバグがあった場合は直ちにアップデートするつもりでいます。
また、Airy関数を試作して遊んでみました。Airy関数の満たす常微分方程式を初期値問題として
解くだけなので、絶対値の小さい範囲しか計算できないと思われます。
Ai(-100) = [0.17675339323949076,0.17675339323961479]
2014/05/27(火)HaswellとFMA
ついでにvmwareを最新にし、vmwareの中で飼ってるubuntuを14.04に。
CPUがSandyBridgeからHaswellになったので、ちょっとFMAで遊んでみました。
Haswellの新機能の一つに、FMA命令の追加があります。
FMA = Fused Multiply Add です。乗算と加算のフュージョンですね。
a * b + cをいっぺんに計算してくれる機能で、内積計算なんかが高速になりそうですが、
精度保証的には「a*bをいったん無誤差で計算してからcを加算し、最後に一回だけ丸める」という
動作が重要です。これを使うとtwoproductが非常に簡単に書けたりします。
で、ちょっと動くかどうか試してみました。
#include <stdio.h>
#include <immintrin.h>
double fma(double a, double b, double c)
{
double d;
__m128d va = _mm_set_sd(a);
__m128d vb = _mm_set_sd(b);
__m128d vc = _mm_set_sd(c);
__m128d vd = _mm_fmadd_sd(va, vb, vc);
_mm_store_sd(&d, vd);
return d;
}
int main()
{
printf("%f\n", fma(1., 2., 3.));
}
こんな感じで動きました。gccだと-mfmaとか-march=nativeとか必要。VCは特に要らない?AVXなので4並列で実行する能力がありますが、このコードは1つだけ計算してます。
1つだけ計算するのに128bitレジスタを使わなきゃいけないのはちょっと違和感?
twoproductに用いたときの性能評価とかはこれからぼちぼちやっていきましょう。