Python で mod の下での逆元を計算するテーブルを作成する

wakabame.hatenablog.com
では, ${}_nC_r$ の素数 $p$ で割った余りを求める方法について解説しました.

実は, この計算アルゴリズムは高速化できます.

${}_nC_r$ の計算結果が何度も呼び出される場合を考えましょう. そのたびに

前者は予め必要な大きさのリストを構築して, 計算結果を最初に持っておけば $O(n)$ で済みます.

一方で, 後者のボトルネックについては, 繰り返し二乗法を全ての整数に対して行わなければならないのが（定数倍とはいえ）しんどいです.

そこで, 次のブログを参考に逆元テーブルの構築を考えます.

以後, $P$ を十分大きな素数とします.

P = 10**9 + 7
N = 10000
inv_t = [0]+[1]
for i in range(2,N):
  inv_t += [inv_t[P % i] * (P - int(P / i)) % P]

これにより, $i=1,2,\cdots , N$ に対して inv_t[i] は" $i$ の法 $P$ の下での逆元"を表します.

$a=0,1, \cdots$ に対して, $a$ の法 $P$ の下での逆元 $a^{-1}$ を求めましょう.
まず,

$x =0,1, \cdots , P-1$ が $a$ の法 $P$ の下での逆元であるとは,
ある $y_a \in \mathbb{N}$ が存在して,
$ax+Py_a =1$
を満たすこと

と定式化できます.
そのような $x$ は $a$ が $P$ の倍数でない限り存在し一意であり, 以後これを $a^{-1}$ と表します.

さて, アルゴリズムを説明していきます.

まず最初に, $0,1$ に対する逆元はそれぞれ $0,1$ と定めてしまいます.

以後, $a$ は $2$ 以上とし, $a$ 未満の逆元は既に得られたとしましょう.

$P$ を $a$ で割り, その商, 余りをそれぞれ $q$ , $r$ とおきます. すなわち

$P = qa+r$

が成り立ちます.

このとき, $r$ は $a$ 未満であることに注意します.

両辺に $a^{-1}$ をかけると,

$P a^{-1} =qa a^{-1} +ra^{-1}$

であり, $a$ の逆元の定義から,

$P a^{-1} =q- Py_a +ra^{-1}$

が成り立ちます.

移項して整理すると,

$ra^{-1} =-q + P(y_a -a^{-1} )$

です.

さらに $r^{-1}$ をかけると,

$a^{-1} -Py_r = -q r^{-1 }+ P(y_a -a^{-1} )r^{-1}$

が成り立ちます.

両辺, $P$ で割った余りに着目すると,

$a^{-1} = -q r^{-1 } \%P$

を得ます.

$r^{-1}$ は既に得られているので, この方法で小さい数から逐次的に逆数を求めればよく, トータルの計算量は $N$ であることがわかりました.