円周率.jp (http://xn--w6q13e505b.jp/method/div.html)

除算

$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 >> kHlafBits, x0 = xl & kHalfMask;
  uint64 y1 = y >> kHlafBits, 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$ ワードの多倍長整数の除算を効率的に行うことができる。