複数のモジュラ逆数を求める効率的な方法

こちらの記事を読んで。

内容まとめ

モジュラ逆数を求めるのに、フェルマーの小定理を用いて $a^{-1} \equiv a^{M-2} \mod M$ ですると、モジュラ乗算が $\log M$ 回必要である。
拡張ユークリッドの互除法を用いるとこれより速いがある程度かかる。その回数を仮に $K$ とする。

$N$個の数 $x_1, \dots, x_N$ に対してモジュラ逆数を求めるのに、単純に繰り返すと乗算は $NK$ 回になる。
ずっと高速に求められる Peter Montgomery のトリックを紹介している。

全ての積 $\prod x_i$ に対するモジュラ逆数 $I$ を求める。$I = 1/x_1 \cdot \ldots \cdot x_N$ 乗算は $N + K$ 回。

この $I$ を用いると、任意の $x_i$ のモジュラ逆数は $x_i^{-1} = I \cdot (x_1 \cdot \ldots \cdot x_{i-1}) \cdot (x_{i+1} \cdot \ldots \cdot x_N)$ である。

これを一つずつ求めるのではなく、もう一工夫する。
最後の項 $x_{i+1} \cdot \ldots \cdot x_N$ は最初に $\prod x_i$ を求めるときの中間値として（後ろから求める）中間値を全て残しておき再利用する。追加の計算量はない。
前2項 $I \cdot (x_1 \cdot \ldots \cdot x_{i-1})$ は、初期値を $I$ として累積積で計算できる。乗算は $N$ 回。
そして対応する左右を掛けることで、$N$ 個全ての値に対するモジュラ逆数が、総乗算回数 $3N+K$ で得られる。

普通 $3 \ll K$ なので、高速化される。

検証

Montgomery のトリックを実装した。

modRecips xs = zipWith mul cs rs
  where
    c : rs = scanr mul 1 xs
    cs = scanl mul (modRecip c) xs

結果を比較して正しいことを確認する。

import Test.QuickCheck

-- フェルマーの小定理による
naive1 = map (\x -> modPower x (modBase - 2))

-- 拡張ユークリッドの互除法による
naive2 xs = map modRecip xs

-- 同値性テスト
prop :: [Positive Int] -> Bool
prop ps = modRecips xs == naive1 xs
  where
    xs = map (reg . getPositive) ps

ghci> quickCheck prop
+++ OK, passed 100 tests.
(0.09 secs, 51,122,488 bytes)

実行時間を比較してみる。

tryit f n = do
  v <- vectorOf n arbitrary :: Gen [Positive Int]
  let xs = map (reg . getPositive) v
  let rs = f xs
  return $ all (1 ==) $ zipWith mul xs rs

なお modBase = 1_000_000_007 である。

ghci> quickCheck $ tryit naive1 10000
+++ OK, passed 100 tests.
(29.29 secs, 20,911,178,216 bytes)
ghci> quickCheck $ tryit naive2 10000
+++ OK, passed 100 tests.
(4.02 secs, 3,493,459,056 bytes)
ghci> quickCheck $ tryit modRecips 10000
+++ OK, passed 100 tests.
(1.78 secs, 1,674,114,352 bytes)

なるほど速いね。

おわりに

どんな値について逆元が必要なのかが事前に分かり、しかも現実的な個数なときにのみ適用できる方法な感じ。
競プロではあまりそういう感じではないからなぁ。