LLVMのstd:: uniform_real_distributionが遅い件


 # 発端

一年ほど前の@kaityo256さんの発言：
<blockquote class="twitter-tweet"><p lang="ja" dir="ltr">うーん、A64fx、やっぱりモンテカルロは苦手っぽいなぁ。</p>&mdash; ロボ太 (@kaityo256) <a href="https://twitter.com/kaityo256/status/1489589319887372293?ref_src=twsrc%5Etfw">February 4, 2022</a></blockquote> <script async src="https://platform.twitter.com/widgets.js" charset="utf-8"></script>

MT19937がコンパイル条件によっては極端に遅いという話だったが、調べようと思ったまま1年近くが経ってしまった。別件で調べていてわかったのは、
 * 遅いのは`std::mt19937`ではなく`std:: uniform_real_distribution`の方（まぁモンテカルロやるなら確実に踏む）
 * LLVMとlibstdc++のときにAarch64であれx86_64であれ発動する、libc++にすることで回避できる

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

https://qiita.com/kojiohta/items/fb6c307365d1db388acc

を参照してください。

# 計測

とりあえずこんなソース：
```c++
#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`を付けてコンパイルしたものを軽く整形）：
```C++
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`で
```c++
        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++を使うにしても標準ライブラリを用いないカーネル部分にとどめるのがよいだろう。

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

Linux環境だとUbuntuでは`libc++`を入れることができた（といっても手動でシンボリックリンクを張る必要があった）がRHELだと簡単には入りそうにない。見苦しいが特殊化してしまうというのは？
```c++
#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++`に移行するというのが正攻法なのだろう。