除算
$n+1$ ワード $\div n$ ワード
割り算での基本は筆算での $n+1$ ワード $\div n$ ワードなので、その演算を考えよう。 つまり $x\div y$ の演算を行いたいとき、 \[ R^{n-1} \leq y \lt R^{n} \leq x \lt R^{n+1} \]
であるとする。また、商が $R$ 以上にならないこと $(yR \gt x)$ も確認が容易なのでその前提も置いておく。 このとき、商 $q$ $(0\le q \lt R)$ を立てたいのだが、一発で求めるのは無理があるので 上位桁から推測し補正することになる。具体的には \[ \hat{q} = \min \left( \left\lfloor \frac{x_{n}R+x_{n-1}}{y_{n-1}} \right\rfloor, R-1 \right) \]
から正しい商 $q$ へ補正する。 このとき、$y_{n-1} \geq \dfrac{R}{2}$ であれば $q \le \hat{q} \le q+2$ となることを保証できるので、その条件をみたすように $x,y$ の両方に $2^k$ をかけて計算する。
$q \leq \hat q$ の確認
前提として $q \leq R-1$ なので $\hat q = R-1$ の時は問題ない。 なので $\hat q = \left\lfloor \dfrac{x_{n}R+x_{n-1}}{y_{n-1}} \right\rfloor$ を考えると \[ \hat qy_{n-1} \ge x_{n}R+x_{n-1} \ge x_{n}R+x_{n-1} - y_{n-1} + 1 \] \[ \begin{eqnarray} \therefore x-\hat qy &\le& x-\hat qy_{n-1}R^{n-1}\\ &\le& x - (x_nR^n + x_{n-1}R^{n-1} - y_{n-1}R^{n-1} + R^{n-1})\\ &=& \left(\sum_{j=0}^{n-2}x_jR^j - R^{n-1}\right) + y_{n-1}R^{n-1}\\ &\lt& y_{n-1}R^{n-1} \le y \end{eqnarray} \]
となるので、$\hat q \ge q$ が保証される。
$\hat q \le q + 2$ の確認
前準備の計算をしておく。 \[ \hat q \leq \frac{x_nR+x_{n-1}}{y_{n-1}} = \frac{x_nR^n+x_{n-1}R^{n-1}}{y_{n-1}R^{n-1}} \le \frac{x}{y_{n-1}R^{n-1}} \lt \frac{x}{y - R^{n-1}} \]
ここで $\hat q \ge q+3$ を仮定すると \[ 3 \le \hat q - q \lt \frac{x}{y - R^{n-1}} - \frac{x}{y} + 1 = \frac{x}{y}\frac{R^{n-1}}{y - R^{n-1}} + 1 \] \[ \therefore \frac{x}{y} \gt 2\frac{y-R^{n-1}}{R^{n-1}} \ge 2(y_{n-1}-1) \]
ということで $2(y_{n-1}-1)\le \lfloor x/y \rfloor = q \le \hat q - 3 \le R-4$ より $y_{n-1} \lt R/2$ ということになるので、その対偶を取って \[ y_{n-1} \ge \frac{R}{2} \Rightarrow \hat q \le q + 2 \]
が導かれる。
2 ワード $\div$ 1 ワードの計算
理論的に求められることが分かったので、次は実践編。 計算機では 1 ワードの範囲を超える計算が直接的にはできないので 各ワードを半分のサイズに分けて 4 ハーフワード $\div$ 2 ハーフワード という形で計算する。
constexpr uint64 kHalfBits = 32;
constexpr uint64 R = 1ULL << kHalfBits;
constexpr uint64 kHalfMask = R - 1;
// 0<=y1,y0,x0<R, 0<=x12<R*R
uint64 DivHalf(uint64 x12, uint64 x0, uint64 y1, uint64 y0) {
uint64 q = x12 / y1;
if (q >= R)
q = R - 1;
uint64 r = x12 - q * y1;
while (q * y0 > r * R + x0) { // Loop at most twice.
--q;
r += y1;
}
return q;
}
uint64 Div(uint64 xh, uint64 xl, uint64 y, uint64* rem) {
uint64 x1 = xl >> kHalfBits, x0 = xl & kHalfMask;
uint64 y1 = y >> kHalfBits, y0 = y & kHalfMask;
uint64 q1 = DivHalf(xh, x1, y1, y0);
uint64 xr = xh * R + x1 - q1 * y;
uint64 q0 = DivHalf(xr, x0, y1, y0);
if (rem)
*rem = xr * R + x0 - q0 * y;
return q1 * R + q0;
}
ちなみにこの計算でのワードを $n$ ワードに置き換え、 途中計算にも多倍長整数を用いて再帰的に適用させると $2n$ ワード $\div n$ ワードの多倍長整数の除算を効率的に行うことができる。