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

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

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

シリーズの記事です。

  1. コンセプト
  2. 正答率
  3. 誤差の補正(メイン記事) ← この記事
  4. 測定

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

目次

概要

$a,b$ を非負の整数とします。

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


a/b=a×(1/b)=a×(x/b)/x

$b$ が定数であれば $x/b$ は事前に計算できます。$x$ が $2$ の累乗であればシフトで割り算が計算できるため、実行時の計算から割り算を排除できます。$x/b$ を切り捨てた値がゼロにならないようにするため条件を付けます。


\tag{1} x = 2^n\quad(x≧b)

この記事では $x/b$ が割り切れない場合に切り捨てによる誤差を補正して $a/b$ を求める方法を考えます。

※ $x/b$ が割り切れる場合、$x$ を割り切る $b$ は $1$ と $2$ の累乗だけです。$1$ での割り算は計算不要で、$2$ の累乗での割り算はシフトで処理できるため、容易に計算可能です。

整数計算において割り算は切り捨てられます。数式で切り捨ては床関数と呼ばれる特殊な括弧 $⌊\ ⌋$ で表します。


\text{【例】}\ \left⌊\frac 5 2\right⌋=2

割り切れる場合

$a/b$ が割り切れる場合について、商を $x$ 倍したときの切り捨てによる誤差 $d$ を考えます。$x/b$ が割り切れるならば $d=0$ となることに注意します。


\tag{2} \frac{ax}{b} = a\left⌊\frac{x}{b}\right⌋+d

\begin{aligned}
\tag{3}
d&=\frac{ax}{b}-a\left⌊\frac{x}{b}\right⌋ \\
&=a\left(\frac{x}{b}-\left⌊\frac{x}{b}\right⌋\right) \\
&=a\left(\frac{x\bmod b}{b}\right) \\
&=\frac{a}{b}(x\bmod b)
\end{aligned}

$d$ を求めるために $a/b$ が必要となりますが、$a/b$ は求めたい値そのものです。このような循環を回避するため $d$ を調べます。

範囲の制限

$d$ は商を $x$ 倍したときの誤差のため、最終的に商を求めるときに $d$ を $x$ で割ります。


\tag{4}
\begin{aligned}
\frac a b
&=\frac{ax}b\frac 1 x \\
&=\underbrace{\left(a\left⌊\frac{x}{b}\right⌋+d\right)}_{(2)}\frac 1 x \\
&=a\left⌊\frac{x}{b}\right⌋\frac 1 x+\frac d x
\end{aligned}

$d/x$ の取り得る範囲を調べます。


\begin{aligned}
0 ≦ x\bmod b&< b \\
0 ≦ \frac 1 b (x\bmod b) &< 1 \\
0 ≦ d = \frac a b(x\bmod b) &< a \\
∴\ 0 ≦ \frac d x &< \frac a x
\end{aligned}

条件を追加して範囲を制限します。


\tag{5}
x > a,\ 0 ≦ \frac d x < 1

誤差の補正

(4) を式変形します。


\begin{aligned}
a\left⌊\frac{x}{b}\right⌋\frac 1 x
&=\frac a b-\frac d x \\
\left⌊a\left⌊\frac{x}{b}\right⌋\frac 1 x\right⌋
&=\left⌊\frac a b-\frac d x\right⌋ \\
\end{aligned}

ここでは $a/b$ が割り切れることを想定しているため $a/b$ を床関数の外に出せます。


\begin{aligned}
\left⌊a\left⌊\frac{x}{b}\right⌋\frac 1 x\right⌋
&=\frac a b+\left⌊-\frac d x\right⌋ \\
\frac a b
&=\left⌊a\left⌊\frac{x}{b}\right⌋\frac 1 x\right⌋-\left⌊-\frac d x\right⌋
\end{aligned}

(5) を更に $d > 0$ に制限すれば $⌊-d/x⌋=-1$ となります。


\tag{6} \frac{a}{b} = \left⌊a\left⌊\frac{x}{b}\right⌋\frac{1}{x}\right⌋+1

(3) に (6) を代入すれば $d$ から $a/b$ の循環が消えます。


\tag{7}
\begin{aligned}
d&=\underbrace{\frac{a}{b}(x\bmod b)}_{(3)} \\
&=\underbrace{\left(\left⌊a\left⌊\frac{x}{b}\right⌋\frac{1}{x}\right⌋+1\right)}_{(6)}(x\bmod b)
\end{aligned}

