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

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

整数の割り算を掛け算に変換 (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;
}

かなり高速化しました。チェック結果が総当りと一致することを確認しました。

残念ながらまだ残り半分の除数が正常に計算できません。引き続き検討します。