kashiの日記

• 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命令のエミュレーションを行うような実装が行われたと思われる。

• テンプレートを排除しました。従って内部型はdoubleのみですが、プログラムが読みやすくなっています。
• テンプレートを排除したおかげでboostに依存しません。単体で使えます。
• 必要なファイルはinterval.hpp一つです。
• 内部型のカスタマイズが出来ないこと以外は、kvライブラリに含まれるinterval.hppと機能は同等です。数学関数も多数あります。
• 丸めの変更は単純にfesetroundを使っています。
• 文字列との方向丸め付き相互変換機能もあります。従って、文字列からその文字列の表す数値を含むように区間を初期化できるし、表示も真の値を含むように外側に丸めるようになっています。
• volatileを使った最適化対策はきちんとやってます。

#include <kv/interval.hpp>
#include <kv/mpfr.hpp>
#include <kv/rmpfr.hpp>

int main()
{
        int i;
        kv::interval< kv::mpfr<17000> > x, y, z;
        std::cout.precision(17);

        x = 1.;
        y = 1.;

        for (i=2; i<=10000; i++) {
                z = (1 + 2 * y) / (x * y * y);
                std::cout << i << " " << z << "\n";
                x = y;
                y = z;
        }
}

2 [3,3]
3 [0.77777777777777777,0.77777777777777778]
4 [1.4081632653061224,1.4081632653061225]
5 [2.4744801512287334,2.4744801512287335]
(中略)
9997 [1.0388807450632487,1.0388807450632488]
9998 [2.9584219778650379,2.958421977865038]
9999 [0.76071510663969069,0.7607151066396907]
10000 [1.4727965250386843,1.4727965250386844]

#include <iostream>
#include "lilib.h"

int main(){
        lilib::setPrecision(5000);

        LongInterval x, y, z;
        int i;

        x = 1.;
        y = 1.;

        for (i=2; i<=10000; i++) {
                z = (1 + 2 * y) / (x * y * y);
                std::cout << i << " " << z << "\n";
                x = y;
                y = z;
        }
}

2 < 3.000000, 2.177677e-5009>
3 < 7.777778e-1, 2.177677e-5009>
4 < 1.408163, 1.328045e-5008>
5 < 2.474480, 7.583944e-5008>
(中略)
9997 < 1.038881, 3.292724e-167>
9998 < 2.958422, 2.823943e-166>
9999 < 7.607151e-1, 2.314535e-166>
10000 < 1.472797, 1.307194e-165>