読者です 読者をやめる 読者になる 読者になる

【お知らせ】プログラミング記事の投稿はQiitaに移行しました。

整数の割り算を掛け算に変換 (4)

C言語 x86

前回まで見てきた計算方法で、どれくらい高速化するか実験します。

サンプルでは定数として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倍ほど遅いことが分かります。誤差補正を工夫すると、少し速くなります。誤差補正は小さい数での掛け算を加減算やシフトに置き換えています。複雑に見えますが、いくつかパターンがありそうです。

除算を乗算に変換して誤差を補正するのは次の式によります。bが定数であるという前提で x / b と x % b を事前に計算してハードコードします。

\left\lfloor\frac{a}{b}\right\rfloor = \left\lfloor\frac{a\left\lfloor\frac{x}{b}\right\rfloor+\left(a\left\lfloor\frac{x}{b}\right\rfloor\frac{1}{x}+1\right)(x\bmod b)}{x}\right\rfloor

i386の乗除算命令の仕様

  • div src = edx(上位):eax(下位) / src → eax(商) … edx(剰余)
  • mul src = eax * src → edx(上位):eax(下位)

x = 1 << 32 とすれば、mul命令で帰ってくるedx(上位32bit)が商となります。

.intel_syntax

.globl _div10
_div10:
	xor edx, edx
	mov eax, [esp + 4]
	mov ecx, 10
	div ecx
	ret

.globl _div10mul
_div10mul:
	mov eax, 0xcccccccd # ((1 << 36) / 10 + 1) >> 1
	mul dword ptr [esp + 4]
	mov eax, edx
	shr eax, 3
	ret

.globl _div97
_div97:
	xor edx, edx
	mov eax, [esp + 4]
	mov ecx, 97
	div ecx
	ret

.globl _div97mul
_div97mul:
	push ebx
	mov eax, 0x2a3a0fd # (1 << 32) / 97
	mul dword ptr [esp + 8]
	mov ebx, eax
	mov ecx, edx
	mov eax, 35 # (1 << 32) % 97
	inc edx
	mul edx
	add eax, ebx
	adc edx, ecx
	mov eax, edx
	pop ebx
	ret

.globl _div97mul2
_div97mul2:
	mov eax, 0x2a3a0fd
	mul dword ptr [esp + 4]
	inc edx
	lea ecx, [edx + edx * 2]
	lea ecx, [ecx + ecx * 2]
	sal ecx, 2
	sub ecx, edx
	add eax, ecx
	adc edx, 0
	lea eax, [edx - 1]
	ret

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;
}

これで今回のシリーズは終了です。