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

引き続き整数の割り算を掛け算に変換する方法について考えます。

シリーズの記事です。

【注】第2回までの計算は不完全なため誤差が出ます。第3回で誤差を補正します。

目次

概要

カウント

切り上げ (1)

精度を上げる

切り上げ (2)

チェックの高速化

概要

単純化のため、計算対象を16bit符号なし整数型に限定します。

除数が最大値の半分より大きい（> 32767）場合、商は0か1で、単純な大小比較で求まります。これは掛け算よりも高速です。そのため対象から除外します。

カウント

第1回は以下の式で計算しました。

a / b ≒ (a * (0x10000 / b)) >> 16

除数 0～32767 のうち、この計算方法で誤差がない除数を数えます。

#include <stdio.h>
#include <stdint.h>

int main(void) {
    uint32_t count = 0;
    uint16_t b, max = UINT16_MAX / 2;
    for (b = 1; b <= max; b++) {
        uint32_t a, f = 0x10000 / b;
        uint16_t ok = 1;
        for (a = 0; a <= UINT16_MAX; a++) {
            uint16_t c1 = a / b;
            uint16_t c2 = (a * f) >> 16;
            if (c1 != c2) {
                ok = 0;
                break;
            }
        }
        count += ok;
    }
    printf("%u / %u (%u%%)\n", count, max, count * 100 / max);
    return 0;
}

実行結果

15 / 32767 (0%)

これは酷いです。正常に通った被除数は1と2の累乗だけです。1は計算不要、2の累乗はシフトで計算可能のため、まったく意味がありません。

切り上げ (1)

乗数を求めるとき切り上げるよう変更します。

        uint32_t a, f = (0x10000 + b - 1) / b;

実行結果

15 / 32767 (0%)

まったく精度が向上しません。

精度を上げる

今は積の上位16bitを取っているだけです。除数が2以上の自然数の場合、商は被除数よりも小さくなります。減少した桁数は上位16bitの中で使われることはないため、乗数を大きくすれば有効活用できそうです。

2の累乗を基準にすると2進数での最大桁数が判別できます。たとえば8での除算は被除数 >> 3 で、16での除算は被除数 >> 4 のため、その間の9～15での除算は8を基準にすると最小でも3桁は減少します。これは被除数を2進数で表記したときの桁数そのものです。

これを実装します。

#include <stdio.h>
#include <stdint.h>

int count2(uint16_t x) {
    int ret = 0;
    for (; x; x >>= 1, ret++);
    return ret;
}

int main(void) {
    uint32_t count = 0;
    uint16_t b, max = UINT16_MAX / 2;
    for (b = 1; b <= max; b++) {
        int shift = 15 + count2(b);
        uint32_t a, f = (1u << shift) / b;
        uint16_t ok = 1;
        for (a = 0; a <= UINT16_MAX; a++) {
            uint16_t c1 = a / b;
            uint16_t c2 = (a * f) >> shift;
            if (c1 != c2) {
                ok = 0;
                break;
            }
        }
        count += ok;
    }
    printf("%u / %u (%u%%)\n", count, max, count * 100 / max);
    return 0;
}

実行結果

15 / 32767 (0%)

まったく精度が向上しません。

切り上げ (2)

再度、切り上げるよう変更します。

        uint32_t a, f = ((1 << shift) + b - 1) / b;

実行結果

24969 / 32767 (76%)

劇的に精度が向上しました。

チェックの高速化

計算結果が正しいかどうか総当りで検算しましたが、かなり時間が掛かります。

掛け算を行っているため、誤差は被除数に比例して拡大します。対象となる型での最大の割り切れる数と、その直前の数だけを検算すれば確認できます。

#include <stdio.h>
#include <stdint.h>

int count2(uint16_t x) {
    int ret = 0;
    for (; x; x >>= 1, ret++);
    return ret;
}

int main(void) {
    uint32_t count = 0;
    uint16_t b, max = UINT16_MAX / 2;
    for (b = 1; b <= max; b++) {
        int shift = 15 + count2(b);
        uint32_t f = ((1 << shift) + b - 1) / b;
        uint16_t n2 = UINT16_MAX / b * b;
        uint16_t n1 = n2 - 1;
        uint16_t c1 = n1 / b;
        uint16_t c2 = n2 / b;
        uint16_t d1 = (n1 * f) >> shift;
        uint16_t d2 = (n2 * f) >> shift;
        if (c1 == d1 && c2 == d2) count++;
    }
    printf("%u / %u (%u%%)\n", count, max, count * 100 / max);
    return 0;
}