最終更新: 2024/4/27

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

柏木雅英

1. はじめに

kvライブラリで区間演算を使うとき、デフォルトで提供されているのは、

interval.hpp + rdouble.hpp で「doubleを両端に持つ区間演算」
interval.hpp + dd.hpp + rdd.hpp で「ddを両端に持つ区間演算」
interval.hpp + mpfr.hpp + rmpfr.hpp で「mpfrを両端に持つ区間演算」

の3通りである。ここで、前2者 (double, dd) については、その丸めモードの変え方が何通りもあり、それらはコンパイルオプションで制御するようになっている。それをここにまとめておく。

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

オプションなし

まず、何もオプションを付けない状態では、

Microsoft Visual C++ の場合、float.hで提供される_controlfp_sで、
それ以外の場合、fenv.hで提供されるfesetroundで、

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

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

Visual C++ 2012以前は、fenv.hを実装していなかった。
Visual C++ 2013, 2015ではfenv.hを実装したが、上向き丸めと下向き丸めの丸め方向が逆になるという致命的なバグがあった。
(Visual C++ 2017でようやく正常になった。)

-DKV_FASTROUND

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

また、fesetroundや_controlfp_sでは、

1. コントロールレジスタを読み込み
2. 読み込んだデータの丸めモードに対応する部分のbitを書き換える
3. コントロールレジスタに書き込み

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

1., 2.の部分はプログラム開始時に実行して値を記憶しておく
3. の部分は記憶しておいたデータを書き込むだけ

とすることによって高速化が期待できる。

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

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

-DKV_NOHWROUND

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

アルゴリズムの詳細は、 emu.pdf (2014/6/9版, 2024/4/27に微修正) を見て欲しい。

実装の詳細は

kv/rdouble-nohwround.hpp (doubleを両端に持つ区間演算の場合)
kv/rdd-nohwround.hpp (ddを両端に持つ区間演算の場合)

を見て欲しい。

[注意] -DKV_NOHWROUNDによる方向付き丸めのエミュレートは、 doubleの演算がIEEE754に準拠していることを前提に設計されている。多くの標準的な環境で動作するが、例えばIntel CPUの32bit環境はIEEE754に従わないことがあるので、本ライブラリはコンパイルできないようになっている (version 0.4.54から)。具体的には、#include <cfloat>した後マクロFLT_EVAL_METHODを調べ、それが0でないならばエラーでコンパイルできないようにした。 Intel CPUの32bit環境でどうしてもKV_NOHWROUNDを使いたい場合は、コンパイルオプションとして -msse2 -mfpmath=sse を付けるなどして SSE2を使わせるようにすると、上記マクロが0になってコンパイルできるようになる。このマクロが定義されないコンパイラ、あるいはSSE2を持たないCPUは、古すぎるためサポートされない。

-DKV_USE_AVX512

Intelのいくつかの最新のCPUでは、AVX-512という機能が搭載されている。2017年11月現在では、

Xeon Phi Knights Landing
Skylake-X (Core i9 7900Xなど)
Skylake-SP (Xeon Platinum/Gold/Silver/Bronze など)

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

実装の詳細は

kv/rdouble-avx512.hpp (doubleを両端に持つ区間演算の場合)
kv/rdd-avx512.hpp (ddを両端に持つ区間演算の場合)

を見て欲しい。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、コンパイラなら大丈夫と思うが、以下の手順でチェックして欲しい。また不具合があれば報告して欲しい。

test/test-fma.cc をコンパイルしてみる。コンパイル出来なければfma命令は使えない。
test/test-fma.cc の実行結果が "fma works!" なら問題ないが、 "fma is broken!"ならそのシステムのfma命令はfmaとしての精度を持っていないので、使うことは出来ない。
適当なddを使ったプログラムについて、-DKV_USE_TPFMA を付けた場合と付けない場合で速度比較をしてみる。fmaがエミュレーションで実行されるシステムがあり、この場合異常な速度低下が見られるはず。

5. ARM CPUへの対応

オプションなし

fenv.hで提供されるfesetroundが使われる。

-DKV_FASTROUND

version 0.4.55で、ARM 32bit/64bitにおいてこのオプションが使用可能になった。 Inline Assemblerを用いて高速な丸め変更を行う。

-DKV_NOHWROUND

ハードウェアによる丸めモードの変更を使わずにtwosum, twoproductを用いて方向付き丸めのエミュレートを行う。

6. ベンチマーク

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

使用した計算環境は、

Intel Core i5 11400 2.6GHz/4.4GHz, Memory 32G
Ubuntu 22.04 LTS
gcc 11.2.0

である。

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

計算精度(端点の型) コンパイルオプション
(-O3 -DNDEBUGに加えて) 計算時間

-DKV_USE_TPFMAなし -DKV_USE_TPFMAあり

double なし 6.65 sec

-DKV_FASTROUND 3.74 sec

-DKV_NOHWROUND 7.88 sec 6.05 sec

-DKV_USE_AVX512 -mavx512f 2.74 sec

dd なし 37.20 sec 33.58 sec

-DKV_FASTROUND 25.08 sec 21.52 sec

-DKV_NOHWROUND 92.24 sec 71.81 sec

-DKV_USE_AVX512 -mavx512f 16.47 sec

計算精度(端点の型)	コンパイルオプション (-O3 -DNDEBUGに加えて)	計算時間
-DKV_USE_TPFMAなし	-DKV_USE_TPFMAあり
double	なし	6.65 sec
-DKV_FASTROUND	3.74 sec
-DKV_NOHWROUND	7.88 sec	6.05 sec
-DKV_USE_AVX512 -mavx512f	2.74 sec
dd	なし	37.20 sec	33.58 sec
-DKV_FASTROUND	25.08 sec	21.52 sec
-DKV_NOHWROUND	92.24 sec	71.81 sec
-DKV_USE_AVX512 -mavx512f	16.47 sec

続いて、M1 MacBook Airにおけるベンチマーク結果を記載する。現状ではARMのCPUではハードウェアによる丸めモード変更の実行時間におけるペナルティが大きいようで、エミュレーションが一番速いという残念な結果になっている。

使用した計算環境は、

Apple M1 MacBook Air 3.2GHz, Memory 16G
macOS Monterey 12.4
gcc 12.2.0

である。

計算精度(端点の型) コンパイルオプション
(-O3 -DNDEBUGに加えて) 計算時間

-DKV_USE_TPFMAなし -DKV_USE_TPFMAあり

double なし 24.45 sec

-DKV_FASTROUND 22.36 sec

-DKV_NOHWROUND 9.08 sec 6.03 sec

dd なし 108.8 sec 100.6 sec

-DKV_FASTROUND 102.4 sec 94.64 sec

-DKV_NOHWROUND 84.60 sec 53.76 sec