コンパイル結果を見ていた時にソースでは除算を使っているのに出力されたコードには使われていなかったことがあり、どういう仕組みかを調べた。
ググったところ、Stack Overflowがヒットした: c - Why does GCC use multiplication by a strange number in implementing integer division? - Stack Overflow。 要するに除算は重いから、定数での除算は最適化で乗算やシフトなどの軽い演算に置き換えられるとか。 そこで上記の回答に挙げられていた論文を読んでみる。
“Division by Invariant Integers using Multiplication”
論文のpdf:“Division by Invariant Integers using Multiplication”
Nビット同士の除算で簡単な順に、(i)符号なしの場合、(ii)符号ありで商を0に向かって丸める、(iii)符号ありで商を−∞に丸める、手順が書かれている。 まずは(i)符号なしの場合を読む。
符号なしの場合、導出
論文の4章に符号なしNビット同士の除算について書かれている。
割る数を \(d\)(不変)、割られる数を \(n\)(可変)とした時の商を \(q = \lfloor \frac{n}{d} \rfloor\) として、\(\lfloor \frac{n}{d} \rfloor = \lfloor \frac{m * n}{2^{N + l}} \rfloor\) (4.1式)となる \(m\) を求めるための不等式の導出がわかりにくい…。
まず \(l\) がなにかがこの時点では書かれてないように見える。
後の図4.1に \(l = \lceil \log_2 d \rceil\) と書かれているので、割る数 \(d\) の log2(=2の何乗以下か)を示しているぽい。
\(n = d\) とすると \(2^{N + l} \le m * d\) となる、というのは(4.1)式に代入すれば左辺が1で、それに等しくなるためには \(m * d\) が \(2^{N+l}\) 以上である必要がある、ということから不等式が導かれる。
次の上限を抑える不等式も間をすっ飛ばしているので導出がわかりにくい…。 その上 \(n = q * d - 1\) とするというのは、もともと \(q = \lfloor \frac{n}{d} \rfloor\) なので成り立つわけがなくおかしいので、 変数名を \(q\) から例えば \(i\) に変えて \(n = i * d - 1\) などとする必要があるように思う。
そうして代入すると左辺は \(\lfloor \frac{i * d - 1}{d} \rfloor = \lfloor i - \frac{1}{d} \rfloor = i - 1\)、 右辺 \(\lfloor \frac{m * (i * d - 1)}{2^{N + l}} \rfloor\) が \(i-1\) となるためには分子が分母*iより小さくなる必要がある、という導出になる。
\(d\) を掛けて〜というのは、先ほどのに乗じた \(2^{N+l} * (i * d) = 2^{N+l} * (i * d - 1 + 1) = 2^{N+l} * (i * d - 1) + 2^{N+l}\) となって、移項して \(i * d - 1\) で括れば導出できる。
そんなこんなで定理4.2で不等式(4.3)とその証明が示されるが、あまり理解できず…。
符号なしの場合、アルゴリズム
図4.1 に必要な係数と演算方法が、またその改良を施した版が図4.2に示される。
割る値 \(d\) とそれから算出した \(m\) によって場合わけされる
(MULUH=Nビット同士の乗算結果2Nビットの上位Nビット(符号なし)、SRL=論理右シフト):
- \(d\) が2のべき乗の場合
- ビットシフトのみ:
q = SRL(n, l);
- ビットシフトのみ:
- \(m \ge 2^N\) の場合
t1 = MULUH(m - 2^N, n); q = SRL(t1 + SRL(n - t1, 1), sh_post - 1);
- その他
q = SRL(MULUH(m, SRL(n, sh_pre)), sh_post);
m, sh_post, l は図6.2のCHOOSE_MULTIPLIERで算出する。
残り
5章符号あり除算(0に向かって丸め)ほか書かれている(が理解はできず…)
“Improved division by invariant integers”
https://gmplib.org/~tege/division-paper.pdf
- 除算を1つの
umulと 1つのumulloで計算する、30%ほど高速化、とか
“Reciprocal Multiplication, a tutorial”
同じくStack Overflowに回答されていたリンクを見てみる: Jones on reciprocal multiplication
- マイクロプロセッサやDSPチップからは整数除算命令が削除されてるとか?
- 乗算命令もなくした方がいいとか?
- 固定小数点での逆数乗算の解説
- 乗算をシフトに置き換える話
- RISCマシンには加算&シフト命令がある
- 符号あり被除数の場合の解説
実際の数値の例が多くてわかりやすい。
ちょっと理解が中途半端になってしまったが、力尽きたのでこの辺で終わり…。