階乗 mod 素数
Many Factorials を準備したので宣伝を兼ねて。
$N! \bmod P$ ($P$ は固定された素数) をたくさんの非負整数 $N$ に対して評価する状況を考える。
$N$ の上限が小さい (例えば $N \leq 10 ^ 7$ が保証される) 場合
$N$ の上限を $M $ とする。
予め $N! \bmod P$ を全ての $N=0,1,\ldots,M $ に対して前計算しておくことで、各クエリを定数時間で処理できる。
$N$ の上限が大きい場合
本記事の本題。
方法① 埋め込み
$N \geq P$ ならば $N! \equiv 0 \bmod P$ なので、$N\lt P$ を仮定する。
適当なブロックサイズ $B$ を決めて $(iB)! \bmod P\ (i=0,1,\ldots,\lfloor (P-1) / B\rfloor$ を時間を掛けて計算し、その結果をソースコードに直接埋め込む。各クエリは前計算の結果として得られている $(\lfloor N/B\rfloor B)!$ からの差分だけを計算することで $O(B)$ 時間で処理できる。
なお、定数倍改善として、$\lfloor N/B\rfloor B$ と $(\lfloor N/B\rfloor + 1) B$ のうちより $N$ に近い方からの差分を計算するという方法が考えられる。$(\lfloor N/B\rfloor + 1) B \geq P$ となるケースに注意する必要があるが、この場合は $(P - 1)! \equiv -1 \pmod{P}$ (ウィルソンの定理 - Wikipedia) を用いて $(P - 1)!$ からの差分を計算するとよい。
この方法は手軽だが、次の点に注意する必要がある。
- $P$ が実行時に与えられる場合は使えない
- ソースコード長の制限により $B$ を小さく取れないことがある
方法② $(iB)! \bmod P$ の実行時計算
本節の内容は(現在は削除されている)min_25さんの記事で紹介されていた手法を元にしている。
$(iB)! \bmod P$ を実行時に計算することで方法①の弱点を克服する。
以下では Shift of Sampling Points of Polynomial (以下評価点シフトと呼ぶ) を道具として認める。詳細については Shift of Sampling Points of Polynomial | cp-library-cpp などを参照。
簡単のため $B = 2 ^ K$ とする。また $\displaystyle f _ i(x) \coloneqq \prod _ {j = 1} ^ {2 ^ i - 1} (2 ^ i x + j) = (2 ^ i x + 1) (2 ^ i x + 2) \cdots (2 ^ i x + 2 ^ i - 1)$ と定める。
$f _ K(0), f _ K(1), \ldots, f _ K(\lfloor P/2 ^ K\rfloor)$ を計算することが目標である。
さて、$f _ i$ と $f _ {i + 1}$ の間には次の関係が成り立つ。
$$\begin{aligned} f _ {i + 1} (x) &= \prod _ {j = 1} ^ {2 ^ {i + 1} - 1} (2 ^ {i + 1} x + j) \\ &= \left(\prod _ {j = 1} ^ {2 ^ i - 1} (2 ^ {i + 1} x + j)\right)\cdot (2 ^ {i + 1} x + 2 ^ i) \cdot \left(\prod _ {j = 1} ^ {2 ^ i - 1} (2 ^ {i + 1} x + 2 ^ i + j)\right) \\ &= \left(\prod _ {j = 1} ^ {2 ^ i - 1} (2 ^ i (2x) + j)\right)\cdot 2 ^ i (2x + 1) \cdot \left(\prod _ {j = 1} ^ {2 ^ i - 1} (2 ^ i (2x + 1) + j)\right) \\ &= f _ i(2x) \cdot f _ i(2x + 1) \cdot 2 ^ i (2x + 1). \end{aligned}$$
つまり $f _ i(0), \ldots, f _ i(2 ^ i - 1)$ が既知であると仮定すれば、次の手続きにより $f _ {i + 1}(0), \ldots, f _ {i + 1}(2 ^ {i + 1} - 1)$ を計算することができる。
- 評価点シフトを用いて $f _ i(2 ^ i), f _ i(2 ^ i + 1), \ldots, f _ i(2\cdot 2 ^ {i + 1} - 1)$ を計算する。
- 評価点シフトを用いるための前提条件として、$f _ i(0), f _ i(1), \ldots, f _ i(\deg f _ i)$ が既知でなければならないことに注意する。今回は $\deg f _ i = 2 ^ i - 1$ より条件を満たしていることを確認できる。
- 各 $j = 0, 1, \ldots, 2 ^ {i + 1} - 1$ に対する $f _ {i + 1}(j)$ を、$f _ {i + 1}(j) = f _ i(2j) \cdot f _ i(2j + 1) \cdot 2 ^ i (2j + 1)$ により得る。
1 は $O(i 2 ^ i)$ 時間で、2 は $O(2 ^ i)$ 時間で可能である。結局 $f _ 0(0) = 1$ から初めて上記の手続きを繰り返すことで $O(\sum _ {i = 0} ^ {K - 1} i 2 ^ i) = O(K 2 ^ K)$ 時間で $f _ K(0), \ldots, f _ K(2 ^ K - 1)$ を得ることができる。
$2 ^ K \leq \lfloor P/2 ^ K\rfloor$ の場合、最後に $f _ K(0), \ldots, f _ K(2 ^ K - 1)$ に対してもう一度評価点シフトを用いることで $f _ K(2 ^ K), f _ K(2 ^ K + 1), \ldots, f _ K(\lfloor P/2 ^ K\rfloor)$ を得ることができる。これは $O((P / 2 ^ K) \log P)$ 時間で可能である。
結局、$2 ^ K$ を $\sqrt{P}$ 付近に取れば、前計算は $O(\sqrt{P} \log P)$ 時間で行うことができる。
実用においては、$K$ をより小さく取って前計算に時間を掛け、$1$ 回のクエリ処理に掛かる時間を短くするとよい (実験結果 を参照)。
方法③ $(iB)! \bmod P$ の実行時計算 + オフライン処理による高速化 (1)
クエリをオフラインで処理してよいと仮定する。
$(iB)! \bmod P$ を前計算しておけば、各クエリは $\displaystyle (\lfloor N/B\rfloor B)! \cdot N(N-1)\cdots(N-(N\bmod B)+1)$ として計算できる。$(\lfloor N/B\rfloor B)!$ は前計算の結果から既知なので、各 $N$ に対する $N(N-1)\cdots(N-(N\bmod B)+1)$ を高速に計算したい。
つまり、$\displaystyle f(a,b)\coloneqq \prod _ {i = 0} ^ {b - 1} (a - i)$ と定義すれば、$f(N, N \bmod B)$ が計算したいものである。
各クエリで $N \bmod B$ が全て等しい場合は比較的容易であり、その値を $R$ として、$\displaystyle \prod _ {i = 0} ^ {R - 1} (x - i)$ を $R$ 点ずつに分けて多点評価すると $O( (R + Q) (\log R) ^ 2)$ 時間となる ($Q $ はクエリの数)。
一般にはこのような性質は成り立たない。しかし、$0\leq i\leq b$ なる整数 $i$ に対して $f(a, b) = f(a, i) \cdot f(a - i, b - i)$ が成り立つことを用いて $N \bmod B$ を $2$ 冪の和に分解することで、$b$ の種類数を高々 $\lfloor \log _ 2 (B-1) \rfloor$ に減らすことができる。
具体的には、例えば $N = 49,\ N \bmod B = 13 = 2 ^ 0 + 2 ^ 2 + 2 ^ 3$ のとき $f(49,13) = f(49, 2 ^ 0) \cdot f(49 - 2 ^ 0, 2 ^ 2) \cdot f(49 - 2 ^ 0 - 2 ^ 2, 2 ^ 3)$ のように分解する。
つまり各 $t = 0, 1, \ldots, \lfloor \log _ 2 (B-1) \rfloor$ に対して $\displaystyle \prod _ {i = 0} ^ {2 ^ t - 1} (x - i)$ を高々 $Q$ 個の点で評価すればよく、計算量は $\displaystyle O\left( \sum _ {i = 0} ^ {\lfloor \log _ 2 (B-1)\rfloor - 1} (2 ^ i + Q) i ^ 2 \right) = O(Q (\log B) ^ 3 + B (\log B) ^ 2)$ となる。
$B$ を $\sqrt{P}$ 付近に取ると、前計算と併せて全体で $O(Q (\log P) ^ 3 + \sqrt{P} (\log P) ^ 2)$ 時間となる。
この方法でも、実用上は $B$ をより小さく取って多点評価の回数を減らした方が速い (実験結果 を参照)。
方法④ $(iB)! \bmod P$ の実行時計算 + オフライン処理による高速化 (2)
方法③では多点評価を $O(\log B)$ 回に分けて行ったが、まとめて $O((Q + \sqrt{P}) (\log P) ^ 2)$ 時間で行うことができる (参考: yukicoder Paint and Fill 解法詳細解説 | Mathenachia)。
この方法では $B$ の大きさによる実行時間への影響は小さい (実験結果 を参照)。
実験
$P = 998244353,\ Q = 5\times 10 ^ 5$ のランダムケースで $K = \log _ 2 B$ を動かして実行時間 (秒) を計測した。なお、③' は除算を用いない高速な多点評価を用いて③を高速化したものである。コードの詳細については factorial mod prime · GitHub を参照。
| 方法＼$K$ | $7$ | $8$ | $9$ | $10$ | $11$ | $12$ | $13$ | $14$ | $15$ |
|---|---|---|---|---|---|---|---|---|---|
| ② | 0.80 | 0.53 | 0.46 | 0.56 | 0.90 | 1.6 | 3.0 | 6.0 | 12 |
| ③ | 1.0 | 1.1 | 1.3 | 1.7 | 2.2 | 2.7 | 3.4 | 4.1 | 4.9 |
| ③' | 0.90 | 0.78 | 0.83 | 1.0 | 1.3 | 1.6 | 2.0 | 2.5 | 3.4 |
| ④ | 1.6 | 1.4 | 1.3 | 1.3 | 1.3 | 1.4 | 1.6 | 1.7 | 2.1 |
$P$ が小さいため方法②が一番速いという結果になった。また、③や④はクエリをオンラインで処理できない一方で②はオンラインで処理ができるため、競技の文脈では (大抵の場合 $P \lt 2 ^ {31}$ なので) 方法②を使用するのが無難だと思われる。
実装
最も高速であった方法②で $K=9$ としたものの実装を以下に示す。
