コンパイラによる整数除算最適化の証明

一般的に、CPUにおける整数の除算（割り算）は遅いです。

そのため、コンパイル時定数で割る処理を書くと、コンパイラが「乗算とビットシフト」を用いた高速な処理へ変換します。

この記事では、その具体的な変換方法と正当性を示します。

はじめに

この記事では上記のような最適化の存在を知っていることを前提とします。 知らない場合は、2や3で割る場合についての以下の記事を先に読むことをおすすめします。

2で割ることと3で割ること - Qiita

概要

3で割るとき

たとえば、以下の定数 $3$ で割って返す関数は、

1
2
3
4
5
6
7

#include<stdint.h>

const uint32_t d = 3;

uint32_t div_by_d(uint32_t x) {
    return x / d;
}

x86_64 gcc 12.2 では次のとおり最適化されます。 こちらのサイトで実際に $d$ を変えて試すことができます。

div_by_d(unsigned int):
        push    rbp
        mov     rbp, rsp
        mov     DWORD PTR [rbp-4], edi
        mov     eax, DWORD PTR [rbp-4]
        mov     edx, eax
        mov     eax, 2863311531
        imul    rax, rdx
        shr     rax, 32
        shr     eax
        pop     rbp
        ret

上記のアセンブリでは、C++で記述すると次のような処理に変換されています。 64ビットで $x$ を扱い、 $2863311531$ をかけて $33$ ビット右へシフトしています。

1
2
3
4
5
6

uint32_t div_by_d(uint32_t x) {
    uint64_t t = x;
    t *= (uint64_t) 2863311531;
    t >>= (uint64_t) 33;
    return (uint32_t) t;
}

そして実は、任意の $d (\neq 0)$ に対して上記のような形に最適化できます。 すなわち、任意の $d (\neq 0)$ に対して、次のような形に最適化できる $m, l$ が存在します。

1
2
3
4
5
6

uint32_t div_by_d(uint32_t x) {
    t = x;
    t *= m;
    t >>= 32 + l;
    return (uint32_t) t;
}

記事の内容

以下ではまず、実際に $d$ に対して $m, l$ を得るアルゴリズムと、その正当性を示します。

また、上記のコードのとおりにプログラムを書くと、 $m$ の値によってはオーバーフローすることがあります。 これを回避する方法も示します。

これらは元論文^[1]を参考に、簡潔な形で整理したものになります。

$m, l$ を得るアルゴリズム

次のとおりです。

$$ l = \lceil \log_2 (d) \rceil, \ m = \left\lceil \frac{2^{32 + l}}{d} \right\rceil. $$

どちらもシンプルなアルゴリズムで十分高速（$O(\log d)$ 時間）に求まります。

正当性

上記の $m, l$ で正しく計算できることを示します。

1
2
3
4
5
6

uint32_t div_by_d(uint32_t x) {
    t = x;
    t *= m;
    t >>= 32 + l;
    return t;
}

上記の関数の返り値は、数式で表すと $\left\lfloor mx / 2^{32+l} \right\rfloor$ です。 これが元の x / d に一致していること、すなわち次の式を示せば十分です。

$$ \begin{align} \left\lfloor \frac{mx}{2^{32+l}} \right\rfloor = \left\lfloor \frac{x}{d} \right\rfloor. \end{align} $$

まず、$m \geq 2^{32 + l} / d$ ですから、次のとおり評価できます。

$$ \begin{align*} \left\lfloor \frac{mx}{2^{32+l}} \right\rfloor &\geq \left\lfloor \frac{x}{2^{32+l}} \cdot \frac{2^{32+l}}{d} \right\rfloor \\ &= \left\lfloor \frac{x}{d} \right\rfloor. \end{align*} $$

一方で、 $x = q d + r \ (0 \leq r < d)$ とすると、次のとおり評価できます。

