128bit除算も定数除算最適化したい

一般に整数除算(ここでは剰余演算も含むこととします)を行う命令は遅いです。 そのため、除算はなるべく回避したいものです。

除数がコンパイル時定数の場合、コンパイラの最適化により、除算が取り除かれることがあります。 これを定数除算最適化と呼びます。

定数除算最適化では、除算はいくつかの軽量な命令に分解されます。 二のべき乗で割る場合、単にシフト命令が使われます。 それ以外の数で割る場合、乗算とシフト命令が使われます。 符号付きの除算の場合、零への丸めの部分の取り扱いが面倒ですが、追加のシフト命令と加算命令だけで(分岐命令を使わずに)なんとかできるようです。

さて、この定数除算最適化は、128bit整数の除算の場合は常には行われないようです。 これを何とかしてみたいと思います。 もちろん、手で最適化されたコードを書けばできますが、コンパイル時計算により除数が決まる場合などに対応するのは困難という問題があるなど、柔軟性に欠けます。 そこで、128bit除算を行う、64bit除算のみを用いたC言語ソースコードを作ることにします。 こうすれば、マジックナンバーの導出はコンパイラにまかせることができます。

前提:コンパイラが最適化する場合としない場合

clangは、二のべき乗の時にシフトにする最適化だけ行うようです。 gccも、二のべき乗の時にシフトにする最適化を行います。 それに加えて、除数が小さい時は最適化してくれることが多いようです。 66以下の自然数で割る場合はすべて最適化してくれました。 67,83,101,……で割る場合は__udivti3 (128bit÷128bitを計算する外部ルーチン)の呼び出しになりました。 大きな数でも最適化してくれる数はいくつか存在するようです。 また、二のべき乗以外の時に最適化で出力されるコードは、25~40命令程度です。

実装

基本的な考え方は、64bit÷64bitを基本演算とした長除算を書けばいい、というものです。 これを実装したのが以下のコードになります。

void partial_divide( __uint128_t* R, __uint128_t* Q, unsigned long long D, int shamt ) {
    assert((unsigned long long)(*R >> shamt) == (*R >> shamt));
    *Q += (__uint128_t)((unsigned long long)(*R >> shamt) / D) << shamt;
    *R = *R - (*R >> shamt << shamt) + ((__uint128_t)((unsigned long long)(*R >> shamt) % D) << shamt);
}

void my_divrem( __uint128_t* R, __uint128_t* Q, __uint128_t X, unsigned D ) {
    *R = X;
    *Q = 0;
    partial_divide(R, Q, D, 64);
    partial_divide(R, Q, D, 32);
    partial_divide(R, Q, D, 0);
}

このコードは、Dが32bit整数しか受け付けないという制約がありますが、かなり高速に動作します。

性能測定

実験に使ったコードは以下です。

__uint128_t S = (__uint128_t)1 << 125; // ※グローバル変数にしないと最適化されてしまうので注意
__uint128_t E = S + 1000 * 1000 * 1000;

int main() {
    auto start = std::chrono::system_clock::now();
    __uint128_t sum = 0;
    for( __uint128_t t = S; t < E; ++t ) {
#ifdef G
        sum += t / DD;
#else
        __uint128_t R, Q;
        my_divrem(&R, &Q, t, DD);
        sum += Q;
#endif
    }
    auto end = std::chrono::system_clock::now();
    std::cout << (unsigned long long)sum << " " << std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() << std::endl;
}

私の書いたコードのコンパイルには、clang++12.0.1をおよびg++12.1.0を使いました。コンパイルオプションは-std=c++20 -O3 -march=nativeです。 コンパイルは、AMD EPYC 7702上で行いました(-march=nativeをつけても、以下で使うどのCPU上でも動くコードができました)。

速度の比較相手は、コンパイラの自動最適化が出力した機械語コードです。 これを作るためには、g++12.1.0を使いました。コンパイルオプションは-std=c++20 -O3 -march=native -DGです。

実行にかかった時間は以下のようになりました。

DD=3の場合↓

