kashiの日記

• 特殊関数の説明文を書いた。特にガンマ関数の計算方法の詳細を書いた。
• 内部で例外として使っていたrange_errorをdomain_errorに変更。
• 数学関数を含む常微分方程式が正しく解けないことがあったバグを修正。

• INTLAB - INTerval LABoratory
• INTLAB - INTerval LABoratory

1. 丸めの向きを下向きに変える
2. C₁=ABをBLASで計算する。
3. 丸めの向きを上向きに変える。
4. C₂=ABをBLASで計算する。
5. (丸めの向きを元に戻す)

BLASの中でマルチスレッドあるいはマルチコア計算を行っており、他のスレッドまたはコアに丸めの向きの情報が伝わらない

BLASの中でStrassenアルゴリズムなどの「減算を含む」計算が使われている

• Smithの定理を使ったDKAの精度保証版を作ってみた。
• gamma関数を少し整備した。
• defint-singular.hppの中にひっそりとある端点特異性を持つ関数の積分に、ステップ幅を自動調節するversionを追加した。

• この積分は全体を超高精度で行えばちゃんとまともな精度保証結果が出ることは分かっているが、それは現実的に無理。
• 高精度が必要なのは全体のある一部分だが、kvライブラリはworking precisionを全体で一斉に変えることを想定していて、一部分だけ計算精度を変える混合精度計算が書きづらい。

#include <stdio.h>
#include <math.h>

int main()
{
        double u, x, y, z;
        int i;

        u = 1;
        for (i=0; i<52; i++) u /= 2;

        x = 1+u;
        y = 1-u;
        z = -1;

        printf("%.15g\n", fma(x, y, z));
}

• FMA命令がちゃんと機能していれば、真の値-2^-104が、
• FMAでなくただの乗算+加算だと、1-2^-104が1に丸められ、それに(-1)を足すので0が、

• gcc4.1では、C99のfma関数のサポートそのものが無い。
• gcc4.4では、C99のfma関数が実装され、libm内にfma関数が存在する。-mfmaによるHaswellのFMA命令の有効化は実装されていない。libm内のfma関数は、CPUがFMA命令を持っていればそれを使い、持っていなければ単にfma(a,b,c)=a*b+cを実行する。
• gcc4.8では、-mfmaによるfma命令の有効化が実装された。これを付けてコンパイルすると、CPUがFMA命令を持つことが前提のコードが生成され、main関数中のfmaはlibm呼び出しではなく直接HaswellのFMA命令にコンパイルされる。
• gcc4.8で更に注目するべきは、libm呼び出しの場合IvyBridgeでも何故か正しい値が計算されたことである。これは普通ではなかなか出来ないはずで、CPUがFMA命令を持たなかった場合は(例えばmpfrを用いた)正確なFMA命令のエミュレーションを行うような実装が行われたと思われる。