KOTET'S PERSONAL BLOG

gccの最適化

ある数が定数で割り切れるかを%を使って判定するコードが

#include "stdint.h"

uint32_t is_div_7(uint32_t k)
{
  return (k % 7)==0;
}

こんな感じに最適化されるようです。

is_div_7(unsigned int):
        imul    edi, edi, -1227133513
        xor     eax, eax
        cmp     edi, 613566756
        setbe   al
        ret

わかりやすくCに再翻訳するとこんな感じでしょうか。

#include "stdint.h"

uint32_t is_div_7(uint32_t k)
{
  return (k * -1227133513) <= 613566756;
}

剰余演算(x86では除算と同じ処理です)がなくなり、乗算と比較だけになっているのがわかります。 除算は他の演算と比べてコストが大きいので、なくすことができると嬉しいです。 しかし-1227133513や613566756などという数字はどこから出てきたのでしょうか？

モジュラ逆数

自然数mで割ったときの余りが同じ整数を同一視することを「法mで考える」といいます。 このとき、mと互いに素な整数aに対して

a−1≡xmodm

となるような整数x、つまりaの逆数が存在します。 このxはオイラー関数 φ(n)(=nと互いに素である1以上n以下の自然数の数)を用いて求めることができます。

a−1≡aφ(m)−1

通常このφ(n)を求めるのは素因数分解などのめんどくさい計算が必要になりますが、 素数pのべき乗の場合は単純な式になります。

φ(pk)=pk−pk−1

p=2の場合は以下のようにもっと単純になります。

φ(2k)=2k−2k−1 =2k−1(2−1) =2k−1

べき乗は簡単な工夫で O(logn)で求めることができます。 というわけで法2kのときの奇数aの逆数(逆数が存在するのはaが2kと互いに素のとき、 つまり奇数のときです)はそれなりに速く求めることができて、それは以下のような式になります。

a−1≡a2k−1−1mod2k

オーバーフローは剰余算

ここでuint32_tの性質について考えてみます。 uint32_tは32ビット整数です。 32ビットで表せる範囲を超えた数値はオーバーフロー(アンダーフロー)を起こし、 32ビットで表せる範囲に戻ってきます。

つまり！uint32_tは！勝手に232を法とする剰余類環になっています！

ということはuint32_tの範囲で任意の奇数の逆数を求めることができて、uint32_tの範囲で除算ができます。

uint32_t inv = 3067833783; // 7の逆数
return 63 * inv; // 9

除算の性質

わかりやすく数字を小さくして考えてみましょう。 法8=23において3の逆元は3です。 0から7までの整数に3をかけてみます。

結果に重複はなく、0から7までの数字が並べ替えられていることがわかります。 ここで右側の列をキーにしてソートし直してみましょう。

3で割ったあまりが0のもの、1のもの、2のものと纏まっていることがわかるでしょうか。 法である8を3で割った数(切り上げ)をcとおくと、余りがkとなる数の結果の値域は[ck,c(k+1))になります。 つまり結果がどの範囲にあるか見れば余りがいくつになるかわかるということです。

完成

これらをすべて組み合わせると最初のコードになります。 任意の奇数で同じ手順での最適化が可能です。

is_div_7(unsigned int):
        imul    edi, edi, -1227133513 ; == 3067833783 == 7^(2^31 - 1) mod 2^32 == 7の逆数
        xor     eax, eax
        cmp     edi, 613566756 ; 2^32 / 7
        setbe   al ; (k * 7^(2^31 - 1) mod 2^32) <= 2^32 / 7
        ret

偶数は？

偶数も少し工夫すると簡単に最適化できます。 まずは実際に最適化されたコードを見てみましょう。

uint32_t is_div_14(uint32_t k)
{
  return (k % 14)==0;
}

is_div_14(unsigned int):
        imul    edi, edi, -1227133513
        xor     eax, eax
        ror     edi
        cmp     edi, 306783378
        setbe   al
        ret

偶数cは奇数nと指数mを使って以下のように表せます。

c=n∗2m

まず奇数部分nの逆数を作って判定対象kを割ります。 あとはk/nが2mで割り切れるかを判定すればいいです。

ある数が2mで割り切れるということは、下位mビットがゼロであるということです。 これは2m−1とのANDがゼロかどうかで判定することもできますが、 右ローテート を使うとnで割り切れるかどうかの判定を一度に行うことができます。

kが2mで割り切れる=下位mビットがゼロであるときにk/nをm回右ローテートすると、 これは右シフト と同じであり値は単に2mで割られます。 このときnで割った余りが0になる値域を2m分の1にした範囲に値があれば、kはnと2mの両方で割り切れるためcで割り切れます。

そうでないときにk/nをm回右ローテートすると、下位ビットの1がぐるっとまわって上位ビットに付加されます。 すると値は非常に大きくなるためnで割った余りが0になる値域を2m分の1にした範囲から外れます。

以上のようにして偶数で割り切れるかどうかの判定を除算なしに行うことができました。

is_div_14(unsigned int):
        imul    edi, edi, -1227133513 ; == 3067833783 == 7^(2^31 - 1) mod 2^32 == 7の逆数
        xor     eax, eax
        ror     edi ; 1ビット右ローテート
        cmp     edi, 306783378 ; (2^32 / 7) / 2
        setbe   al
        ret

感想

符号なし32ビット整数を232で割った余りと捉えてオーバーフローを活用するという発想がなかったので勉強になりました。 そのうちなにかに使えるかもしれない……

参考

• Optimizing is-multiple checks with modular arithmetic
• モジュラ逆数 - Wikipedia
• 演算子強度低減 - Wikipedia