使ったCPU 自動最適化 my_divrem(clang) my_divrem(gcc)
AMD EPYC 7702 3.0s 3.5s 6.2s
Intel Xeon E5-2643 v3 3.1s 3.4s 5.9s
Intel Core i9-12900KF 1.3s 1.6s 1.9s

DD=67の場合↓

使ったCPU 自動最適化 my_divrem(clang) my_divrem(gcc)
AMD EPYC 7702 27s 3.7s 6.3s
Intel Xeon E5-2643 v3 33s 3.5s 6.3s
Intel Core i9-12900KF 3.9s 1.7s 1.9s

DD=3の場合、my_divrem(clang)は自動最適化よりやや遅い程度となりました。 「自動最適化」はさすがにコンパイラが出力したコードだけあって速いです。 とはいえ、私の強引な実装もそれほど遅くはないようです。 my_divrem(gcc)はかなり遅いですが、これはスピルが大量に発生してしまったのが原因です。 __uint128_tを取り扱うコードはレジスタ圧が高いため、レジスタ割り付けが難しいようです。

DD=67の場合、my_divremはclangもgccも両方とも自動最適化より速くなりました。 「自動最適化」は定数除算最適化をしていないため、除算命令を実行する必要があり低速です。 最近のIntelのCPUでは除算命令が高速化されたのか、除算命令にコンパイルされてもある程度の速さですが、依然として定数除算最適化をした場合よりも二倍以上遅いです。

もっと大きな数での剰余

Dが32bitよりも少し大きいくらいの場合であれば、partial_divideを呼ぶ回数を増やすことで対処可能です。 例えば、Dが42bitに収まる整数であれば、次にようにすればよいです。

void my_divrem( __uint128_t* R, __uint128_t* Q, __uint128_t X, unsigned long long D ) {
    *R = X;
    *Q = 0;
    partial_divide(R, Q, D, 66);
    partial_divide(R, Q, D, 44);
    partial_divide(R, Q, D, 22);
    partial_divide(R, Q, D, 0);
}

つまり、Dがkビットであれば、(64-k)ビットずつずらしながら除算を行えばよいです。 どこまでやればいいかというと、(64-k)*Nが64以上になるところまでです。

剰余環における演算がしたい場合

剰余環における乗算を実装するために128bit整数が必要となっている場合、さらに最適化することができます。 法をMとして、Mが42bitに収まる整数であれば、次のようにして問題ありません。

unsigned long long mul_with_mod( unsigned long long x, unsigned long long y, unsigned long long M ) { // x,yはM未満とする
    __uint128_t Q, R = (__uint128_t)x * y; // Rは84bitに収まる
    partial_divide(&R, &Q, M, 22); // 84-22=64なので64bit除算でOK
    partial_divide(&R, &Q, M, 0); // 上の除算の余りRは42+22=64で64bitに収まっている
    return R;
}

Mがkビットである時、(64-k)ビットずつずらしながら除算をするという点は先ほどと同じです。 一方、どこまでやるかが違っていて、(64-k)*Nが2k-64以上になるところまでです。

例えば、法Mが48bitに収まる整数であれば、partial_divideを三回に増やします。このとき、shamtは32, 16, 0となります。 法Mが51bitに収まる整数であれば、partial_divideを四回に増やします。このとき、shamtは39, 26, 13, 0となります。

まとめ

128bit整数を、32bitに収まる整数定数で割るプログラムを最適化しました。 具体的には、コンパイラが最適化しやすいようなコード(partial_divide)から構成される除算関数(my_divrem)を作りました。 その結果、gccが最適化できないような定数除算(例えば、67で割る)プログラムにおいて除算命令を回避し、大幅に高速化できることが分かりました。 gccが最適化する定数除算(例えば、3で割る)プログラムでも、速度低下はわずかでした。 定数除算最適化を行わないclang向けにコンパイルする場合に利用価値が高そうです。

また、32bitより少し大きいくらい(例えば、42bit)の整数で割る場合にも応用できることを示しました。 さらに、重要な応用として、剰余環における乗算の実装の場合、一般の場合より演算回数を少なくできることを示しました。