この記事は？

二項係数 mod 合成数を列挙します。素因数ごとの結果を並列処理することで高速化する方法を紹介します。
並列処理をする部分以外は suisenさんの記事と同様です。 suisen さんの記事とかぶるところは本記事では省略している部分もあるので、詰まったら参照してください。

二項係数 mod 合成数

やりたいこと

正整数 $N$ および $M$ が与えられたとき、 $_{N} C_{i} \ (i=0,1,\cdots,N)$ を $\bmod M$ で列挙する

競技プログラミングで現実的に問われる範囲として、 $M$ は $10^9$ 程度以下（より正確には $M^2 \le 2^{63}-1$ ）、 $N$ は $10^5$ から $10^7$ 程度を想定しています。以下、「上述の範囲」と言うとこの範囲を指します。

概要

suisen さんの方法

$M$ の（$N$ 以下の）素因数の個数を $L$ とする。
$M$ と互いに素な部分は$O(N)$ で計算できる。
$M$ と共通素因数を持つ部分については、各素因数 $p$ について、 $_{N} C_{i}$ が $p$ で何回割れるかを調べることで（素因数ごとに） $O(N)$ でできる
全体計算量は $O(NL)$ ¹
$L$ は上述の範囲で $9$ 以下

本記事の方法

$L$ 個の素因数を $K$ 個のグループに分けることで、共通素因数を持つ部分を高速化する。
$K=\lceil \displaystyle\frac{L}{\lfloor \log_{\lfloor\log_2 N\rfloor + 1} N \rfloor} \rceil$ とすると、全体計算量は $O(N(K+\log\log L))$
$K$ は上述の範囲では（前計算がボトルネックにならないことを満たしつつ） $K\le 2$ が達成可能。

指数の上限について

各素数 $p$ について、 $p$ が $_{N} C_{i}$ を割り切る回数の最大値はどれぐらいでしょうか。 $p$ が $x$ を割り切る回数を $f(x)$ と書くと

$f( _{N}C_{i})$
$=f(\displaystyle\frac{N!}{i!(N-i)!})$
$=f(N!)-f(i!)-f((N-i)!)$
$=\displaystyle\sum_{e=1}^{\lfloor \log_p N\rfloor} \bigg( \bigg\lfloor \frac{N}{p^e}\bigg\rfloor
-\bigg\lfloor \frac{i}{p^e}\bigg\rfloor
-\bigg\lfloor \frac{N-i}{p^e}\bigg\rfloor \bigg)$
$\le \lfloor \log_p N\rfloor$

であることから、 $\le \lfloor \log_p N\rfloor$ 回以下であることが分かります（最後の不等号は、カッコ内が $1$ 以下であることから分かります ² ）。

p 進法で a-b と b を足した時の繰り上がりの回数と等しくなります。log_p(a) で抑えられます。
— 熨斗袋 (@noshi91) April 2, 2023

計算ステップ

ここから計算方法の詳細を説明します。

前計算

$M$ の素因数のうち $N$ 以下のもの（ $p_1,\cdots, p_L$ とする）をすべて列挙する
各 $p_i$ について、 $p_i^e\ (e=0,1,\cdots, \lfloor \log_2 N \rfloor)$ を $\bmod M$ で求める ³
$P=\{p_1,\cdots,p_L\}$ の分割 $S = \{S_k\} _{k=1}^{K}$ ⁴ を良い感じに決める
各 $S_k = {\{p_{k,1},\cdots,p_{k,l_k}\}}$ について、 $E_k=\displaystyle\max_{p\in S_k}(\lfloor\log_p N\rfloor+1)$ とする。各 $\{e_{k,i}\}_{i=1}^{l_k}\ (0\le e_{k,i} \lt E_k)$ について、 $\displaystyle\prod_{i=1}^{l_k}{p_{k,i}}^{e_{k,i}}$ を $\bmod M$ で求め、指数の辞書順にテーブルに入れておく（つまり大きさ ${E_k}^{l_k}$ のテーブルになる） ⁵ 。

テーブルの構築

$_{N} C_{i}$ は $\displaystyle\frac{N-i}{i+1}\ (i=0,\cdots)$ の累積積なので各素因数 $p$ について、分子・分母が $p$ で割り切れる部分にメモを付けていけばよいです。ただし $N\times L$ のテーブルを作りたくないので、各グループごとにまとめて処理を行います。具体的には次のようにします。

