任意 mod で二項係数を列挙する


(追記 2020/11/30) :

計算量解析がより厳密になりました. 以前は $\displaystyle O\left(\frac{N(\log N)^2}{\log\log N}\right)$ と書いていましたが, $O(N(\log N)(\log\log N))$ であることが分かりました.

また, 本記事の方法を改善し, $\displaystyle O\left(\frac{N\log M}{\log\log M}\right)$ での列挙について suisen-kyopro.hatenablog.com に書いています. 合わせて読んでみて下さい.

(追記終わり)

$N=10^5$, $M=10^9$ くらいの場合に, 二項係数 ${}_NC_0, {}_NC_1, \dots, {}_NC_N \bmod M $ を $O(N(\log N)(\log\log N))$ で列挙します. 以降, 表記の簡単のため, 二項係数を $C(N,i)$ のように表すことにします.

\(M\) が素数の場合は, 「1000000007 で割ったあまり」の求め方を総特集! 〜 逆元から離散対数まで 〜 - Qiita で紹介されているように, 任意の \(1\leq x\lt M\) に対して \(x\cdot x'\equiv 1 \pmod M\) を満たす $x'$ ($x$ の \(\bmod M\) での乗法逆元) が存在するので, 計算は容易です.

しかし, \(M\) が素数でない場合は, 逆元が存在するとは限りません. そこで, 別の方法を考えます.

まず, $C(N,i)$ は式 $(1)$ の漸化式で求まるので, 初期値を $X=1$ として, $i=1,2,\dots N$ の順に $X=C(N,i)$ を計算します.

$$
C(N,i)=\begin{cases}\displaystyle 1& \text{(if $i=0$)} \\\displaystyle\frac{N-i+1}{i}\cdot C(N,i-1)& \text{(otherwise)}\end{cases}\tag{1}
$$

ここで, 素因数分解の形で $X$ を管理することにします. つまり, 素数 $p_1, p_2, \dots$ を用いて, $X={p_1}^{q_1}{p_2}^{q_2}\cdots$ の形で保持します. 素数は無限個ありますが, $C(N,i)$ が持つ素因数は $N$ 以下なので, $N$ 以下の素数を列挙しておけば十分です. 以下では, $N$ 以下の素数の個数を $K$ とします.

まず, 素因数分解の指数 $q_1, q_2,\dots, q_K$ の更新を考えます. $N-i+1={p_1}^{a_1}{p_2}^{a_2}\cdots {p_K}^{a_K}$, $i={p_1}^{b_1}{p_2}^{b_2}\cdots {p_K}^{b_K}$とすると, $q_j\leftarrow q_j+a_j-b_j$ により更新すればよいです.

ただし, $\displaystyle K=O\left(\frac{N}{\log N}\right)$ なので, 全ての $j$ に対して毎回更新していると計算量が $\displaystyle O\left(\frac{N^2}{\log N}\right)$ となってしまいます. そこで, $a_j\neq 0$ または $b_j\neq 0$ であるような $j$ のみを選んで更新することにします. 定数倍高速化をしているだけに見えるかもしれませんが, 実はこのような $j$ の個数は合計で $\displaystyle\sum_{j=1}^K\frac{N}{p_j}=O(N\log\log N)$ *1 程度しかありません.

エラトステネスの篩において, 自身を割り切る素数を覚えておくようにすれば素因数分解は $O(\log N)$ で可能 *2 なので, 結局 $q$ の更新は全体 $O(N\log N)$ で出来ます.

素因数分解が求まっても, 元の値が高速に求まらなければ意味がありません. 例えば, 毎回愚直に $X={p_1}^{q_1}{p_2}^{q_2}\cdots {p_K}^{q_K}$ を計算していると, 累乗が $O(1)$ で計算できたとしても全体で $\displaystyle O\left(\frac{N^2}{\log N}\right)$ となってしまいます.

いま, 求めたいのは $X={p_1}^{q_1}{p_2}^{q_2}\cdots {p_K}^{q_K}$ です. $p_i$ は変化せず, $q_i$ は高速に更新できることが分かっています. そこで, $r_i={p_i}^{q_i}$ とし, これを \(\bmod M\) 上の乗算 $\times$ を二項演算とするセグメント木に載せます. すると, \(X=r_1\times r_2\times \cdots \times r_K\) は $O(\log K)=O(\log N)$ で計算できます. 更新回数は合計で $O(N\log\log N)$ 回だったので, セグメント木パートの計算量は全体で $O(N(\log N)(\log\log N))$ です.

最後に, $r_i={p_i}^{q_i}$ の更新にかかる計算量を確認します. 二分累乗法により $O(\log q_i)$ で計算できますが, $q_i$ はどれくらい大きくなるでしょうか. $X$ は $N!$ を割り切るので, $N!$ を素因数分解した場合の $q_i$ の大きさを評価してみます. $N!$ が持つ素因数 $p_i$ の個数は, $\displaystyle \sum_{j=1}^{\infty} \left\lfloor\frac{N}{{p_i}^j}\right\rfloor$ 個であることが知られており *3,
$$\sum_{j=1}^{\infty} \left\lfloor\frac{N}{{p_i}^j}\right\rfloor\leq\sum_{j=1}^{\infty}\frac{N}{{p_i}^j}=\frac{N}{p_i-1}=O(N)$$
より, 累乗は $O(\log N)$ で計算できます. $r_i$ の更新回数は $O(N\log\log N)$ 回なので, 累乗計算パートも全体で $O(N(\log N)(\log\log N))$ であることが示せました.

以上より, 全体 $O(N(\log N)(\log\log N))$ で任意 mod での $C(N,i)$ を列挙することが出来ました.

実装↓

github.com

*1:エラトステネスの篩の計算量解析と同様です. $N$ 以下の素数の逆数和は $O(\log\log N)$ であることから導かれます

*2:全ての素因数は $2$ 以上なので, 割るごとに値は半分以下になります. 従って, 結局この割り算は高々 $\lceil \log_2{N} \rceil$ 回で終わります.

*3:参考 : ルジャンドルの定理(階乗が持つ素因数のべき数) | 高校数学の美しい物語