$\binom{n}{m} \bmod p^e$ を高速に解く。（$n, m, p^e \lt 10^{300}$, $pe \lt 10^6$, $p$: 素数）
$\newcommand{\stirlingfirst}[2]{\genfrac{[}{]}{0pt}{}{#1}{#2}}$Andrew Granville さんの Binomial coefficients modulo prime powers (BinCoeff.pdf) とは違った計算方法。
前計算後の $\binom{n}{m} \bmod{p^e}$ の計算量（乗算・剰余算の回数）を、$O\left(e^2 \left(\frac{\log{n}}{\log{p}} + \min(\log{n}, p \log {p}) \right) \right)$ (?) から $O(e \log{n})$ にできる。
（$\genfrac{[}{]}{0pt}{}{n}{m}$: 第 1 種 Stirling 数）
計算方法
$(n!)_p$ を $i \not\equiv 0 \pmod{p}$ を満たす $i$ ($1 \le i \le n$) の積と定義するとき、$(n!)_p \bmod p^e$ が求められれば十分。
まず $u, v$ ($0 \le v < p$) を用いて $n = up + v$ とかく。このとき、$$ \begin{aligned} ((up + v)!)_p &= \left(\prod_{i=0}^{u-1} \prod_{j=1}^{p-1} (ip + j)\right) \cdot \prod_{j=1}^{v} (up+j) \\ &\equiv \left(\prod_{i=0}^{u-1} \left( \sum_{k=0}^{e-1} (ip)^k \genfrac{[}{]}{0pt}{}{p}{k+1} \right) \right) \cdot \left(\sum_{k=0}^{e-1} (up)^k \genfrac{[}{]}{0pt}{}{v+1}{k+1}\right) \pmod{p^e} \\ &= \genfrac{[}{]}{0pt}{}{p}{1}^u \left(\prod_{i=0}^{u-1} \left(1 + \sum_{k=1}^{e-1} \frac{\genfrac{[}{]}{0pt}{}{p}{k+1}}{\genfrac{[}{]}{0pt}{}{p}{1}} (ip)^k \right) \right) \cdot \left(\sum_{k=0}^{e-1} (up)^k \genfrac{[}{]}{0pt}{}{v+1}{k+1}\right) \\ &\equiv \genfrac{[}{]}{0pt}{}{p}{1}^u f_{p,e}(u) \left( \sum_{k=0}^{e-1} (up)^k \genfrac{[}{]}{0pt}{}{v+1}{k+1} \right) \pmod{p^e} \end{aligned} $$ が成り立つ。ここで、$f_{p,e}$ は適当な $(2e-2)$ 次以下の多項式である（後で証明）。よって、$f_{p,e}(0), \ldots, f_{p,e}(2e-2)$ を計算しておくことで $O(e \log{p})$ で $f_{p,e}(u)$ の値を補間できる。
したがって、適当な前計算を行っておけば、$(n!)_p \bmod p^e$ は $O(e \log{p})$ で計算できる。
実装方法と計算量
1. $\genfrac{[}{]}{0pt}{}{n}{m}$ ($1 \le n \le p$, $1 \le m \le \min(p, e)$) の前計算: $O(p \cdot \min(p, e))$
2. $f_{p,e}(0), ..., f_{p,e}(2e-2)$ や多項式補間用の前計算: $O(e \cdot \min(p, e) + e \log{p})$
3. 必要な $(n_i!)_p$ の計算: $O(\log_p{n} \cdot (e \log{p})) = O(e \log{n})$
合計で $O(pe + e \log n)$ 回の乗算・剰余算で、$\binom{n}{m} \bmod p^e$ を計算できる。
多項式となる証明
$a_k := \frac{\genfrac{[}{]}{0pt}{}{p}{k+1}}{\genfrac{[}{]}{0pt}{}{p}{1}} \mod p^e$ とし、$$g(u) := \prod_{i=0}^{u-1} \left(1 + \sum_{k=1}^{e-1} a_k (ip)^k\right)$$ とおくと、$$\begin{aligned} s_0(u) &= 1, \\ s_k(u) &= \sum_{j=1}^k \left( \sum_{i=0}^{u-1} s_{k-j}(i) \cdot i^j\right) a_j \,\,\,(k \ge 1) \end{aligned}$$ なる $2k$ 次の多項式 $s_k$ ($0 \le k < e$) を用いて、
$$g(u) \equiv \sum_{k=0}^{e-1} s_k(u) p^k \pmod{p^e}$$ とかける。
その他
Andrew さんの方法は多倍長整数なしで実装しようとすると手間がかかるのだけど、この方法だと多倍長整数を持たなくても実装しやすい。
- $p$ が素数のとき、$\genfrac{[}{]}{0pt}{}{p}{i}$ ($2 \le i < p$) は $p$ で割り切れる
- $p \ge 5$ のとき $\genfrac{[}{]}{0pt}{}{p}{2}$ は $p^2$ で割り切れる
などといった性質を用いれば、もう少し $f_{p,e}(x)$ の次数を落とせるのだと思う。
（実験的には、$e - 1 + \lfloor(p-1)/(e-1)\rfloor$ 次まで落とせそう(?)）
