Int 平方根アルゴリズム

概要

整数 $n \geq 0$ に対して、$s = \lfloor \sqrt{n} \rfloor$ を計算する。 すなわち、$s^2 \leq n < (s+1)^2$ を満たす最大の整数 $s$ を求める。

内部は mpn レベルで次のように構成される:

$$\boxed{\text{Int::sqrtRem}} \;\xrightarrow{\text{mpn}}\; \begin{cases} \text{sqrtrem1} & \text{(1 ワード)} \\ \text{sqrtrem2} & \text{(2 ワード)} \\ \text{dc\_sqrtrem} & \text{(3+ ワード)} \end{cases}$$

API の使い方については Int API リファレンス — 平方根 を参照。

アルゴリズム選択

calx は入力のワード数に応じてアルゴリズムを自動選択する。 すべてのサイズで Zimmermann の分割統治法 (dc_sqrtrem) がベースとなっており、 1–2 ワードの特殊ケースは専用の高速ルーチンで処理する。

ここで $M(n)$ は $n$ ワードの乗算コスト、1 ワード = 64 bit である。

sqrtrem1 (1 ワード平方根)

正規化された 1 ワード ($a_0 \geq 2^{62}$) の平方根と余りを計算する。 除算を一切使わず、逆平方根テーブル引き + Newton 反復 2 回で根を求める。

アルゴリズム

1. テーブル引き: 384 エントリの事前計算済みテーブル invsqrttab から $1/\sqrt{a_0}$ の 8 ビット近似値を取得する。
2. Newton 反復 1: 逆平方根の精度を $\sim$16 ビットに倍増する。
3. Newton 反復 2: 逆平方根の精度を $\sim$32 ビットに倍増する。
4. 正方向の根: $s \approx a_0 \cdot (1/\sqrt{a_0})$ で平方根を得て、余りから最終補正を行う。

逆平方根 $1/\sqrt{a_0}$ を求めてから $s = a_0 \cdot (1/\sqrt{a_0})$ で変換するのは、 Newton 反復に除算が不要になるためである ($x_{k+1} = x_k (3 - a x_k^2) / 2$ は乗算のみ)。

sqrtrem2 (2 ワード平方根)

正規化された 2 ワード ($n_1 \geq 2^{62}$) の平方根を計算する。 sqrtrem1 を上位ワードに適用し、1 回の除算ステップで精度を拡張する。

1. sqrtrem1(n_1) で上位ワードの平方根 $s_0$ と余り $r_0$ を得る
2. $[r_0, n_0]$ を $2s_0$ で割って商 $q$ を得る
3. $s = s_0 \cdot 2^{32} + q$ が 2 ワードの平方根候補
4. 余りから $q^2$ を引き、負なら $s$ を 1 減らして補正する

結果は 64 ビットの根であり、in-place 計算をサポートする。

Karatsuba 平方根 (Zimmermann)

Paul Zimmermann (1999) が発表した再帰的平方根アルゴリズムである。 平方根の計算を乗算 1 回と同じオーダー $O(M(n))$ に帰着させる。

基本的なアイデア

$n$ ワードの入力 $N$ を 4 つのブロックに分割する:

$$N = a_3 B^3 + a_2 B^2 + a_1 B + a_0, \quad B = 2^{n/4 \cdot 64}$$

アルゴリズムは以下の 7 ステップで進む:

1. 再帰: 上位半分 $[a_3, a_2]$ に対して再帰的に dc_sqrtrem を呼び出し、 平方根 $s'$ と余り $r'$ を計算する
2. 余りキャリー調整: 再帰の余りキャリーを処理する
3. 除算: 余り $r'$ と $a_1$ を連結した値を $2s'$ で除算 (Burnikel-Ziegler) し、商 $q$ と余り $u$ を得る
4. 商の抽出: $q$ の上位ビットを処理し、$s = s' \cdot B + q_{\text{low}}$ を構成する
5. 奇数商の補正: $q$ が奇数の場合、余りを調整する
6. 自乗減算: 余り $r = u \cdot B + a_0 - q_{\text{low}}^2$ を計算する
7. 負余り補正: $r < 0$ なら $r \mathrel{+}= 2s - 1$、$s \mathrel{-}= 1$ で補正する

sqrtRem の役割

Zimmermann アルゴリズムは平方根 $s$ と余り $r$ を同時に返す:

$$N = s^2 + r, \quad 0 \leq r \leq 2s$$

余りを返すことで、再帰呼び出しの結果を外側のステップに直接利用できる。 この「平方根と余り」のペアを返す設計が再帰を効率的にする鍵である。

計算量の直観

再帰の各レベルでは、サイズ $n$ の除算 ($O(M(n))$) とサイズ $n/2$ の自乗 ($O(M(n/2))$) を行う。 再帰のコスト $T(n)$ は:

$$T(n) = T(n/2) + O(M(n))$$

$M(n) \geq cn$ (少なくとも線形) であれば、Master theorem により $T(n) = O(M(n))$ となる。 すなわち、平方根は乗算 1 回と同じコストで計算できる。

isSquare (完全平方数判定)

$n$ が完全平方数 ($n = k^2$ となる整数 $k$ が存在する) かどうかを判定する。 単純に sqrtRem を呼んで余りがゼロか確認するのは $O(M(n))$ のコストがかかるため、 calx では GMP 方式の3 層フィルターで非平方数を高速に棄却する。

フィルター構成

第 1 層: mod 256 フィルター

平方数 $k^2$ を $256$ で割った余りは、$256$ 通りのうち $44$ 通りしか取りえない。 入力の最下位バイトをビットマスクテーブルで引くだけで、 非平方数の $82.81\%$ を $O(1)$ で棄却できる。

第 2 層: 小素数剰余フィルター

mod_34lsub1 により $n \bmod (2^{48} - 1)$ を高速に計算し、 $3, 5, 7, 13, 17, 97$ 等の小さい素数に対する二次剰余判定を行う。 これらは $2^{48} - 1$ の因数であるため、追加の剰余演算なしに判定できる。

第 3 層: 完全検算

フィルターを通過した入力に対してのみ sqrtRem を呼び出し、 余り $r = 0$ かどうかを確認する。 全 3 層を合わせた棄却率は $99.62\%$ であり、 非平方数に対する平均コストは $O(n)$ に抑えられる。

参考文献

• Zimmermann, P. (1999). "Karatsuba Square Root". INRIA Research Report RR-3805.
• Brent, R. P. and Zimmermann, P. (2010). Modern Computer Arithmetic. Cambridge University Press. Chapter 1.5: Square Root.
• GMP Documentation. "Integer Roots". https://gmplib.org/manual/