Int 平方根・n 乗根

概要

Int クラスは以下の根号関連の演算を IntSqrt の静的メソッドとして提供する:

• 整数平方根: $\lfloor \sqrt{n} \rfloor$ を計算する
• 整数 n 乗根: $\lfloor n^{1/k} \rfloor$ を計算する
• 完全平方数判定: $n$ が完全平方数かどうかを高速に判定する (isSquare)

いずれの演算もオペランドのサイズに応じて最適なアルゴリズムが自動選択される。 ユーザーが明示的にアルゴリズムを選ぶ必要はない。

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

平方根アルゴリズム選択

calx は入力のサイズに応じて 3 段階のアルゴリズムを自動選択する。

ここで $M(n)$ は $n$ ワードの乗算コストを表す。 1 ワード = 64 bit であり、17 ワードは約 1,088 bit (約 327 桁) に相当する。

52 bit 以下の入力はハードウェアの倍精度浮動小数点演算で正確に計算できるため、 std::sqrt を直接使用する。

Babylonian 法 (Newton 法)

古代バビロニアに起源を持つ、平方根の古典的な反復法である。 これは $f(x) = x^2 - n$ に対する Newton 法そのものであり、次の漸化式で収束する:

$$x_{k+1} = \frac{x_k + n / x_k}{2}$$

2 次収束: 各反復で正しい桁数が倍になる。 たとえば 4 桁の精度があれば、次の反復で 8 桁、その次で 16 桁の精度が得られる。

初期近似値

収束速度は初期近似値の精度に大きく依存する。 calx では入力の最上位ワードから浮動小数点の std::sqrt で初期近似値を計算する。 これにより 52 bit 程度の精度を持つ初期値が得られ、反復回数を最小化する。

整数平方根の終了条件

実数の Newton 法は $|x_{k+1} - x_k|$ が十分小さくなれば停止するが、 整数平方根では異なる終了条件が必要である。 $x_{k+1} \geq x_k$ となった時点で反復を停止し、$x_k$ が $\lfloor \sqrt{n} \rfloor$ となる。 これは列 $\{x_k\}$ が $\lfloor \sqrt{n} \rfloor$ に対して単調減少することに基づく。

計算量

計算量: $O(n \cdot M(n))$

一見すると各反復は $n$ ワードの除算 1 回であり、 反復回数は $O(\log n)$ 回なので $O(M(n) \log n)$ に見える。 しかし除算のコストが $M(n)$ であり、しかも中間値のサイズが反復ごとに変化するため、 厳密な解析では $O(n \cdot M(n))$ となる。 これは Karatsuba 平方根の $O(M(n))$ と比べて漸近的に遅いが、 小さい入力ではアルゴリズムの単純さと定数倍の軽さが有利に働く。

Karatsuba 平方根 (Zimmermann)

Paul Zimmermann が 1999 年に発表した再帰的平方根アルゴリズムである。 平方根の計算を乗算サイズの演算に帰着させ、$O(M(n))$ という最適な漸近計算量を達成する。 これは乗算 1 回と同じオーダーで平方根が計算できることを意味する。

基本的なアイデア

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

$$n = a_3 \cdot B^3 + a_2 \cdot B^2 + a_1 \cdot B + a_0 \quad (B = 2^{16n})$$

上位半分 $a_3 \cdot B + a_2$ に対して再帰的に平方根を計算し、 その結果を用いて乗算サイズの除算 1 回で全体の平方根に精密化する。

sqrtRem

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

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

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

計算量

計算量: $O(M(n))$ -- 乗算 1 回と同じ漸近コスト

17 ワード以上の入力に対して使用される。 Babylonian 法の $O(n \cdot M(n))$ と比較すると、入力サイズが大きくなるほど差は広がる。

完全平方数判定 (isSquare)

$n$ が完全平方数かどうかを判定する。 素朴な方法は平方根を計算して検証することだが、 大きな整数の平方根は高コストな演算である。 calx ではモジュラーフィルタリングにより、 平方根を計算することなく非平方数の大部分を棄却する。