答えを入れる長さ $N$ のテーブル $ans$ を作り、最初すべての要素を $1$ で初期化する。
各 $1\le k\le K$ について次を行う。
1. 長さ $N$ のテーブル $T$ を作り、最初すべての要素を $0$ で初期化する。
2. $S_k = {\{p_{k,1},\cdots,p_{k,l_k}\}}$ に含まれる素因数 $p_{k,i}$ について、分子が $p_{k,i}$ で割れる位置の $T[i]$ に ${E_k}^{i-1}$ を足し、分母が $p_{k,i}$ で割れる位置の $T[i]$ から ${E_k}^{i-1}$ を引く。
3. $T$ を累積和にする
4. 前計算で作ったテーブルにより、このグループの結果が分かるので、 $ans$ を更新する
可逆部分の計算を行う

計算量

オーダー評価

素因数の列挙は $O(N)$ ⁶、指数の前計算は $O\big(\displaystyle\sum_{k=1}^{K}{E_k}^{l_k}\big)$ 、各素因数で割れる位置をマークする部分が $O\big(\displaystyle\sum_{p\in P}\frac{N}{p-1}\big)\subset O(N\log\log L)$ 、グループごとの累積和の計算と $ans$ の更新が $O(NK)$ 、可逆部分が $O(N)$ です。このうち指数の前計算と $O(NK)$ 部分がトレードオフになっています（分割を粗くすると $O(NK)$ 部分が速くなるが、前計算の時間が増える）。

前計算がボトルネックにならないように、グループごとに作る前計算テーブルの大きさを $N$ 以下にすることを考えます。各グループの要素数 $l$ を $\lfloor \log_{\lfloor\log_2 N\rfloor + 1} N \rfloor$ 以下にすれば十分です。これは $K=\lceil \displaystyle\frac{L}{\lfloor \log_{\lfloor\log_2 N\rfloor + 1} N \rfloor} \rceil$ でできます。これにより全体の時間計算量

$$O\bigg(N\Big(K+\log\log L\Big)\bigg) \subset O\bigg(N\Big(\displaystyle\frac{L\log\log N}{\log N}+\log\log L\Big)\bigg) \subset O\bigg(N\Big(\displaystyle\frac{\log M \log \log N}{\log \log M\log N}+\log\log \log M\Big)\bigg)$$

を達成できます。空間計算量はグループごとのテーブルを使いまわせば $O(N)$ です。

最悪ケースの評価と分割方法

上述の範囲での最悪ケースとして $P=\{2,3,5,7,11,13,17,19,23\}$ の場合 ⁷ を考えると、分割を例えば $S=\{\{2,3,5,7\},\{11,13,17,19,23\}\}$ とすると $K=2$ が達成できます ⁸ 。これ以外の場合でも、分割数が $2$ 以下で最も計算量の合計が小さくなるものを選べば良いです ⁹ ¹⁰ 。

問題・実装例

問題
 AC コード（PyPy 3）

ちなみに PyPy だと $N=5\times 10^6$ ぐらいが（競技プログラミングで使うには）限界っぽいです。

できた。 PyPy だと 500 万で 1.6 秒ぐらい
— きり (@kiri8128) April 5, 2023

厳密には $L=0$ の可能性があるので $O(N(L+1))$ と書くべきだが、煩雑になるので省略する。後の $O(NK)$ も同様。 ↩
$p$ 進法で $i$ と $N-i$ を足したときに、下から $e$ 桁目で繰り上がりが発生すれば $1$ 、そうでなければ $0$ ↩
本質的に必要なのは $e\le \lfloor \log_{p_i} N \rfloor$ の範囲ですが、後の実装の都合で少し多めにしています。 ↩
$S = \{S_k\} _{k=1}^{K}$ が $P$ の分割であるとは、 $\displaystyle\bigcup_{k=1}^{K}S_k=P$ かつ $k\neq l$ について $S_k \cap S_l=\phi$ を満たすことをいう ↩
データ保持の方法として $E_k$ 進法のようにしていますが、各素因数ごとに必要な指数を考慮することでテーブルをより小さくする（ $\displaystyle\prod_{p\in S_k}(\lfloor\log_p N\rfloor+1)$ ）こともできます。ただしこの部分はボトルネックになりづらいので後のメインの処理を簡単にするためこのようにしています。 bit 演算で行うこともできます。 ↩
$N$ が大きいときは $O(\sqrt{M})$ やこちらの記事の方法でもできます。 ↩
$M=223092870=2\cdot 3\cdot 5\cdot 7\cdot 11\cdot 13\cdot 17\cdot 19\cdot 23$ のときなど ↩
このとき $10^5\le N\le 10^7$ の範囲で $\displaystyle\sum_{k=1}^{K}{E_k}^{l_k}\le N$ が達成できることが確認できます。 ↩
$L$ 通りを全部試すなどで十分です。 ↩
$N$ が $10^5$ より小さいオーダーだと分割数を増やした方が効率が良くなることもあります。 ↩