LoginSignup
21
9

More than 1 year has passed since last update.

@k_nitadori

LLVMのstd:: uniform_real_distributionが遅い件

Last updated at Posted at 2023-02-05

発端

一年ほど前の@kaityo256さんの発言：

MT19937がコンパイル条件によっては極端に遅いという話だったが、調べようと思ったまま1年近くが経ってしまった。別件で調べていてわかったのは、

  • 遅いのはstd::mt19937ではなくstd:: uniform_real_distributionの方（まぁモンテカルロやるなら確実に踏む）
  • LLVMとlibstdc++のときにAarch64であれx86_64であれ発動する、libc++にすることで回避できる

ということ。libstdc++とlibc++については

を参照してください。

計測

とりあえずこんなソース：

#include <iostream>
#include <iomanip>
#include <random>
#include <chrono>

int main(){
        std::mt19937_64 gen(334);
        std::uniform_real_distribution<double> uni(-0.5, 0.5);

        auto lrand = [&gen](){
                return gen();
        };
        auto drand = [&gen, &uni](){
                return uni(gen);
        };
        auto inv = [&gen](){
                return 1.0 / gen();
        };

        auto benchmark = [](auto kernel, long count){
                auto beg = std::chrono::system_clock::now();

                decltype(kernel()) sum = 0;
                for(long i=0; i<count; i++){
                        sum += kernel();
                }

                auto end = std::chrono::system_clock::now();

                std::chrono::duration<double> elapse = end - beg;

                std::cout << "elapse: " << elapse.count() << " sec" << std::endl;
                std::cout  << 1.e9 * elapse.count() / count << " nsec/call" << std::endl;
                std::cout << std::hex << std::hexfloat;
                std::cout << "sum : " << sum << std::endl;
                std::cout << std::dec << std::defaultfloat;
        };

        puts("MT");
        benchmark(lrand, 100'000'000);
        puts("U01");
        benchmark(drand, 100'000'000);
        puts("inv");
        benchmark(inv, 100'000'000);

        return 0;
}

型変換や除算が遅いのかと疑った名残でinv()というカーネルが残ってる。

とりあえずM1 MacのUbuntu 22 (UTM)：

MT
elapse: 0.277519 sec
2.77519 nsec/call
sum : 4b8c1836f9f50d02
U01
elapse: 41.4284 sec
414.284 nsec/call
sum : 0x1.817dd4e4a651fp+10
inv
elapse: 0.273696 sec
2.73696 nsec/call
sum : 0x1.a767ba2e0358ap-34

100倍以上遅い。clang++-14 -O3のときの結果だが、-ffast-=mathを付けると発動しなかった。
わざわざ仮想マシンのLinuxでやったのだが、Macにbrewで入れたLLVMではそもそもlibstdc++は使われずにこれは発動しないらしい。

SkylakeのXeonマシン：

MT
elapse: 0.342702 sec
3.42702 nsec/call
sum : 4b8c1836f9f50d02
U01
elapse: 6.33091 sec
63.3091 nsec/call
sum : 0x1.817dd4e4a651fp+10
inv
elapse: 0.380439 sec
3.80439 nsec/call
sum : 0x1.a767ba2e0358ap-34

こちらは20倍程度のオーバーヘッドであるが、同様に遅い。こちらも-ffast-mathで回避はできた。
なおこういったfpの危ないオプションを使わなければx86でもARMでもビットレベルで結果が一致している。

原因

おそらく原因となったソースコード（-Eを付けてコンパイルしたものを軽く整形）：

template<typename _RealType, size_t __bits, typename _UniformRandomNumberGenerator>
_RealType generate_canonical(_UniformRandomNumberGenerator& __urng)
{
        static_assert(std::is_floating_point<_RealType>::value, "template argument must be a floating point type");

        const size_t __b = std::min(static_cast<size_t>(std::numeric_limits<_RealType>::digits), __bits);
        const long double __r = static_cast<long double>(__urng.max()) - static_cast<long double>(__urng.min()) + 1.0L;
        const size_t __log2r = std::log(__r) / std::log(2.0L);
        const size_t __m = std::max<size_t>(1UL, (__b + __log2r - 1UL) / __log2r);
        _RealType __ret;
        _RealType __sum = _RealType(0);
        _RealType __tmp = _RealType(1);
        for (size_t __k = __m; __k != 0; --__k)
        {
                __sum += _RealType(__urng() - __urng.min()) * __tmp;
                __tmp *= __r;
        }
        __ret = __sum / __tmp;
        if (__builtin_expect(__ret >= _RealType(1), 0))
        {
                __ret = std::nextafter(_RealType(1), _RealType(0));
        }
        return __ret;
}

