Int 冪乗・モジュラ演算
概要
calx の Int は以下の冪乗・モジュラ演算を提供する:
pow(base, exp)— 冪乗 $a^e$ を計算するpowMod(base, exp, mod)— 冪剰余 $a^e \bmod m$ を計算するinverseMod(a, m)— モジュラ逆元 $a^{-1} \bmod m$ を計算する
いずれも繰り返し二乗法 (binary exponentiation) を核とするアルゴリズムである。 指数 $e$ のビット表現を走査しながら「二乗」と「乗算」を組み合わせることで、 $O(\log e)$ 回の乗算で冪乗を計算する。
calx では operator^ は冪乗演算子として定義されている。
C++ 標準では ^ はビット XOR だが、
calx の Int ではビット XOR には bitwiseXor() を使用する。
API の使い方については Int API リファレンス — 冪乗 および Int API リファレンス — モジュラ演算 を参照。
繰り返し二乗法 (Binary Exponentiation)
繰り返し二乗法は冪乗計算の基本アルゴリズムである。 指数 $e$ を 2 進表現し、各ビットに対して「二乗」と「乗算」を繰り返す。
基本アイデア
指数 $e$ を 2 進数で $e = (b_{k-1}\, b_{k-2}\, \cdots\, b_1\, b_0)_2$ と表すと、 冪乗は次のように分解できる:
$$a^e = a^{\sum_{i} b_i \cdot 2^i} = \prod_{b_i = 1} a^{2^i}$$
$a^{2^i}$ は $a$ を $i$ 回自乗するだけで得られるため、 $e$ を素朴に $e$ 回乗算する $O(e)$ の方法に比べて $O(\log e)$ 回の乗算で済む。
Left-to-right 走査
pow() は left-to-right (MSB から LSB へ) の走査を採用している。
アルゴリズムは以下の通り:
- $r \leftarrow 1$ で初期化
- 指数 $e$ の最上位ビットから最下位ビットへ順に走査する
- 各ビットについて: $r \leftarrow r^2$ (二乗)
- そのビットが 1 なら: $r \leftarrow r \times a$ (底を乗算)
計算例: $2^{10}$
$10 = (1010)_2$ なので、4 ビットを MSB から走査する:
| ステップ | ビット | 操作 | $r$ の値 |
|---|---|---|---|
| 1 | 1 | $r = 1^2 \times 2$ | $2$ |
| 2 | 0 | $r = 2^2$ | $4$ |
| 3 | 1 | $r = 4^2 \times 2$ | $32$ |
| 4 | 0 | $r = 32^2$ | $1024$ |
結果: $2^{10} = 1024$。乗算回数は最大 $2 \lfloor \log_2 e \rfloor$ 回である。
計算量
計算量: $O(\log e \cdot M(n))$
ここで $e$ は指数、$M(n)$ は $n$ ワードの乗算コストである。 冪乗では中間結果が急速に膨張する点に注意が必要である。 $a^e$ の桁数は概ね $e \times \log_{10}(a)$ 桁であり、 たとえば $a = 10^{100}$ (100 桁)、$e = 1000$ ならば結果は約 100,000 桁になる。 各ステップで中間値のサイズが増大するため、後半のステップほど乗算コストが重くなる。
冪剰余 (Modular Exponentiation)
powMod(a, e, m) は $a^e \bmod m$ を計算する。
暗号・数論アルゴリズムで最も頻繁に使われる演算の一つである。
アルゴリズム: Montgomery 冪剰余
$m$ が奇数の場合、calx は Montgomery 乗算ベースの冪剰余を使用する。 Montgomery 乗算は除算を使わずに剰余乗算を実現する手法で、 通常の $a \times b \bmod m$ に必要な除算を、シフトと加減算に置き換える。
Montgomery 乗算の核心は、$R = 2^{64n}$（$n$ はワード数）として 「Montgomery 表現」$\tilde{a} = a \cdot R \bmod m$ を使うことである。 Montgomery 乗算 $\text{Mont}(\tilde{a}, \tilde{b})$ は $\tilde{a} \cdot \tilde{b} \cdot R^{-1} \bmod m$ を除算なしで計算する (REDC)。
Sliding window 法
指数の走査には left-to-right sliding window 法を使用する。 指数のビット列をウィンドウ幅 $w$ で走査し、 奇数べきのテーブル $\{a^1, a^3, a^5, \ldots, a^{2^{w-1}-1}\}$ を事前計算する。 ウィンドウ幅は指数のビット長に応じて動的に選択される:
- $\leq 64$ ビット: $w = 3$
- $\leq 256$ ビット: $w = 4$
- $\leq 1024$ ビット: $w = 5$
- $\leq 4096$ ビット: $w = 6$
- $> 4096$ ビット: $w = 7$
MASM 特化カーネル
$n = 2, 4, 8$ ワードの Montgomery 乗算・二乗には専用の MASM アセンブリカーネル
(mpn_mont_mul_4, mpn_mont_sqr_4 等) を使用する。
MULX/ADCX/ADOX 命令でデュアルキャリーチェーンを形成し、
レジスタ常駐でメモリアクセスを最小化する。
$n \leq 8$ では全バッファをスタック上に配置し、ヒープ割り当てゼロで動作する。
中間値の爆発を防ぐ
冪剰余の核心は、各ステップで $\bmod m$ を取ることにある。 素朴に $a^e$ を先に計算してから $\bmod m$ を取ると、 中間値が天文学的な桁数に膨れ上がってしまう。 Montgomery 乗算では REDC により各ステップで自動的に $m$ 未満に縮小される。
負の指数
$e < 0$ の場合、powMod は以下のように計算する:
$$a^{-|e|} \bmod m = (a^{|e|} \bmod m)^{-1} \bmod m = \text{inverseMod}(a^{|e|} \bmod m,\; m)$$
すなわち、まず $a^{|e|} \bmod m$ を計算し、その結果のモジュラ逆元を求める。
応用
- RSA 暗号: 暗号化 $c = m^e \bmod n$、復号 $m = c^d \bmod n$
- Diffie-Hellman 鍵交換: $g^a \bmod p$ の計算
- Miller-Rabin 素数判定: $a^d \bmod n$ の計算 (calx の
isPrime()で使用)
計算量: $O(\log e \cdot M(n))$、ここで $n$ は $m$ のワード数
モジュラ逆元 (inverseMod)
inverseMod(a, m) は $a \cdot x \equiv 1 \pmod{m}$ を満たす $x$ を求める。
逆元が存在するための必要十分条件は $\gcd(a, m) = 1$ である。
戦略 1: 拡張ユークリッド互除法 (デフォルト)
is_prime = false (デフォルト) の場合、拡張ユークリッド互除法を用いる。
拡張ユークリッド互除法は $\gcd(a, m) = a \cdot x + m \cdot y$ となる $x, y$ を求めるアルゴリズムである。 $\gcd(a, m) = 1$ のとき $a \cdot x + m \cdot y = 1$ であり、 両辺を $\bmod m$ すると $a \cdot x \equiv 1 \pmod{m}$ が得られる。
この方法は $m$ が合成数でも素数でも動作し、最も汎用的である。 計算量は $O(\log^2 m)$ (GCD と同じ) であり、小さな $m$ では非常に高速である。
戦略 2: フェルマーの小定理 (is_prime = true)
is_prime = true を指定すると、フェルマーの小定理を用いる:
$$a^{p-1} \equiv 1 \pmod{p} \quad \Rightarrow \quad a^{-1} \equiv a^{p-2} \pmod{p}$$
$m$ が素数であることが既知の場合、$a^{m-2} \bmod m$ を powMod で計算するだけでよい。
GCD 計算を回避でき、非常に大きな $m$ では NTT 乗算を活用した powMod のほうが
拡張ユークリッド互除法より高速になる場合がある。
注意: $m$ が素数でないのに is_prime = true を指定すると、
誤った結果が返される。呼び出し側が素数性を保証する必要がある。
前提条件
$\gcd(a, m) = 1$ でなければ逆元は存在しない。 この条件が満たされない場合の動作は未定義である。
関連記事: 多倍長整数の高度なアルゴリズム
mod の正規化
C++ の % 演算子は被除数の符号を保存する。
たとえば -7 % 5 = -2 である。
これは数学的な剰余の定義とは異なる。
calx の IntModular::mod は常に非負の剰余を返す:
$$\text{mod}(a, m) = r \quad \text{ただし} \quad r \in [0, |m|)$$
| 式 | C++ の % | calx の mod |
|---|---|---|
| $-7 \bmod 5$ | $-2$ | $3$ |
| $7 \bmod 5$ | $2$ | $2$ |
| $-7 \bmod{-5}$ | $-2$ | $3$ |
数学では剰余は常に $[0, |m|)$ の範囲に正規化される。
calx はこの数学的定義に従うことで、
モジュラ演算の結果が常に非負であることを保証する。
これにより powMod や inverseMod の結果も常に非負になる。
設計判断
pow で left-to-right を採用する理由
pow() は left-to-right 走査を使用する。
left-to-right では結果が大きい値から構築されるため、
メモリアクセスパターンがやや局所的になり、キャッシュ効率が向上する。
また、right-to-left では底と結果の両方を逐次更新する必要があるが、
left-to-right では結果のみを更新すればよい。
powMod で Montgomery + sliding window を採用する理由
powMod() は Montgomery 乗算と left-to-right sliding window 法を使用する。
Montgomery 乗算は除算を回避し、REDC (Montgomery 縮小) により各ステップの剰余演算を
シフトと加減算で実現する。
Sliding window 法は固定ウィンドウ法やナイーブな繰り返し二乗法より乗算回数が少ない。
$m$ が偶数の場合は Montgomery 表現が使えないため、 従来の繰り返し二乗法 (right-to-left) にフォールバックする。
n=2,4,8 特化カーネルの効果
256 bit (n=4) の冪剰余では Montgomery ループが数百回実行されるため、 1 回あたりの mont_mul/mont_sqr のコストが直接的に全体性能に影響する。 n=4 特化の MASM カーネルは汎用パスより関数呼び出しオーバーヘッドが小さく、 レジスタ常駐で L1 キャッシュミスを回避する。
inverseMod で 2 つの戦略を提供する理由
拡張ユークリッド互除法（HGCD ベース）は汎用的だが、
フェルマーの小定理による方法は $m$ が素数のとき GCD 計算を完全に回避できる。
特に非常に大きな素数 $m$ に対しては、
powMod が NTT 乗算を活用するため効率的である。
両方を提供することで、用途に応じた最適な選択が可能になる。
参考文献
- Montgomery, P. L. (1985). "Modular multiplication without trial division". Mathematics of Computation, 44(170), 519–521.
- Knuth, D. E. The Art of Computer Programming, Volume 2: Seminumerical Algorithms. 3rd ed. Addison-Wesley, 1997. §4.6.3.
- Cohen, H. A Course in Computational Algebraic Number Theory. Springer, 1993.