コンパイラの出力を見ていると、除数が定数の場合に整数の割り算を掛け算に変換していることに気付きます。これがどんな原理で行われているのか考えてみました。
シリーズの記事です。
- コンセプト ← この記事
- 正答率
- 誤差の補正(メイン記事)
- 測定
【注】第2回までの計算は不完全なため誤差が出ます。第3回で誤差を補正します。
目次
概要
割り算は除数を分母とする掛け算に変換できます。
- a / b → a * (1 / b)
整数だけで計算する場合、商は切り捨てられるため、1 / b は b > 1 のときゼロとなってしまいます。これでは意味がありません。これを回避するため、乗数の分子を巨大数として後で割る方法が考えられます。巨大数をxとします。
- a / b → a * (x / b) / x
x / b をコンパイル時に計算して定数として埋め込むことで、実行時の割り算を排除します。しかし後処理の / x で割り算が発生すると意味がありません。割り算を排除するためシフトで処理するなら、xは2の冪乗である必要があります。
x86の乗算命令(mul)は、被乗数と乗数のサイズは同じですが、積は倍のサイズで提供されます。16bitを例に取ると以下の通りです。
- 16bit(AX) * 16bit = 32bit: 上位16bit(DX), 下位16bit(AX)
つまり16bitの割り算は、計算過程の許容ビット数が32bitとなります。単純に上位16bitを商として取り出すならxは16bit最大値+1で0x10000となります。
- a / b ≒ (a * (0x10000 / b)) >> 16
実装
C言語で16bit符号なし整数型の割り算を実装して、適当にテストします。
- div.c
#include <stdio.h> #include <stdint.h> int main(void) { uint16_t a, b; for (b = 2; b < 100; b++) { uint32_t f = 0x10000 / b; for (a = 0; a < 100; a++) { uint16_t c1 = a / b; uint16_t c2 = (a * f) >> 16; if (c1 != c2) printf("[NG] %d / %d = %d != %d; 0x%x\n", a, b, c1, c2, f); } } return 0; }
実行結果
[NG] 3 / 3 = 1 != 0; 0x5555 [NG] 6 / 3 = 2 != 1; 0x5555 [NG] 9 / 3 = 3 != 2; 0x5555 [NG] 12 / 3 = 4 != 3; 0x5555 [NG] 15 / 3 = 5 != 4; 0x5555 [NG] 18 / 3 = 6 != 5; 0x5555 (以下略)
それっぽい値は出ていますが、割り切れるはずの計算がうまくできないようです。0x10000 / b を切り捨てているのが原因です。コンパイラの出力を見ると単純な掛け算ではなく、複雑な後処理をしています。誤差についてきちんと検討する必要がありそうです。