最終更新: 2017/11/8

丸めモードの変え方とコンパイルオプションまとめ

柏木 雅英

1. はじめに

kvライブラリで区間演算を使うとき、デフォルトで提供されているのは、 の3通りである。ここで、前2者 (double, dd) については、 その丸めモードの変え方が何通りもあり、それらはコンパイルオプションで制御するようになっている。 それをここにまとめておく。

2. 丸めモードの変え方とコンパイルオプション

オプションなし

まず、何もオプションを付けない状態では、 丸めモードを変更している。実装の詳細はkv/hwround.hppを見て欲しい。 コンパイラがC99準拠を謳っている場合、fenv.hが使えることになっているので、多くの環境でこのまま動作するはずである。

なお、Visual C++の場合のみ実装が異なるのは、以下の理由による。

-DKV_FASTROUND

IntelのCPUはFPUとSSE2の2種類の演算器を持っており、また丸めモード 指定のためのレジスタは個別に持っているため、上のfesetroundや _controlfp_sでは両方のレジスタの丸めモードを変更している。 ところが、一般的にIntelの64bitモードの場合、SSE2の演算しか使わないため、 SSE2のレジスタのみを変更するようにすると高速化が期待できる。

また、fesetroundや_controlfp_sでは、

という動作を行っている。ここで、コントロールレジスタの丸めモードに関する 部分以外はプログラム実行中に変化しないと仮定すると、 とすることによって高速化が期待できる。

-DKV_FASTROUNDを付けると、この2つを行なうことによって高速化する。ただし、 IntelのCPUで64bitモードで動作している場合しか使うことは出来ない

実装の詳細はkv/hwround.hppを見て欲しい。

-DKV_NOHWROUND

丸めモードの変更を使わずにtwosum, twoproductを用いて 方向付き丸めのエミュレートを行うようになる。 ハードウェアによる丸めモードの変更一切行わずデフォルトの 最近点丸めモードのみで動作するので、非常に幅広い環境で動くことが期待される。 速度はあまり速くないはず。

アルゴリズムの詳細は、 emu.pdf (2014/6/9版) を見て欲しい。

実装の詳細は

を見て欲しい。

[注意] -DKV_NOHWROUNDによる方向付き丸めのエミュレートは、 doubleの演算がIEEE754に 準拠していることを前提 に設計されているため、多くの標準的な環境で動作するが、 IntelのCPUの32bitモードでは正しく動作しない可能性がある。 IntelのCPUの32bitモードはFPUによって内部80bitで計算されてしまうため、 IEEE754に準拠した計算結果にならないことがあるためである。

-DKV_USE_AVX512

Intelのいくつかの最新のCPUでは、AVX-512という 機能が搭載されている。2017年11月現在では、 に搭載されている。512bitのレジスタによるベクトル化の機能だが、 このAVX-512では、 命令そのものに丸めモード指定を含んでいるような加減乗除と平方根 が追加された。これを使うと、丸めモードの変更のオーバーヘッドなしに区間演算が 実現できる。 当然、AVX-512を持つ最新のCPUでしか動作しない。 また、コンパイル時にAVX-512命令を生成させるにはコンパイルオプションが 必要なことも多く、 2017年11月現在、gccでは -mavx512f が必要

実装の詳細は

を見て欲しい。intrinsic命令を使って実装している。本来8並列の能力を 持つAVX-512だが、その機能は使っていない。

test-rounding.ccによる動作チェック

これらの、オプション無し、-DKV_FASTROUND, -DKV_NOHWROUND, -DKV_USE_AVX512 について、test/test-rounding.cc をそのオプションでコンパイル、実行すると、 正しく動作しているかどうか簡易チェックできる。 正しく丸めの向きが変更されていれば、
add_up: OK
add_down: OK
sub_up: OK
sub_down: OK
mul_up: OK
mul_down: OK
div_up: OK
div_down: OK
sqrt_up: OK
sqrt_down: OK
のように表示される。

3. 丸めモードの変え方の実装コード

"add_up" (加算の上向き丸め) を例にとって、 それぞれのコンパイルオプションの場合の実装コードのイメージを記しておく。

オプションなし

#include <fenv.h>

static double add_up(const double& x, const double& y) {
    volatile double r, x1 = x, y1 = y;
    fesetround(FE_UPWARD);
    r = x1 + y1;
    fesetround(FE_TONEAREST);
    return r;
}

-DKV_FASTROUND

#include <emmintrin.h>

unsigned long int reg = _mm_getcsr();
unsigned long int nearest = (reg & ~0x00006000) | 0x00000000;
unsigned long int down = (reg & ~0x00006000) | 0x00002000;
unsigned long int up =   (reg & ~0x00006000) | 0x00004000;
unsigned long int chop = (reg & ~0x00006000) | 0x00006000;

static double add_up(const double& x, const double& y) {
    volatile double r, x1 = x, y1 = y;
    _mm_setcsr(up);
    r = x1 + y1;
    _mm_setcsr(nearest);
    return r;
}

-DKV_NOHWROUND

#include <cmath>
#include <limits>