ここで先ずは$[0,1)$の疑似乱数を所与の実数型に対して作っている。
53-bitの倍精度に対して作るなら、32-bitのmt19937なら2回、64-bitのmt19937_64なら1回呼べばいいよというのがsize_t __mという変数の意味。その1とか2を計算するために、long doubleで

        const long double __r = static_cast<long double>(__urng.max()) - static_cast<long double>(__urng.min()) + 1.0L;
        const size_t __log2r = std::log(__r) / std::log(2.0L);

を計算しているというのが問題の箇所。ここはもちろんコンパイル時の定数にできる箇所なので（とはいってもconstexprのような厳格さはない）、GNUでは定数最適化されるコードがLLVMでは違ったように解釈されてしまっていた模様。
x86では次のようなコードになっている：

        flds    .LCPI0_0(%rip)
        fstpt   (%rsp)
        callq   logl
        fstpt   32(%rsp)                        # 10-byte Folded Spill
        flds    .LCPI0_1(%rip)
        fstpt   (%rsp)
        callq   logl
        fldt    32(%rsp)                        # 10-byte Folded Reload
        fdivp   %st, %st(1)

型変換の部分が長くて複雑だったので省いたが、x87でのlog, log, divが発生していることがわかる。loglについては関数コールになっているが、x87には $\mathtt{FYL2X}(x,y):=y \log_2 x$ のような計算をしてくれる命令があるのでハードウェアで処理されることにはなるが、それでも自然対数を2回計算して除算をするのは勿体ない、、、
Aarch64の場合：

        ldr     q0, [x8, :lo12:.LCPI1_0]
        bl      logl
        adrp    x8, .LCPI1_1
        str     q0, [sp]                        // 16-byte Folded Spill
        ldr     q0, [x8, :lo12:.LCPI1_1]
        bl      logl
        mov     v1.16b, v0.16b
        ldr     q0, [sp]                        // 16-byte Folded Reload
        bl      __divtf3
        bl      __fixunstfdi

こちらではlong doubleは128-bitの四倍精度となっているらしい。四倍精度のlog, log, divとunsigned longへの型変換はソフトウェアエミュレーションになっている。この結果80-bitをハードウェアで処理できたx86に比べて被害が大きな（そのかわり発見されやすい）ものとなっている。

対処法

当初の問題は富士通コンパイラのclang modeで発現したのだった。-Nclang -Kfastでも踏んでしまう。こちらはベンダがlibc++も準備してくれているのでそちらを使うのがいいだろう。ただしtrad modeで生成したオブジェクトとリンクできなくなるといった制限があったかもしれない1。Trad modeを使うのはC言語で書いたカーネル部分にとどめるか、C++を使うにしても標準ライブラリを用いないカーネル部分にとどめるのがよいだろう。

Linux環境だとUbuntuではlibc++を入れることができた（といっても手動でシンボリックリンクを張る必要があった）がRHELだと簡単には入りそうにない。見苦しいが特殊化してしまうというのは？

#include <random>
template<>
double std::generate_canonical<double, 53, std::mt19937_64>(std::mt19937_64& __urng)
{
        using _RealType = double;

        const long double __r = static_cast<long double>(__urng.max()) - static_cast<long double>(__urng.min()) + 1.0L;
        const size_t __m = 1;
        _RealType __ret;
        _RealType __sum = _RealType(0);
        _RealType __tmp = _RealType(1);
        for (size_t __k = __m; __k != 0; --__k)
        {
                __sum += _RealType(__urng() - __urng.min()) * __tmp;
                __tmp *= __r;
        }
        __ret = __sum / __tmp;
        if (__builtin_expect(__ret >= _RealType(1), 0))
        {
                __ret = std::nextafter(_RealType(1), _RealType(0));
        }
        return __ret;
}

__m = 1、これを計算するのに莫大な命令数を消費していた。
__rは$2^{64}$なので除算もlong doubleの乗算も発生していない。

まとめ

原因は判明してみれば富士通マターというよりはLLVMマターだった。とはいえモンテカルロ法等で利用頻度の高い関数なので使う人は注意して欲しい。
原因をまとめると：

  • GNUのソースコードが厳格なコンパイル時評価を要求していなかったがGCCなら問題なかった
  • long doubleの演算が多少発生してもx86では大きな被害とはならなかったがAarch64では四倍精度ライブラリの呼び出しとなってしまい非常に遅くなった

というもの。ライブラリにperformance portabilityを期待するのは難しく、根本的な解決をめざすのならLLVMのclang++を使うときはlibc++に移行するというのが正攻法なのだろう。

  1. 調べたらtradと結合するときこそlibc++でなければならないとのことでした。こっちをデフォルトにしてくれてもいい気がする、、、 ↩

21
9
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
21
9