2で割ることと3で割ること

この記事でお題にするのはCPUレジスタ上の整数除算です。以下、単に除算とも書きます。

除算は非常に高コストな演算なため、コンパイラは最適化によって、できるだけ整数除算を別の計算に置き換えようとします。

最適化ができる場合の一つとして、割る数が定数である場合があります。頭のいいコンパイラは、除算を乗算とビットシフト等を駆使した演算に置き換えます。この記事では、そういった最適化の背景にある理屈を部分的に解説します。

計算機環境としてはモダンなx86 CPUを仮定します。したがってレジスタは32/64ビットであり、負数は2の補数表現になっています。ある程度は他の命令セットでも通用する話になっているかもしれません。

そもそも整数の除算とは

プログラミングにおける整数の除算の定義について確認します。整数$n$を整数$d$で割るとき

$$
n = q \times d + r
$$

が成り立つように除算$n / d$を定めます。商$q$と余り$r$の選び方は明らか……に見えて負の数が絡むと曖昧で、プログラミング言語の仕様によります。

いくつかの言語を確認してみます。

C99

#include <stdio.h>
int main(void) {
	printf("%d %d\n",  5/ 2,  5% 2); // 2 1
	printf("%d %d\n",  5/-2,  5%-2); // -2 1
	printf("%d %d\n", -5/ 2, -5% 2); // -2 1
	printf("%d %d\n", -5/-2, -5%-2); // 2 -1
}

※ C99より前のC（C89など）では、この結果通りになるかは処理系依存でした。

Ruby

puts  5 /  2,  5 %  2 # 2 1
puts  5 / -2,  5 % -2 # -3 -1
puts -5 /  2, -5 %  2 # -3 1
puts -5 / -2, -5 % -2 # 2 -1

Clojure

(println (quot  5  2) (rem  5  2)) ; 2 1
(println (quot  5 -2) (rem  5 -2)) ; -2 1
(println (quot -5  2) (rem -5  2)) ; -2 -1
(println (quot -5 -2) (rem -5 -2)) ; 2 -1

有理数の範囲では$-5/2 = -2.5$ですが、商をどちら方向に丸めるかで差が生じている事がわかります。
Rubyは商が小さくなる方向に丸めるようです。一方、CやClojureで採用されているルールは、商の絶対値が小さくなる方向（0方向）への丸めです。

後者を床関数・天井関数を使った数式で書くとこんな感じになります。

$$
n / d = \begin{cases}
\lfloor \frac{n}{d} \rfloor & (x \geq 0) \
\lceil \frac{n}{d} \rceil & (x \lt 0)
\end{cases}
$$

多くのプログラミング言語や命令セット（x86含む）で採用されている除算はこれが定義になります。本記事でも、以降は0方向へ丸める除算のみを扱います。

除算は「ものすごく」遅い

CPU上における（レジスタ上の）演算は、その計算結果が得られるまで何クロックかかるかというレイテンシの観点で命令の「重さ」を確認できます。
例えば現代的なx86 CPUにおいての64ビット整数演算だと、加算は1クロックです。シフト演算も1クロックです。乗算は4クロックぐらいです。

じゃあ除算はというと、Intel CPUでいうSkylake世代では、なんと最大 97クロック かかります¹。メモリアクセスの方が速くなりかねないぐらいの遅さです。

つまり、四則演算などよく使われる演算の中で、除算は圧倒的に遅いのです。

コンパイラは除算命令を回避しようとする（定数２で割る場合）

それほどまでに遅い除算なので、コンパイラは極力除算命令の生成を避けようとします。

有名なのは2の累乗で割る場合でしょうか。Clangに-S -O2 -masm=intelのようにオプションを渡すことでアセンブリを生成し、最適化したらどうなるかを教えてもらいます。

div2u.c

uint32_t div2u(uint32_t x) {
    return x / 2;
}

div2u(unsigned int):
        mov     eax, edi
        shr     eax
        ret

$2^n$で除算するとは、もし$n \geq 0$ならば、$n$だけ右算術シフトさせるのに等しいことは有名です。符号付き整数を扱う場合にはシフトだけだと期待した値にはならないので、ちょっとした補正が必要です。

div2.c

uint32_t div2(int32_t x) {
    return x / 2;
}

div2(int):
        mov     eax, edi
        shr     eax, 31   ; 符号ビットを取り出し
        add     eax, edi  ; それを元の値に加え
        sar     eax       ; 右算術シフト
        ret

コンパイラがやっていることは、負の数の場合に1を足してからシフトすることです。これでなぜ上手くいくかですが、右シフトで捨てられる最下位ビットは、もとの値の正負に関わらず、重み+1を持っていました。つまりシフト量1での右算術シフトは、負の無限大への丸めをする割る2となります。

今回実現したい除算は0方向への丸めなので、正に関してはシフト結果と一致していますが負については1のズレが起こります。そこで負の場合に予め1を足して補正するというわけです。

定数3で割る問題

さて、実は定数2で割る場合だけでなく、定数3で割る場合についてもコンパイラは頑張ってくれます。これについても見ていきます。