$$ \begin{align*} \frac{mx}{2^{32+l}} - \left\lfloor \frac{x}{d} \right\rfloor &= \left\lceil \frac{2^{32 + l}}{d} \right\rceil \frac{x}{2^{32+l}} - q \\ &\leq \left( \frac{2^{32 + l}}{d} + 1\right) \frac{x}{2^{32+l}} - q \\ &= \frac{x}{d} + \frac{x}{2^{32+l}} - q \\ &= \frac{x}{d} + \frac{x}{2^{32 + \lceil \log_2 (d) \rceil}} - \frac{x-r}{d} \\ &= \frac{r}{d} + \frac{x}{2^{32 + \lceil \log_2 (d) \rceil}} \\ &\leq \frac{d - 1}{d} + \frac{x}{2^{32} \cdot d} \\ &\leq \frac{1}{d} \left( d - 1 + \frac{2^{32} - 1}{2^{32}} \right) \\ &< 1. \end{align*} $$

以上より、(1)式が示せました。

オーバーフローを避ける方法

以上より、以下の計算を期待通り行うことができれば、 $d$ で割った値が得られます。

1
2
3
4
5
6

uint32_t div_by_d(uint32_t x) {
    t = x;
    t *= m;
    t >>= 32 + l;
    return t;
}

一方で、 $m$ は次のとおり $2^{32}$ を超える可能性があります。 そのような場合、 $x$ が大きいと $64$ ビットにおさまらずオーバーフローしてしまいます。

$$ \begin{align*} m &= \left\lceil \frac{2^{32 + l}}{d} \right\rceil \\ &= \left\lceil \frac{2^{32 + \lceil \log_2 (d) \rceil}}{d} \right\rceil \\ &\leq \left\lceil \frac{2^{33} \cdot d}{d} \right\rceil \\ &= 2^{33}. \end{align*} $$

したがって、$m$ が $2^{32}$ 以上となる場合、以下のように分割して計算します。

$$ \begin{align*} mx = (m - 2^{32}) x + 2^{32} x. \end{align*} $$

$(m - 2^{32}) \leq 2^{33} - 2^{32} = 2^{32}$ですから、$(m - 2^{32}) x, 2^{32} x$ はそれぞれオーバーフローせずに計算できます。

そして、以下の理由により、$(m - 2^{32}) x, 2^{32} x$ をそれぞれ $32$ ビット右にシフトしてから計算しても、値が変わらないことがわかります。

• $mx$ の上位 $32 - l$ ビットにのみ興味がある
• $2^{32} x$ の下位 $32$ ビットが $0$ である

したがって、以下のように計算することで、オーバーフローを避けることができます。

1
2
3
4
5
6
7
8

uint32_t div_by_d(uint32_t x) {
    uint64_t t = x;
    t *= (uint64_t) m - ((uint64_t) 1 << (uint64_t) 32);
    t >>= (uint64_t) 32;
    t += x;
    t >>= l;
    return t;
}

補足

最適化の気持ち

実数の除算で $x / d$ という計算を考えます。 これは、$(1 / d)$ をコンパイル時に計算しておくことで、 $x \cdot (1 / d)$ と乗算に変換できます。

実態としては、これを $m, l$ と追加の$32$ビットを用いて近似計算しています。 具体的には、

$$ \frac{m}{2^{32+l}} = \frac{1}{d} + \varepsilon $$

となるように $(1 / d)$ を固定小数点数として近似します。 このとき、 $x$ に $m$ をかけて$32 + l$ ビット右シフトすることで、

$$ \left\lfloor \frac{mx}{2^{32+l}} \right\rfloor = \left\lfloor \frac{x}{d} + \varepsilon x \right\rfloor $$

が得られます。 上手に $m, l$ を選んであげることで、誤差の $\varepsilon x$ が $x < 2^{32}$ において十分小さくなり、$\left\lfloor x / d \right\rfloor$ が得られます。 証明では、これを示していたのでした。

$N$ ビット整数での除算

今回の証明において、$32$ ビットを $N$ ビットとしても問題ありません。 したがって、たとえば $64$ ビットの整数同士で演算したい場合は、 $128$ ビット用いることで同様に計算できます。

$d$ が特殊な場合

先述の記事にも書かれていますが、$d = 2^k$ の場合は右に$k$シフトさせるとさらに高速になります。 また、それ以外の場合でも、より高速にする手法が色々あるようです。

符号付きの場合

符号付きの場合も同様な手法が適用できるようですが、符号なしを理解して満足してしまったので勉強していません。 気になる方は元論文¹を参照ください。