モジュラーフィルタ

完全平方数は任意の法 $m$ に対して平方剰余でなければならない。 この性質を利用して、$n \bmod m$ が平方剰余でないことを検出すれば、 平方根を計算せずに非平方数と判定できる。

calx では GMP の perfsqr.c に学んだ 3 層フィルタを適用する:

1. Filter 1: mod 256 — 最下位 1 バイトの平方剰余テーブルで判定。 $O(1)$ で約 82% の非平方数を棄却する。
2. Filter 2: mod_34lsub1 — $2^{24} - 1 = 16{,}777{,}215$ のリム和を計算し、 小素数 (3, 5, 7, 11, 13 等) で平方剰余をチェック。 $O(n)$ で追加 ~97.8% を棄却する。
3. Filter 3: sqrtRem — フィルタを通過した候補のみ、 実際の平方根を計算して $s^2 = n$ を検証する。

3 層のフィルタを組み合わせることで、非平方数の 99.9% 以上が 平方根を計算することなく棄却される。

mod 256 テーブルは 32 バイト (ビットテーブル) でキャッシュ内に常駐し、 mod_34lsub1 は多倍長整数の全リムを 1 回走査するだけで複数の小素数による判定を同時に実行する。

n 乗根

整数 $A$ の $n$ 乗根 $\lfloor A^{1/n} \rfloor$ を計算する。 $f(x) = x^n - A$ に対する Newton 法を用いる:

$$x_{k+1} = \frac{(n-1) \cdot x_k + A / x_k^{n-1}}{n}$$

精度倍増法

Newton 法の 2 次収束の性質を活用し、反復ごとに作業精度を倍増させるテクニックを適用する。 初期の反復は低い精度で実行し、収束に従って精度を上げていく。

具体的には:

1. 最終的に必要な精度 $P$ を決定する
2. 必要な反復回数を $\lceil \log_2 P \rceil$ と見積もる
3. 最初の反復は精度 $P / 2^k$ で実行し、次は $P / 2^{k-1}$、...、最後は精度 $P$ で実行する

初期の反復は小さい数を扱うため軽く、 計算コストの大部分は最後の数回の反復に集中する。 この方式により、全反復を最終精度で行う場合と比べて総コストが大幅に削減される。

初期近似値

入力の最上位ワードから浮動小数点の $n$ 乗根を計算し、初期近似値とする。 これにより約 52 bit の精度を持つ初期値が得られ、 精度倍増法の最初のステップで十分な精度を確保できる。

計算量

計算量: $O(M(n) \cdot \log n)$

精度倍増法により、各反復のコストは幾何級数的に増加するが、 総コストは最後の反復のコストに支配されるため、 $O(M(n) \cdot \log n)$ に収まる。

nthRootRem

平方根と同様に、$n$ 乗根と余りのペアを返す nthRootRem も提供する:

$$A = s^n + r, \quad 0 \leq r \leq (s+1)^n - s^n - 1$$

設計判断

Zimmermann 閾値が 17 ワードの理由

Karatsuba 平方根は再帰的なアルゴリズムであり、分割・結合のオーバーヘッドがある。 17 ワード未満の入力では、Babylonian 法の単純な反復のほうが このオーバーヘッドを上回る速度を発揮する。 17 ワード (約 1,088 bit) は実測ベンチマークに基づいて決定された閾値である。

isSquare にモジュラーフィルタを使う理由

平方根の計算は $O(M(n))$ 以上のコストがかかるが、 モジュラーフィルタは $O(1)$ (ワード数に依存しない一定時間) で実行できる。 非平方数が圧倒的に多い一般的なユースケースでは、 フィルタにより平方根計算の呼び出し回数が 1000 分の 1 以下に削減されるため、 全体の性能が劇的に向上する。

参考文献

• Zimmermann, P. (1999). "Karatsuba Square Root". INRIA Research Report.
• Bernstein, D. J. "Detecting perfect powers in essentially linear time".
• Crandall, R. and Pomerance, C. Prime Numbers: A Computational Perspective. Springer, 2005.