div3.c

int32_t div3(int32_t x) {
    return x / 3;
}

div3(int):
        movsx   rax, edi
        sar     edi, 31
        imul    rax, rax, 1431655766
        shr     rax, 32
        sub     eax, edi
        ret

div3が何をやっているか解読します。

1. 元の値に定数1431655766を掛ける
2. 乗算結果の上32ビットを取り
3. 元の値の符号ビットを足す
4. 計算結果は下32ビットとする

かなり驚くべきことだと思うのですが、割る3が消えてしまいました。代わりに1431655766 という謎のマジックナンバーを使った実装が現れています。ここからは、コンパイラが採用しているこの技法を見ていきます。

割る3を乗算でやる技法について

改めて問題を確認します。話を単純にするため、入力は符号付き32ビット整数ということにしておきます²。

• 割る3を除算命令を使わずにやりたい
• 除算以外の整数演算とビットシフトは比較的高速にできる
  • ビットシフトなどを使えば$\div 2^n$は実現できる
• プログラム中に何らかの定数を使ってもよい

ここで思い浮かぶのは、除算は逆数の乗算に変換できることです。例えば$\div 3$は$\times 1/3$にできます。しかし、整数の範囲で$1/3$を精度よく表現することはできません。

発想は「もう3で割られている」大きな整数$M$を作ることです。なぜ$M$を大きくしたいかというと、その方が誤差が小さくなるからです（後述）。

$$
M = \frac{？}{3}
$$

割る対象の数に$M$をかけ、かけた値に右ビットシフトをして$？$の部分を除去すると、割る3部分だけが残るはずだ、というのが方針です。ビットシフトで消えるように、$？$には$2^n$になる大きな値を選びます。ただし$2^n$だけだと3の倍数にならないので、ほんの少し調整用の定数を足します。$M$は乗算に使うので、符号付き32ビット整数に収まらなければならない点に注意してください。

割る3における謎のマジックナンバー1431655766の正体は、このような発想から作られた整数$M$です。

$$
M = \frac{2^{32} + 2}{3} = 1431655766
$$

div3で入力に$M$をかけ、シフト量32で右シフトしているのは、これによって分母の３だけを残そうという方針なわけです。

いやいや、勝手に$+ 2$してるけど大丈夫なの？ と思われるかもしれませんが、実はこれは「誤差」として許容できます。すなわち、足してるけど最終的な結果に影響がないということです。その点も含めて、何を計算しているのかもう少し詳しく見てみます。

入力の整数を$n$と書くことにし、説明を簡単にするため$n \ge 0$のときのみを考えます。元の手順では$M$をかけていたので$n \times M$です。これにシフト量32で右シフトしていたので、結局は

$$
q = \left\lfloor \frac{n \times M}{2^{32}} \right\rfloor
$$

なる$q$を計算していたことになります。というのも、割る2の節で確認したように、$n$が非負のとき右シフトは除算と解釈できるのです。式を整理してみます。

$$
\begin{aligned}
q &= \left\lfloor \frac{n \times M}{2^{32}} \right\rfloor \
&= \left\lfloor \frac{n \times (2^{32} + 2)}{2^{32} \times 3} \right\rfloor \
&= \left\lfloor \frac{2^{32}n + 2n}{3 \times 2^{32}} \right\rfloor \
&= \left\lfloor \frac{n}{3} + \frac{2n}{3 \times 2^{32}} \right\rfloor
\end{aligned}
$$

さて、この$q$は、3で割る場合と等しいでしょうか？

一般に、床関数に関して次の補題が知られています。整数$n, d$（ただし$d \not = 0$）および実数$x$に対して、$0 \le x < \rvert \frac{1}{d} \lvert $ならば以下の式が成り立ちます。

$$
\left\lfloor \frac{n}{d} + x \right\rfloor = \left\lfloor \frac{n}{d} \right\rfloor
$$

式はちょっと複雑ですが、要するに「誤差」項$x$が十分小さいなら無視して（削除して）かまわないということです。

補題の$x$は、$q$でいう項 $\frac{2n}{3 \times 2^{32}}$ にあたります。$n$は32ビット符号付き整数と仮定しているので、$n < 2^{31}$であることを使って項の大きさを評価すると

$$
\frac{2n}{3 \times 2^{32}} < \frac{2 \times 2^{31}}{3 \times 2^{32}} = \frac{1}{3}
$$

したがって$q$は補題の条件にあてはまっているので

$$
q = \left\lfloor \frac{n}{3} + \frac{2n}{3 \times 2^{32}} \right\rfloor = \left\lfloor \frac{n}{3} \right\rfloor
$$

が成り立ちます。すなわち$q$は除算の定義そのものなので、最初の手続きによる計算式は、除算に一致することが示せました。
この結果は3の倍数を作るために$+2$してもよかったか、という疑問の答えにもなっています。式の評価を進めると実は問題にならないことが分かりました。

nが負の場合

