[追記: 2021/04/19]
この記事が分かりにくいと思ったので Qiita の方に書き直しました．一応記事は残しますが，deprecated 的な感じでお願いします．
qiita.com
[追記終わり]
suisen-kyopro.hatenablog.com これを $\displaystyle O\left(\frac{N\log M}{\log\log M}\right)$ にします.
基本的なアイデアは素因数分解であり, 前記事と同様です. 式 $(1)$ を用いるのも同じです. しかし, 今回セグメント木は使いません.
$$
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}
$$
前回の記事では, \(M\) の素因数以外の素数についても指数を管理して累乗を計算していました. しかし, これには無駄な点が多いです. というのも, \(M^2\leq 2^{63}-1\) (符号付き 64 bit 整数の最大値) の範囲では, 実は \(M\) の素因数の個数は高々 $9$ つしかありません. *1
\(M\) と互いに素であれば \(\bmod M\) 上の乗算に関して逆元が存在するので, 除算が可能です. そこで, 「除算が可能な部分は逆元を使って楽に計算しよう」というのが今回のメインテーマです.
以下では, \(M\) の素因数は小さい順に \(p_1, \dots, p_K\) であるとします.
除算が可能な部分と不可能な部分を分けるため, \(X=C(N,i)\) を次のように分解します. $$X=Y\cdot {p_1}^{q_1}{p_2}^{q_2}\cdots {p_K}^{q_K}.$$ また, 上式において $Z={p_1}^{q_1}{p_2}^{q_2}\cdots {p_K}^{q_K}$ とします.
その他の整数 $A$ についても同様に $$A=(M \text{と素因数を共有しない部分}) \times (M \text{と素因数を共有する部分})$$ のように分割して考えます. 表記の簡単のため, $$\begin{align}S_M(A)&=(M \text{と素因数を共有しない部分}),\\ T_M(A)&=(M \text{と素因数を共有する部分})\end{align}$$ と定義しておきます *2. つまり, $S_M(X)=Y$, $T_M(X)=Z$ ということになります.
前置きはこれくらいにして, 本計算に入ります. $Z$ の指数部分 $q_i$ の更新方法は前記事と基本的に同様です. しかし, 今回は $p_1,\dots ,p_K$ 以外の素因数の指数は管理する必要がありません.
そこで, 式 $(2)$ を満たすような新しい配列 $d$ を用意します *3.
$$
d_i=\max \{p_j\mid p_j \text{は $i$ を割り切る} \} \cup\{0\}\quad (1\leq i\leq N) \tag{2}
$$
この $d_i$ は, 擬似的には以下のように求められます.
// primeFactors(M) は M の素因数を小さい順に列挙する for (int prime : primeFactors(M)) { for (int i = prime; i <= N; i += prime) { d[i] = prime; } }
計算回数は \(\displaystyle\sum_{p_i}\frac{N}{p_i}\) 回程度になりますが, これは次のように評価されます. $$\displaystyle\sum_{p_i}\frac{N}{p_i}\leq\sum_{p_i}N=NK=O\left(\frac{N\log M}{\log\log M}\right).\tag{3}$$
いま求めた $d$ を用いることで, 更新式に現れる $N-i+1$ および $i$ に対して, $S_M(N-i+1)$, $T_M(N-i+1)$, $S_M(i)$, $T_M(i)$ を効率的に求めることが出来ます. 以下にコードのイメージを示します.
int s1 = N - i + 1; // div は p_1,...,p_K のいずれか for (int div = d[s1]; d[s1] > 0; div = d[s1]) { int count = 0; do { s1 /= div; count++; } while (s1 % div == 0); // この時点で count は N-i+1 が div で割り切れる回数 (T_M(N-i+1) における div の指数) となっているので, Z の指数を更新する. // indexOf() は p_i -> i を表す仮の関数 q[indexOf(div)] += count; } // この時点で s1=S_M(N-i+1) となっているので, Y を更新する. Y *= s1; int s2 = i; for (int div = d[s2]; d[s2] > 0; div = d[s2]) { int count = 0; do { s2 /= div; count++; } while (s2 % div == 0); q[indexOf(div)] -= count; } // 実際は mod M での s2 の逆元を掛ける Y /= s2;
計算量を考えましょう. Y/=s2 の部分で逆元が必要になっていますが, 合成数 $\bmod$ の場合でも $1,\dots,N$ に対する逆元の列挙は $O(N)$ で可能です*4.
割り算パートの解析は難しそうですが, 全体で $\displaystyle O\left(\frac{N\log M}{\log\log M}\right)$ であることが示せます. 各 $p_i$ について, $1,2,\dots,N$ の中に素因数として合計いくつ含まれるかが分かれば良いです. これは, 有名な 「$N!$ は素数 $P$ で何回割り切れるか?」と同じ問題です. そして, この問題の答えは $\displaystyle \sum_{i=1}^{\infty}\left\lfloor\frac{N}{P^i}\right\rfloor $ です *5. 従って, 計算量は $$\sum_{p_i}\sum_{j=1}^{\infty}\left\lfloor\frac{N}{{p_i}^j}\right\rfloor\leq\sum_{p_i}\sum_{j=1}^{\infty}\frac{N}{{p_i}^j}=\sum_{p_i}\frac{N}{p_i-1}=O\left(\frac{N\log M}{\log\log M}\right) \tag{4}$$ です (最後の変形は式 $(3)$ と同じ理屈です).
最後に素因数分解形 ${p_1}^{q_1}{p_2}^{q_2}\cdots {p_K}^{q_K}$ から $Z$ の値を復元するパートを考えます. 基本的には以下に示すコードのように愚直に計算します.
long Z = 1; for (int j = 0; j < K; j++) { Z *= pow(p[j], q[j]); }
しかし, この場合 pow(p[j], q[j]) の部分が $O(1)$ でなければ \(\displaystyle O\left(\frac{N\log M}{\log\log M}\right)\) は達成できません. そこで, $\displaystyle K=O\left(\frac{\log M}{\log\log M}\right)$ なので全て前計算してしまいましょう. $Z$ は $N!$ を割り切るので, 各 $p_i$ に対して $\displaystyle \sum_{j=1}^{\infty}\left\lfloor\frac{N}{{p_i}^j}\right\rfloor\leq\frac{N}{p_i-1}$ まで計算しておけば十分です. そして, この前計算に掛かる計算量は $(4)$ より \(\displaystyle O\left(\frac{N\log M}{\log\log M}\right)\) です.
以上より, 全体 \(\displaystyle O\left(\frac{N\log M}{\log\log M}\right)\) で $C(N,i)$ を列挙することが出来ました.
実装↓ (篩パートで計算量に $\log\log N$ を付けたくなかったので, 線形篩を用いています)
*1:$2\cdot 3\cdot 5\cdot 7\cdot 11\cdot 13\cdot 17\cdot 19\cdot 23\cdot 29=6469693230$ で, \(6469693230^2\geq 2^{63}\) です.
*2:この辺りの分かりやすい/一般的な記号があれば教えて下さい...
*3:\(\max\) を付けているのは, 尺取りをするためです
*4: 参考 : 線形篩 – 37zigenのHP
