引き続き整数の割り算を掛け算に変換する方法について考えます。
シリーズの記事です。
【注】第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%)
問題なく計算できています。
