ここまで見てきた計算方法で、どれくらい高速化するか実験します。
シリーズの記事です。
- コンセプト
- 正答率
- 誤差の補正（メイン記事）
- 測定 ← この記事
【注】第2回までの計算は不完全なため誤差が出ます。第3回で誤差を補正します。
目次
測定結果
サンプルでは定数として10と97を使用します。Athlon 64 X2 4600+ 2.4GHzでの測定結果を掲載します。
| div10 | 除算命令 div | 80.66s |
| div10mul | 乗算命令 mul | 16.14s |
| div97 | 除算命令 div | 82.46s |
| div97mul | 乗算命令 mul を2回（後者は誤差補正） | 21.67s |
| div97mul2 | 乗算命令 mul と算術命令による誤差補正 | 17.16s |
除算命令が乗算命令の4～5倍ほど遅いことが分かります。誤差補正を工夫すると、少し速くなります。誤差補正は小さい数での掛け算を加減算やシフトに置き換えています。複雑に見えますが、いくつかパターンがありそうです。
計算式
除算を乗算に変換して誤差を補正するのは、第3回の記事で導出した次の式によります。a,b は非負の整数です。bが定数であるという前提で x / b と x % b を事前に計算して埋め込むことを想定しています。
\begin{aligned}
x&=2^n\quad(x>a,\ x≧b) \\
c&=a\left⌊\frac{x}{b}\right⌋ \\
f&=c+\left(\left⌊\frac{c}{x}\right⌋+1\right)(x\bmod b) \\
\left⌊\frac{a}{b}\right⌋&=\left⌊\frac f x\right⌋
\end{aligned}
いくつか実例で確認します。
b=10
\begin{aligned}
x&=2^{32}\quad(x>a,\ x≧10) \\
c&=a\underbrace{\left⌊\frac{2^{32}}{10}\right⌋}_{429496729} \\
f&=c+\left(\left⌊\frac{c}{2^{32}}\right⌋+1\right)\underbrace{(2^{32}\bmod 10)}_{6} \\
\left⌊\frac{a}{10}\right⌋&=\left⌊\frac f{2^{32}}\right⌋
\end{aligned}
b=97
\begin{aligned}
x&=2^{32}\quad(x>a,\ x≧97) \\
c&=a\underbrace{\left⌊\frac{2^{32}}{97}\right⌋}_{44278013} \\
f&=c+\left(\left⌊\frac{c}{2^{32}}\right⌋+1\right)\underbrace{(2^{32}\bmod 97)}_{35} \\
\left⌊\frac{a}{97}\right⌋&=\left⌊\frac f{2^{32}}\right⌋
\end{aligned}
実装
i386の乗除算命令の仕様を確認します。
mul src # eax * src → edx(上位):eax(下位) div src # edx(上位):eax(下位) / src → eax(商) ... edx(剰余)
x = 1 << 32 のとき、mul命令で帰ってくるedx（上位32bit）が商となります。bが32bitのため x > a, x > b が保証されます。
b = 10 と b = 97 での計算をアセンブラで実装します。
.intel_syntax .globl _div10 _div10: xor edx, edx # edx = 0 mov eax, [esp + 4] # eax = a mov ecx, 10 # ecx = 10 div ecx # eax ... edx = edx:eax / ecx ret # return eax .globl _div10mul _div10mul: mov eax, 0xcccccccd # eax = (2^36 / 10 + 1) / 2 mul dword ptr [esp + 4] # edx:eax = eax * a mov eax, edx shr eax, 3 ret # return (edx / 8) .globl _div97 _div97: xor edx, edx # edx = 0 mov eax, [esp + 4] # eax = a mov ecx, 97 # ecx = 97 div ecx # eax ... edx = edx:eax / ecx ret # return eax .globl _div97mul _div97mul: push ebx mov eax, 0x2a3a0fd # eax = 2^32 / 97 mul dword ptr [esp + 8] # edx:eax = eax * a mov ebx, eax mov ecx, edx # c = edx:eax mov eax, 35 # eax = 2^32 % 97 inc edx mul edx # edx:eax = eax * (edx + 1) add eax, ebx adc edx, ecx # edx:eax += c mov eax, edx pop ebx ret # return edx .globl _div97mul2 _div97mul2: mov eax, 0x2a3a0fd # eax = 2^32 / 97 mul dword ptr [esp + 4] # edx:eax = eax * a inc edx # ++edx lea ecx, [edx + edx * 2] lea ecx, [ecx + ecx * 2] sal ecx, 2 sub ecx, edx add eax, ecx adc edx, 0 # edx:eax += edx * (2^32 % 97) lea eax, [edx - 1] ret # return edx - 1
32bit総当りの計算速度を測定します。
#include <stdio.h> #include <stdint.h> #include <time.h> extern uint32_t div10(uint32_t); extern uint32_t div10mul(uint32_t); extern uint32_t div97(uint32_t); extern uint32_t div97mul(uint32_t); extern uint32_t div97mul2(uint32_t); int main(void) { clock_t t1, t2; uint32_t i; i = 0; t1 = clock(); do div10(i++); while (i != 0); t2 = clock(); printf("div10: %.2fs\n", ((double)t2 - t1) / CLOCKS_PER_SEC); i = 0; t1 = clock(); do div10mul(i++); while (i != 0); t2 = clock(); printf("div10mul: %.2fs\n", ((double)t2 - t1) / CLOCKS_PER_SEC); i = 0; t1 = clock(); do div97(i++); while (i != 0); t2 = clock(); printf("div97: %.2fs\n", ((double)t2 - t1) / CLOCKS_PER_SEC); i = 0; t1 = clock(); do div97mul(i++); while (i != 0); t2 = clock(); printf("div97mul: %.2fs\n", ((double)t2 - t1) / CLOCKS_PER_SEC); i = 0; t1 = clock(); do div97mul2(i++); while (i != 0); t2 = clock(); printf("div97mul2: %.2fs\n", ((double)t2 - t1) / CLOCKS_PER_SEC); return 0; }
これで今回のシリーズは終了です。
