扩展 Lucas 定理
\(\newcommand{\stirlingI}[2]{\genfrac{[}{]}{0pt}{}{#1}{#2}}\)Lucas 定理想必大家都知道了：令 \(p\) 为一个质数，\(n = u_1 p + u_2\)，\(m = v_2 p + v_2\)，其中 \(0 \leq u_2, v_2 < p\)，则有
\[
\binom{n}{m} \equiv \binom{u_1}{v_1} \binom{u_2}{v_2} \pmod{p}.
\]
但是 Lucas 定理只对模为质数的情况有效。虽然对于合数我们有 CRT，但是还是得处理模为质数幂的情况。这篇文章主要讲模为质数幂该如何处理。文章非常多的地方参考翻译了 prabowo 的博客（虽然他也是参考了 min_25 的博客），但是也有一些小修改（虽然他也有一些小修改），主要就是改进了一下 bound。
之前 Andrew Granville 也有一篇 关于拓展 Lucas 定理的文章，但是里面特判太多了……复杂度也高，我写过一次之后不想再写了……
拆开 \(n!\)
令 \(E(n) = \max_k \{p^k | n!\}\) 表示 \(p\) 进制表示中 \(n!\) 末位 0 的个数，显然 \(E(n)\) 非常好求。我们在定义 \[ (n!)_p = \prod_{i=1, p \nmid i}^n i, \] 即 \([n]\) 中所有非 p 的倍数的乘积，于是 \(n!\) 可以表示为 \[ n! = p^{E(n)} \prod_{k=0} \left( \left\lfloor n/p^k \right\rfloor !\right)_p, \] 就可以做乘除法了。我们的任务即为快速计算 \((n!)_p \bmod p^e\) 。
不妨令 \(n = up + v\)，我们可以将 \((n!)_p\) 拆成两部分： \[ (n!)_p = \left( \prod_{i=0}^{u-1} \prod_{j=1}^{p-1} (ip + j) \right) \left( \prod_{i=1}^v (up + j) \right), \] 注意到两者均有 \(\prod_{j=1}^r (ip + j)\) 这么一个结构，可以将其视作关于 \(p\) 的多项式并展开，而我们只需要这个乘积 \(\bmod p^e\) 之后的结果。具体说来，我们有 \[ \prod_{i=0}^{k-1} (n + i) = \sum_{i=0}^k \stirlingI{k}{i} n^i, \] 其中 \(\stirlingI{k}{i}\) 为 第一类 Stirling 数，递推为 \[ \stirlingI{n + 1}{k} = n \stirlingI{n}{k} + \stirlingI{n}{k-1}, \quad \stirlingI{0}{0} = 1, \stirlingI{0}{n} = \stirlingI{n}{0} = 0. \] 这样 \((n!)_p\) 即可化为： \[ \begin{aligned} (n!)_p & = \left( \prod_{i=0}^{u-1} \prod_{j=1}^{p-1} (ip + j) \right) \left( \prod_{i=1}^v (up + j) \right) \\ & = \left( \prod_{i=0}^{u-1} \sum_{k=1}^{p} \stirlingI{p}{k} (ip)^{k-1} \right) \left( \sum_{k=1}^{v} \stirlingI{v+1}{k} (up)^{k-1} \right) \\ & \equiv \left( \prod_{i=0}^{u-1} \sum_{k=1}^{e} \stirlingI{p}{k} (ip)^{k-1} \right) \left( \sum_{k=1}^{e} \stirlingI{v+1}{k} (up)^{k-1} \right) \pmod{p^e} \\ & = \stirlingI{p}{1}^u \left( \prod_{i=0}^{u-1} \sum_{k=1}^{e} \frac{\stirlingI{p}{k}}{\stirlingI{p}{1}} (ip)^{k-1} \right) \left( \sum_{k=1}^{e} \stirlingI{v+1}{k} (up)^{k-1} \right) \end{aligned} \]
我们令 \[ f_{p, e}(u) := \prod_{i=0}^{u-1} \sum_{k=1}^{e} \frac{\stirlingI{p}{k}}{\stirlingI{p}{1}} (ip)^{k-1} \bmod p^e = \prod_{i=0}^{u-1} \left( 1 + \sum_{k=1}^{e-1} \frac{\stirlingI{p}{k+1}}{\stirlingI{p}{1}} (ip)^k \right) \bmod{p^e}, \]
现在 \((n!)_p\) 就被我们分成了三部分：\(\stirlingI{p}{1}^u\)，\(f_{p, e}(u)\) 和 \(\sum_{k=1}^{e} \stirlingI{v+1}{k} (up)^{k-1}\)，第一项和第三项都很好算，压力就传到 \(f_{p, e}(u)\) 上了。
现在给出一个大结论：\(f_{p, e}(u)\) 是一个关于 \(u\) 的多项式。
证明 \(f_{p, e}\) 为关于 \(u\) 的多项式
我们证明一个更一般的结论。对于 \(1 \leq k < e\)，令 \(w_k(i)\) 为一个关于 \(i\) 的多项式，再令 \[ g(x, u) = \prod_{i=0}^{u-1} \left( 1 + \sum_{k=1}^{e-1} w_k(u) x^k \right) \mod x^e, \] 则有 \(g\) 是一个关于 \(u\) 的多项式，其次数不超过 \[ \deg_u g(x, u) \leq \lambda (e-1), \quad \lambda = \max_{1 \leq k < e} \frac{\deg w_k + 1}{k}. \]
证明如下：令 \(0 \leq k < e\)，\(h_k(u)\) 为 \(g(x, u)\) 中 \(x^k\) 前的系数，我们对 \(k\) 做归纳，证明 \(h_k(u)\) 是一个关于 \(u\) 的多项式，其次数不超过 \(\lambda k\)。
- 若 \(k = 0\)，则有 \(h_0(u) = 1\)，显然满足要求；
- 若 \(0 < k < e\)，则 \[ h_k(u) = h_{k-1}(u-1) + \sum_{t=0}^{k-1} h_t(u-1) w_{k-t}(u - 1) = h_{k-1}(0) + \sum_{t=0}^{k-1} \sum_{u' \geq 0}^{u-1} h_t(u') w_{k-t}(u'), \] 容易看出 \(h_t(u') w_{k-t}(u')\) 为一个关于 \(u'\) 的 \(\deg h_t + \deg w_{k-t}\) 次多项式，故有 \(h_k\) 也为一个关于 \(u\) 的多项式，其次数不超过 \[ \deg h_k \leq \max_{0 \leq t < k} \left( \deg h_t + \deg w_{k-t} + 1 \right), \] 这个 1 就是因为有个对于 \(u'\) 的求和。由 \(\lambda\) 定义知 \[ \deg w_{k - t} + 1 \leq \lambda (k - t), \] 再加上对 \(k\) 做的归纳，有： \[ \deg h_k \leq \max_{0 \leq t < k} \left( \lambda t + \lambda (k - t) \right) = \lambda k, \] 归纳成立。熟悉的朋友很快就可以看出来，这个 \(\deg h_k\) 的递归不就是做一个 \(\deg w_k + 1\) 的背包嘛，上界自然就是“密度”最大的多项式了。
剩下的就很简单了，我们有 \[ g(x, u) = \sum_{k=0}^{e-1} h_k(u) x^k \quad \Longrightarrow \quad \deg_u g \leq \max_{0 \leq k < e} \deg h_k \leq \lambda (e-1). \]
\(f_{p, e}\) 的次数上界
现在还剩下一个问题：既然 \(f_{p, e}\) 是个多项式，那么次数最高多少？
如果我们按照上文的定义来， \[ f_{p, e}(u) = \prod_{i=0}^{u-1} \left( 1 + \sum_{k=1}^{e-1}\frac{\stirlingI{p}{k+1}}{\stirlingI{p}{1}} (ip)^k \right) \mod p^e, \] 故对于所有 \(1 \leq k < e\)，令 \(w_k(i) = \frac{\stirlingI{p}{k+1}}{\stirlingI{p}{1}} i^k\) 即可，此时有 \[ \lambda = \max_{1 \leq k < e} \frac{\deg w_k + 1}{k} = \max_{1 \leq k < e} \frac{k + 1}{k} = 2, \quad \Longrightarrow \quad \deg_u f_{p, e} \leq 2 (e - 1). \]
但是聪明一点看了 prabowo 的提示我们可以发现：当 \(p\) 为素数且 \(1 < k < p\) 时，\(p\) 能整除 \(\stirlingI{p}{k}\)。所以如果 \(e < p\)，那涉及到所有的 \(\stirlingI{p}{k+1}\) 均为 \(p\) 的倍数，我们就可以把这个 \(p\) 的因子提出来：
\[
\begin{aligned}
f_{p, e} & = \prod_{i=0}^{u-1} \left( 1 + \sum_{k=1}^{e-1}\frac{\stirlingI{p}{k+1}}{\stirlingI{p}{1}} (ip)^k \right) \mod p^e \\
& = \prod_{i=0}^{u-1} \left( 1 + \sum_{k=1}^{e-1}\frac{\stirlingI{p}{k+1} / p}{\stirlingI{p}{1}} i^k p^{k+1} \right) \mod p^e.
\end{aligned}
\]
这样我们就可以对于 \(k > 1\)，令 \(w_k(i) = \frac{\stirlingI{p}{k} / p}{\stirlingI{p}{1}} i^{k-1}\) 以及 \(w_1(i) = 0\)，此时的 \(\lambda\) 就可以更小：
\[
\lambda = \max_{1 \leq k < e} \frac{\deg w_k + 1}{k} = 1, \quad \Longrightarrow \quad \deg_u f_{p, e} \leq e - 1.
\]
这里其实不太严谨，因为 \(\bmod p^e\) 下不能进行 \(/p\)，但是我们可以现在实数范围内做。这里的除法只是为了证明，代码里面不会用到。
另外一种情况是 \(e \geq p\)。这里就多涉及到了一个数 \(\stirlingI{p}{p} = 1\)，所以就只能整理得到 \[ \begin{aligned} f_{p, e} & = \prod_{i=0}^{u-1} \left( 1 + \sum_{k=1}^{e-1}\frac{\stirlingI{p}{k+1}}{\stirlingI{p}{1}} (ip)^k \right) \mod p^e \\ & = \prod_{i=0}^{u-1} \left( 1 + \sum_{k=1}^{p-2}\frac{\stirlingI{p}{k + 1} / p}{\stirlingI{p}{1}} i^k p^{k+1} + \frac{1}{\stirlingI{p}{1}} (ip)^{p-1} \right) \mod p^e. \end{aligned} \] 故有对于 \(1 < k < p\)，\(w_k(i) = \frac{\stirlingI{p}{k} / p}{\stirlingI{p}{1}} i^k\)，以及 \(w_{p-1} = \frac{\stirlingI{p}{p-1}/p}{\stirlingI{p}{1}} {i^{p-2}}+ \frac{1}{\stirlingI{p}{1}} i^{p-1}\) 。此时的 \(\lambda\) 为： \[ \lambda = \max_{1 \leq k < e} \frac{\deg w_k + 1}{k} = \frac{p}{p - 1}, \quad \Longrightarrow \quad \deg_u f_{p, e}(u) \leq \frac{p}{p - 1} (e - 1). \] 由于 \(p \geq 2\)，所以这个 bound 总是比 \(2(e-1)\) 好的。
综上所述，我们有 \[ \deg_u f_{p, e}(u) \leq d^* := \left\lfloor \frac{p (e - 1)}{p-1} \right \rfloor. \]
Stirling 数整除性
上文用到了一个性质：当 \(p\) 为素数且 \(1 < k < p\) 时，\(p\) 能整除 \(\stirlingI{p}{k}\)。这里来稍微证一下。
我们考虑 \(\mathbb{F}_p\) 上的多项式 \(f(x) = \prod_{i=1}^{p-1} (x+i) \in \mathbb{F}_p[x]\)，由 Stirling 数性质可得： \[f(x) = \prod_{i=1}^{p-1} (x + i) = \sum_{i=1}^p \stirlingI{p}{i} x^{i-1},\] 注意到 \(f(0) = (p-1)! \equiv -1 \pmod{p}\)，对于 \(x \neq 0\)，则有 \(f(x) = 0\)；另一方面，多项式 \(g(x) = x^{p-1} - 1\) 在这 \(p\) 个点上的值和 \(f\) 一样，而 \(f\) 和 \(g\) 次数均为 \(p - 1\) 次。由多项式的唯一性可得 \(f = g\)，对比系数即有： \[ \stirlingI{p}{i} \equiv \begin{cases} 1, & i = p, \\ 0, & 1 < i < p, \\ -1, & i = 1.\end{cases} \pmod{p} \]
算法
- 对于 \(0 \leq v \leq p, 0 \leq k \leq e\)，预处理出 \(\stirlingI{v}{k}\)，复杂度为 \(O(pe)\)；
- 预处理出 \(f_{p, e}\) 。由于 \(d^*\) 次多项式需要 \(d^* + 1\) 个点来确定，暴力点的做法就是对于 \(0 \leq u \leq d^*\)，直接暴力求出 \(f_{p, e}(u)\)，然后 Lagrange interpolate 一下，复杂度为 \(O(e^2)\)。
- 询问就直接套公式： \[ (n!)_p \equiv \stirlingI{p}{1}^u f_{p, e}(u) \left( \sum_{k=1}^{e} \stirlingI{v+1}{k} (up)^{k-1} \right) \pmod{p^e} \] 单次询问复杂度为 \(O(e)\)。
如果 \(d^* \geq p\)，那么 Lagrange interpolation 需要小心一点：我们直接有公式 \[ f_{p, e}(u) = \sum_{i=0}^{d^*} f_{p, e}(i) \prod_{j \neq i} \frac{u - j}{i - j} = \sum_{i=0}^{d^*} f_{p, e}(i) \binom{u}{j}\binom{u-j-1}{k-j} (-1)^{k-j}, \] 所以得到的的一定是整数。但是 \(\prod_{j \neq i} \frac{u - j}{i - j}\) 这个式子不是很好求，原因就是 \(i - j\) 在 \(\bmod{p^e}\) 下不一定有逆元……所以我们就只能把所有数写成 \(r \equiv \frac{n}{p^a} \pmod {p^e}\) 这种形式后再做乘除法了……
代码有 prabowo 的 C++ 版本，我写了一个 Python 版本。之前还有一个 min_25 的版本，但是可惜 min_25 博客现在已经无了，代码也已经删库跑路了，我记得里面有好多有趣的东西呢……
Hexo
一开始写的时候因为用了 \newcommand{\aaa}[1]{#1} 这样的命令，而 Hexo 把这个 # 认成了 SwigTags 了，所以渲染一直有问题。解决方法很简单，在 front-matter 里面加一句
1 | disableNunjucks: true |
就可以了。参考 这个 issue （注意这个 comment 里面的 disableNunjucks 打错字了）和 这个 issue。