(6) は $d > 0$ に制限して得られた式ですが、$x/b$ が割り切れるならば $x \bmod b=0$ より $d=0$ となるため、(7) は (5) の範囲すべてで成り立ちます。

計算式の導出

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


\tag{8}
\begin{aligned}
\frac{ax}{b}
&=\underbrace{a\left⌊\frac{x}{b}\right⌋+d}_{(2)} \\
&=a\left⌊\frac{x}{b}\right⌋+\underbrace{\left(\left⌊a\left⌊\frac{x}{b}\right⌋\frac{1}{x}\right⌋+1\right)(x\bmod b)}_{(7)}
\end{aligned}

(8) は複雑なので関数とします。


f(a) = a\left⌊\frac{x}{b}\right⌋+\left(\left⌊a\left⌊\frac{x}{b}\right⌋\frac{1}{x}\right⌋+1\right)(x\bmod b)

$x$ で割れば掛け算によって割り算をする式が得られます。


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

(9) は $a/b$ が割り切れることを前提としています。

割り切れない場合

ここまではすべて $a/b$ が割り切れる場合を考えてきました。割り切れない場合にも $f$ が適用できるかを調べます。

順序の保存

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


\begin{aligned}
&f(a+1)-f(a) \\
&=\left\{(a+1)\left⌊\frac{x}{b}\right⌋+\left(\left⌊(a+1)\left⌊\frac{x}{b}\right⌋\frac{1}{x}\right⌋+1\right)(x\bmod b)\right\} \\
&\quad-\left\{a\left⌊\frac{x}{b}\right⌋+\left(\left⌊a\left⌊\frac{x}{b}\right⌋\frac{1}{x}\right⌋+1\right)(x\bmod b)\right\} \\
&=\underbrace{\left⌊\frac{x}{b}\right⌋}_{≧1}+\underbrace{\left\{\left⌊(a+1)\left⌊\frac{x}{b}\right⌋\frac{1}{x}\right⌋-\left⌊a\left⌊\frac{x}{b}\right⌋\frac{1}{x}\right⌋\right\}(x\bmod b)}_{≧0} ≧ 1
\end{aligned}

よって $a$ が増加すれば $f(a)$ も増加します。つまり $f$ は順序を保存します。


\tag{10} f(a) < f(a+1)

(10) は $a/b$ が割り切れることを仮定していないため、割り切れなくても成立します。

はさみうち

$a/b$ が割り切れない場合、$b$ で割り切れる $a$ 以下の最大の数を $a'$ とします。


\tag{11} \frac{a'}b = \left⌊\frac{a}{b}\right⌋

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


a' < a < a'+b

(10) より $f$ は順序を保存するため


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

式を $x$ で割ります。


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

左側と右側は割り切れるため (9) と (11) より


\underbrace{\frac{f(a')}{x}=\frac{a'}{b}}_{(9)} < \frac{f(a)}{x} < \underbrace{\frac{f(a'+b)}{x}=\frac{a'+b}{b}}_{(9)}

\underbrace{\frac{a'}{b}=\left⌊\frac{a}{b}\right⌋}_{(11)} < \frac{f(a)}{x} < \underbrace{\left⌊\frac{a}{b}\right⌋}_{(11)}+1

$f(a)/x$ を切り捨てれば $a/b$ を切り捨てた値と一致します。


\tag{12} \left⌊\frac{a}{b}\right⌋ = \left⌊\frac{f(a)}{x}\right⌋

$a/b$ が割り切れれば $f(a)/x$ も割り切れますが、これは (12) に含まれます。

実装

$f$ をコード化します。数式とコードとの対応を示します。


\begin{aligned}
x&=\underbrace{2^n}_{\text{1 << n}}\quad(x>a,\ x≧b) \\
c&=a\left⌊\frac{x}{b}\right⌋ \\
f(a)&=c+\left(\underbrace{\left⌊\frac{c}{x}\right⌋}_{\text{c >> n}}+1\right)(\underbrace{x\bmod b}_{\text{x \% b}}) \\
\left⌊\frac a b\right⌋&=\underbrace{\left⌊\frac{f(a)}{x}\right⌋}_{\text{f >> n}}
\end{aligned}

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%)

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