gccの最適化
GCC is now transforming: ((x % CONSTANT) == 0)
— Marc B. Reynolds (@marc_b_reynolds) May 5, 2019
to mod inverse and using rotate when it's even.
Not in clang yet thoughhttps://t.co/nfywTIrTe0
ある数が定数で割り切れるかを%を使って判定するコードが
#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} \equiv x \mod m $$
となるような整数$x$、つまり$a$の逆数が存在します。 この$x$はオイラー関数 $\varphi(n)$(=$n$と互いに素である$1$以上$n$以下の自然数の数)を用いて求めることができます。
$$ a^{-1} \equiv a^{\varphi(m) - 1} $$
通常この$\varphi(n)$を求めるのは素因数分解などのめんどくさい計算が必要になりますが、 素数$p$のべき乗の場合は単純な式になります。
$$ \varphi(p^k) = p^k - p^{k-1} $$
$p=2$の場合は以下のようにもっと単純になります。
$$\begin{align} \varphi(2^k) &= 2^k - 2^{k-1} \\\ &= 2^{k-1}(2-1) \\\ &= 2^{k-1} \end{align}$$
べき乗は簡単な工夫で $\mathcal{O}(\log{n})$で求めることができます。 というわけで法$2^k$のときの奇数$a$の逆数(逆数が存在するのは$a$が$2^{k}$と互いに素のとき、 つまり奇数のときです)はそれなりに速く求めることができて、それは以下のような式になります。
$$ a^{-1} \equiv a^{2^{k-1}-1} \mod 2^k $$
オーバーフローは剰余算
ここでuint32_tの性質について考えてみます。
uint32_tは32ビット整数です。
32ビットで表せる範囲を超えた数値はオーバーフロー(アンダーフロー)を起こし、
32ビットで表せる範囲に戻ってきます。
つまり！uint32_tは！勝手に$2^{32}$を法とする剰余類環になっています！
ということはuint32_tの範囲で任意の奇数の逆数を求めることができて、uint32_tの範囲で除算ができます。
uint32_t inv = 3067833783; // 7の逆数
return 63 * inv; // 9
除算の性質
わかりやすく数字を小さくして考えてみましょう。 法$8 = 2^3$において$3$の逆元は$3$です。 0から7までの整数に3をかけてみます。
| x | x*3%8 |
|---|---|
| 0 | 0 |
| 1 | 3 |
| 2 | 6 |
| 3 | 1 |
| 4 | 4 |
| 5 | 7 |
| 6 | 2 |
| 7 | 5 |
結果に重複はなく、0から7までの数字が並べ替えられていることがわかります。 ここで右側の列をキーにしてソートし直してみましょう。
| x*3%8 | x |
|---|---|
| 0 | 0 |
| 1 | 3 |
| 2 | 6 |
| 3 | 1 |
| 4 | 4 |
| 5 | 7 |
| 6 | 2 |
| 7 | 5 |
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 * 2^m $$
まず奇数部分$n$の逆数を作って判定対象$k$を割ります。 あとは$k/n$が$2^m$で割り切れるかを判定すればいいです。
ある数が$2^m$で割り切れるということは、下位$m$ビットがゼロであるということです。 これは$2^m-1$とのANDがゼロかどうかで判定することもできますが、 右ローテート を使うと$n$で割り切れるかどうかの判定を一度に行うことができます。
$k$が$2^m$で割り切れる=下位$m$ビットがゼロであるときに$k/n$を$m$回右ローテートすると、 これは右シフト と同じであり値は単に$2^m$で割られます。 このとき$n$で割った余りが0になる値域を$2^m$分の1にした範囲に値があれば、$k$は$n$と$2^m$の両方で割り切れるため$c$で割り切れます。
そうでないときに$k/n$を$m$回右ローテートすると、下位ビットの1がぐるっとまわって上位ビットに付加されます。 すると値は非常に大きくなるため$n$で割った余りが0になる値域を$2^m$分の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ビット整数を$2^{32}$で割った余りと捉えてオーバーフローを活用するという発想がなかったので勉強になりました。 そのうちなにかに使えるかもしれない……