static void twosum(const double& a, const double& b, double& x, double& y) {
    double tmp;

    x = a + b;
    if (std::fabs(a) > std::fabs(b)) {
        tmp = x - a;
        y = b - tmp;
    } else {
        tmp = x - b;
        y = a - tmp;
    }
}

static double succ(const double& x) {
    static const double th1 = std::ldexp(1., -969);
    static const double th2 = std::ldexp(1., -1021);
    static const double c1 = std::ldexp(1., -53) + std::ldexp(1., -105);
    static const double c2 = std::ldexp(1., -1074);
    static const double c3 = std::ldexp(1., 53);
    static const double c4 = std::ldexp(1., -53);

    double a, c, e;

    a = std::fabs(x);
    if (a >= th1) return x + a * c1;
    if (a < th2) return x + c2;
    c = c3 * x;
    e = c1 * std::fabs(c);
    return (c + e) * c4;
}

static double add_up(const double& x, const double& y) {
    double r, r2;

    twosum(x, y, r, r2);
    if (r == std::numeric_limits<double>::infinity()) {
        return r;
    } else if (r == -std::numeric_limits<double>::infinity()) {
        if (x == -std::numeric_limits≤double>::infinity() || y == -std::numeric_limits::infinity()) {
            return r;
        } else {
            return -(std::numeric_limits≤double>::max)();
        }
    }

    if (r2 > 0.) {
        return succ(r);
    }

    return r;
}

-DKV_USE_AVX512

#include <immintrin.h>

static double add_up(const double& x, const double& y) {
    volatile __m128d x1, y1, z1;
    double r;
    x1 = _mm_set_sd(x);
    y1 = _mm_set_sd(y);
    z1 = _mm_add_round_sd(x1, y1, _MM_FROUND_TO_POS_INF | _MM_FROUND_NO_EXC);
    _mm_store_sd(&r, z1);
    return r;
}

4. twoproductでのFMA命令の使用

方向付き丸めのエミュレートやdouble-double計算ではtwoproductと呼ばれる アルゴリズムが使われているが、このアルゴリズムはFMA命令を使うと 高速化できることが知られている。 -DKV_USE_TPFMAオプションを使うと、twoproductを FMA命令を用いた高速な実装に置き換えることが出来る。 実装コードは以下の通り。

-DKV_USE_TPFMAなし

static void split(const double& a, double& x, double& y) {
    static const double sigma = (double)((1L << 27) + 1);
    double tmp;

    tmp = a * sigma;
    x = tmp - (tmp - a);
    y = a - x;
}

static void twoproduct(const double& a, const double& b, double& x, double& y) {
    static const double th = std::ldexp(1., 996);
    static const double c1 = std::ldexp(1., -28);
    static const double c2 = std::ldexp(1., 28);
    static const double th2 = std::ldexp(1., 1023);

    double na, nb, a1, a2, b1, b2;

    x = a * b;
    if (std::fabs(a) > th) {
        na = a * c1;
        nb = b * c2;
    } else if (std::fabs(b) > th) {
        na = a * c2;
        nb = b * c1;
    } else {
        na = a;
        nb = b;
    }
    split(na, a1, a2);
    split(nb, b1, b2);
    if (std::fabs(x) > th2) {
        y = ((((a1 * 0.5) * b1 - (x * 0.5)) * 2 + a2 * b1) + a1 * b2) + a2 * b2;
    } else {
        y = (((a1 * b1 - x) + a2 * b1) + a1 * b2) + a2 * b2;
    }
}

-DKV_USE_TPFMAあり

static void twoproduct(const double& a, const double& b, double& x, double& y) {
    x = a * b;
    y = fma(a, b, -x);
}
FMA命令の使用については、IntelのHaswell以降のCPUとあまり古くないOS、 コンパイラなら大丈夫と思うが、以下の手順でチェックして欲しい。 また不具合があれば報告して欲しい。

5. ベンチマーク

interval<double>, interval<dd>を多用したプログラムとして、 example/test-nishi.cc を題材として、上記4種のコンパイルオプションにおける 計算速度を調べてみた。 なおtest-nishi.ccは、-DTEST_DDで全ての計算がddで行われる。

使用した計算環境は、

Intel Core i9 7900X 3.3GHz/4.3GHz, Memory 128G
Ubuntu 16.04 LTS
gcc 5.4.0
である。

なお、前節で述べたように -DKV_USE_TPFMA オプションを使うと 計算速度が改善する場合があり、その場合は両者の計算時間を記した。

計算精度(端点の型) コンパイルオプション 計算時間
-DKV_USE_TPFMAなし -DKV_USE_TPFMAあり
double -O3 10.81 sec
-O3 -DKV_FASTROUND 8.40 sec
-O3 -DKV_NOHWROUND 19.91 sec 12.68 sec
-O3 -DKV_USE_AVX512 -mavx512f 5.55 sec
dd -O3 54.10 sec 44.94 sec
-O3 -DKV_FASTROUND 46.08 sec 36.72 sec
-O3 -DKV_NOHWROUND 185.7 sec 118.9 sec
-O3 -DKV_USE_AVX512 -mavx512f 21.42 sec

計算時間は3回計測した平均を記した。