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

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

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

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

割り算を掛け算に変換するのは逆数を利用します。整数だけの計算では2以上の逆数がゼロになってしまうため、事前にx倍して後で割ります。

  • a / b = a * (1 / b) = a * (x / b) / x

x / b を事前に計算して、xが2の累乗ならシフトで計算できるため、実行時の計算から割り算を排除できます。ただし事前に計算した x / b を切り捨てて整数値にすることで誤差が生じます。

x = 2^n \ \ \cdots\ (1)

割り切れる場合

a / b が割り切れる場合、x倍した状態での誤差dを考えます。

\frac{ax}{b} = a\left\lfloor\frac{x}{b}\right\rfloor+d \ \ \cdots\ (2)
d = \frac{ax}{b}-a\left\lfloor\frac{x}{b}\right\rfloor = a\left(\frac{x}{b}-\left\lfloor\frac{x}{b}\right\rfloor\right) \ \ \cdots\ (3)

帯分数

x / b を帯分数で表記します。

\frac{x}{b} = \left\lfloor\frac{x}{b}\right\rfloor+\frac{x\bmod b}{b} \ \ \cdots\ \ (4)

(3)に(4)を代入します。

d = a\left(\left\lfloor\frac{x}{b}\right\rfloor+\frac{x\bmod b}{b}-\left\lfloor\frac{x}{b}\right\rfloor\right) = a\cdot\frac{x\bmod b}{b} = \frac{a}{b}(x\bmod b) \ \ \cdots\ (5)

x % b がゼロ、つまり x / b が割り切れる場合、誤差が発生しないことが分かります。xは2の累乗のため、この条件を満たすbは1と2の累乗だけです。1での除算は計算不要で、2の累乗での除算はシフトで処理できるため、これらは除外して考えます。

矛盾

dを求めるために a / b が必要となっていますが、a / b は求める解そのものです。このような循環は矛盾しています。どうしたものでしょうか。

最大の誤差

aとbは同じサイズの整数型です。x / b を切り捨てた値がゼロにならないようにするため、b < x という条件が必要です。誤差の補正を前提にすれば、xの値をbに応じて決めずに一律にしても問題ありません。xはaやbの型の 最大値 + 1 とすれば、aの最大値は x - 1 となります。

これにより最大の誤差を求めます。

d_{max} = \frac{x-1}{b}(x\bmod b)

dはx倍されているため、商の誤差は1未満です。

\frac{d_{max}}{x} = \frac{x-1}{bx}(x\bmod b) \approx \frac{x\bmod b}{b} < 1

xは十分に大きいため近似しました。

\frac{x-1}{x} \approx 1

bの剰余の最大値は b - 1 です。

x\bmod b < b
\frac{x\bmod b}{b} < 1

結論

商は切り捨てられるため、誤差が1未満であっても、ゼロでなければ表面化します。つまり誤差が発生すれば商は1減少します。

x / b が割り切れなければ誤差が発生するため、a / b が割り切れる場合、誤差は1です。

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

(5)に(6)を代入します。

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

(2)に(7)を代入します。

\frac{ax}{b} = 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) \ \ \cdots\ (8)

(8)は a / b が割り切れる場合は常に成り立ちます。誤差は完全に補正されます。

関数

(8)の右辺が長いため関数にします。

f(a) = 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)

定義より

\frac{ax}{b} = f(a)

両辺をxで割ります。

\frac{a}{b} = \frac{f(a)}{x} \ \ \cdots\ (9)

ただし a / b が割り切れるときに限ります。

増加

aが増えたときにf(a)の変化を調べます。aは整数のため1で比較します。

f(a+1)-f(a) = \left\lfloor\frac{x}{b}\right\rfloor\left(\frac{x\bmod b}{x}+1\right) > 0
f(a) < f(a+1) \ \ \cdots\ (10)

aが増加すればf(a)も増加します。

割り切れない場合

a / b が割り切れない場合、a以下で最大の割り切れる数a'は、次のように求められます。

a' = \left\lfloor\frac{a}{b}\right\rfloor b

a'の次の割り切れる数は a' + b のため、次の関係が成り立ちます。

a' < a < a'+b

(10)より

f(a') < f(a) < f(a'+b)

式をxで割ります。

\frac{f(a')}{x} < \frac{f(a)}{x} < \frac{f(a'+b)}{x}

(9)より

\frac{a'}{b} < \frac{f(a)}{x} < \frac{a'+b}{b}
\frac{a'}{b} < \frac{f(a)}{x} < \frac{a'}{b}+1

f(a) / x を切り捨てると直前の整数と等しくなります。

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

よって、a / b が割り切れない場合でも f(a) / x は商と一致します。

実装

(8)をコード化します。nの定義は(1)です。

uint16_t div(uint16_t a, uint16_t b) {
    int n = 16;
    uint32_t x = 1 << n;
    uint32_t c = a * (x / b);
    uint32_t f = c + ((c >> n) + 1) * (x % b);
    return f >> n;
}

数式を弄り回した割にコードはあっさりしています。

掛け算を2回しています。bは定数で x / b と x % b は事前に決まります。除数が小さければ後者はあまり大きな数にはならないため、シフトや足し算で表現するなど最適化を工夫する余地があります。

テスト

総当りでテストします。

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

#define N 16

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

実行結果

32767 / 32767 (100%)

問題なく計算できています。