前節の議論では$n < 0$の場合について確認していません。しかし同じような計算をすると同様に除算と一致します。
$n$が負の場合、div3では符号ビットを足しているので、$q$の式が$+1$された形に変化します。1は整数なので床関数の中にそのまま入れて問題ない点に注意しつつ、式を整理していきます。

$$
\begin{aligned}
q &= \left\lfloor \frac{n \times M}{2^{32}} \right\rfloor + 1 \
&= \left\lfloor \frac{2^{32}n + 2n}{3 \times 2^{32}} \right\rfloor + 1 \
&= \left\lfloor \frac{2^{32}n + 2n}{3 \times 2^{32}} + 1 \right\rfloor \
&= \left\lfloor \frac{2^{32}n + 2n + 3 \times 2^{32}}{3 \times 2^{32}} \right\rfloor \
\end{aligned}
$$

床関数と天井関数には次のような関係が知られています。整数$n, d( \not = 0)$について

$$
\left\lfloor \frac{n}{d} \right\rfloor = \left\lceil \frac{n - d + 1}{d} \right\rceil
$$

つまり、床関数の中身が分数の形をしているとき、分子に$-d + 1$を加えれば床関数を天井関数に置き換えられるという定理です。定理を用いて$q$の式変形を続けます。

$$
\begin{aligned}
q &= \left\lfloor \frac{2^{32}n + 2n + 3 \times 2^{32}}{3 \times 2^{32}} \right\rfloor \
&= \left\lceil \frac{2^{32}n + 2n + 3 \times 2^{32} + (-3 \times 2^{32} + 1)}{3 \times 2^{32}} \right\rceil \
&= \left\lceil \frac{2^{32}n + 2n + 1}{3 \times 2^{32}} \right\rceil \
&= \left\lceil \frac{n}{3} + \frac{2n + 1}{3 \times 2^{32}} \right\rceil \
\end{aligned}
$$

$n \ge 0$の場合と同様に$q$の「誤差」の大きさを評価します。$n$は32ビット符号付き整数でありかつ負なので、$-2^{31} \le n \le -1$を使って「誤差」の項を挟み込んでみます（実質的には$n$に$-2^{31}, -1$をそれぞれ代入して計算します）。

$$
-\frac{1}{3} < -\frac{1}{3} + \frac{1}{3 \times 2^{32}} \le \frac{2n + 1}{3 \times 2^{32}} \le - \frac{1}{3 \times 2^{32}} < 0
$$

すなわち「誤差」は0よりは小さいし、$n$が最小の時を見ても$-1/3$よりは大きいと分かりました。

先に使った補題には、実は天井関数バージョンのものがあります。整数$n, d$（ただし$d \not = 0$）および実数$x$に対して、$-\rvert \frac{1}{d} \lvert < x \le 0 $ならば以下の式が成り立ちます。

$$
\left\lceil \frac{n}{d} + x \right\rceil = \left\lceil \frac{n}{d} \right\rceil
$$

今回はこちらが適用できるので「誤差」部分を除去します。

$$
\begin{aligned}
q &= \left\lceil \frac{n}{3} + \frac{2n + 1}{3 \times 2^{32}} \right\rceil \
&= \left\lceil \frac{n}{3} \right\rceil
\end{aligned}
$$

この$q$もまた、除算の定義と一致していることが分かります。

これらの議論は一般化できるか？

割る2と割る3について触れました。それなら一般化して割るdについても最適化が考えられるのだろうか、となるのはもっともなことです。

この記事では時間の関係でまとめきれなかったのですが、割るdに対して除算命令を避けた命令生成を与えるアルゴリズムは存在します。単に定数Mをうまく選ぶだけではなく、dによって命令列を少しずつ変える必要もあります。

参考資料のリンク先にはより完全な解説が載っています。

まとめ

除算をコンパイラが最適化する場合の一部について、背後にあるアイデアから確認してみました。2と3の場合の一部だけでこれだけの分厚さになってしまったので、一般的なdについて解説するのはちょっと難しすぎました……。

とはいえ、除算を乗算に変換するという面白い技芸について紹介はできたかなと思います。

以上の議論は、実は『ハッカーのたのしみ―本物のプログラマはいかにして問題を解くか』に載っている内容です。ただ、筆者にとって解説があまりにも簡潔すぎて難しかったので、今回Qiita記事の形でまとめ直してみることにしました。

参考資料

• ハッカーのたのしみ―本物のプログラマはいかにして問題を解くか
• https://ridiculousfish.com/blog/posts/labor-of-division-episode-i.html
• https://uops.info/table.html
• Torbjörn Granlund and Peter L. Montgomery. 1994. Division by invariant integers using multiplication. SIGPLAN Not. 29, 6 (June 1994), 61–72. DOI:https://doi.org/10.1145/773473.178249

1. Cannonlake世代からidivは高速化されたため、ここまで悲惨ではなくなっていますhttps://www.anandtech.com/show/13699/intel-architecture-day-2018-core-future-hybrid-x86/2 ↩

2. 乗算が出現する技法なので、64ビット整数が入力だとオーバーフローについて考えざるを得ず、説明がもう少しややこしくなってしまいます。 ↩