组合数求模(续)

0. 写在前面

本文是组合数求模的续篇，是对前一篇文章里第 $3$ 节内容的拓展。在本文中，一些前置讨论将简要带过，如果需要更细致的描述，请参考前文。

本文介绍的是利用 $\mathcal{O}(pe+e{\log }_{p}n)$ 次大数乘法计算 $((\genfrac{}{}{0ex}{}{n}{m})\mathrm{m}\mathrm{o}\mathrm{d}{p}^e)$ 的技巧，其中 $p$ 是质数，$n,m,e$ 是非负整数，并且默认 $p$ 和 $e$ 是小数字，$n$ 和 $m$ 是与 $p^e$ 同阶或比 $p^e$ 更高阶的大数字。

这里没有将大数运算的底层代价算入复杂度，主要是考虑到具体实现时存在优化，例如可以通过改用 $p$ 进制存储消除大数取模运算，也可以通过压缩存储进制降低常数。

有一定组合数学基础的读者可以先粗略阅读 Andrew Granville 关于组合数模质数幂次的研究综述， 但是笔者不建议细读这篇综述，因为它涉及太多其他方面的介绍，而且由于它的综述性质，很多定理都是直接引用（或略作改动）并不加证明的，这在一定程度上导致这篇综述里出现了很多明显的错误 。

本文可以看作是关于 Min_25 在前者的研究基础上进行的改进 的一篇翻译文章。（行文风格可能十分不同？）可惜笔者不才，只对 Min_25 的改进做出一些证明和常数优化，并仅对于 $p$ 较大的情况分析了其他可能的更优做法（详见 3.4 节）。

注：Min_25 似乎在计算操作数时漏掉了 $\frac{1}{\log p}$ 的因子，导致分析结果变大，实际算出的结果应该和本文一致。

author: skywalkert
original article: http://blog.csdn.net/skywalkert/article/details/104681101
last update time : 2020-03-14

1. 预备知识

本节介绍本文需要用到的新增知识及证明，只需要了解做法的读者可以先跳过此节 。（太长不看 ）如遇到未知的定义或定理，请参考前一篇文章或留言询问。

1.1 库默尔定理 (Kummer’s theorem)

令 ${\nu }_p(x)$ 表示最大的非负整数 $k$ 使得 $p^k\mid x$，其中 $p$ 是质数。库默尔定理表明：${\nu }_p\left(\left(\genfrac{}{}{0ex}{}{n}{m}\right)\right)$ 等于 $m$ 和 $(n-m)$ 在 $p$ 进制下做加法的进位次数。

以下是笔者的证明方法。

• 对于任意整数 $u\ge 0$, $0\le v<p$，显然有 ${\nu }_p((up+v)!)={\nu }_p(u!)+u$，故可以得到 ${\nu }_p(n!)={\sum }_{k\ge 1}\lfloor \frac{n}{p^k}\rfloor$。
• 若令 $n={\sum }_{i\ge 0}{n}_i{p}^i$ $(0\le n_i<p)$ 为 $n$ 在 $p$ 进制下的表示，且 $S_p(n)={\sum }_{i\ge 0}{n}_i$，那么有 ${\nu }_p(n!)={\sum }_{i\ge 0}{n}_i\frac{p^i-1}{p-1}=\frac{n-S_p(n)}{p-1}$。
• 对于组合数 $\left(\genfrac{}{}{0ex}{}{n}{m}\right)$，我们可以得到 ${\nu }_p\left(\left(\genfrac{}{}{0ex}{}{n}{m}\right)\right)={\nu }_p(n!)-{\nu }_p(m!)-{\nu }_p((n-m)!)=\frac{S_p(m)+S_p(n-m)-S_p(n)}{p-1}$。
• 令 $r=n-m$，考虑 $m$ 与 $r$ 在 $p$ 进制下进行加法的过程，令 $c_i$ 表示该过程中第 $i$ 位的结果向第 $(i+1)$ 位进位的次数，则有 $n_i=m_i+r_i+c_{i-1}-p{c}_i$，其中 $c_{-1}=0$。根据定义，我们可以得到 ${\nu }_p\left(\left(\genfrac{}{}{0ex}{}{n}{m}\right)\right)=\frac{S_p(m)+S_p(n-m)-S_p(n)}{p-1}={\sum }_{i\ge 0}{c}_i$，即库默尔定理。

顺带一提，对于多元组合数 $\left(\genfrac{}{}{0ex}{}{n}{m_1,m_2,\dots ,m_k}\right)$，库默尔定理依旧可以得到 ${\nu }_p\left(\left(\genfrac{}{}{0ex}{}{n}{m_1,m_2,\dots ,m_k}\right)\right)=\frac{1}{p-1}\left(-S_p(n)+{\sum }_{j=1}^k{S}_p(m_j)\right)={\sum }_{i\ge 0}{c}_i$，这里的 $c_i$ 由定义 $n_i=\left({\sum }_{j=1}^k{m}_j\right)+c_{j-1}-p{c}_j$ 得到。

1.2 上升幂 (Rising factorial) 与无符号第一类斯特林数 (Unsigned stirling numbers of the first kind)

定义关于 $x$ 的 $n$ 阶上升幂为 $x^{\overline{n}}={\prod }_{i=0}^{n-1}(x+i)$，显然这是一个 $n$ 次多项式，故可以定义 $x^{\overline{n}}={\sum }_{k=0}^n{\mathrm{c}\mathrm{o}\mathrm{e}\mathrm{f}\mathrm{f}}_{n,k}{x}^k$。

由 $x^{\overline{n+1}}=x^{\overline{n}}(x+n)$ 可以得到 ${\mathrm{c}\mathrm{o}\mathrm{e}\mathrm{f}\mathrm{f}}_{n+1,k}={\mathrm{c}\mathrm{o}\mathrm{e}\mathrm{f}\mathrm{f}}_{n,k-1}+n{\mathrm{c}\mathrm{o}\mathrm{e}\mathrm{f}\mathrm{f}}_{n,k}$，而初始情况为 ${\mathrm{c}\mathrm{o}\mathrm{e}\mathrm{f}\mathrm{f}}_{0,0}=1$, ${\mathrm{c}\mathrm{o}\mathrm{e}\mathrm{f}\mathrm{f}}_{0,k}=0$ $(k>0)$。

无符号第一类斯特林数 $\left[\genfrac{}{}{0ex}{}{n}{k}\right]$ 的定义有两种，一种就是 ${\mathrm{c}\mathrm{o}\mathrm{e}\mathrm{f}\mathrm{f}}_{n,k}$ 的定义，另外一种定义是大小为 $n$ 且能恰好划分为 $k$ 个环的置换数量。根据维基百科的相关词条，笔者认为第一种定义是最先出现的，也即 $x^{\overline{n}}={\sum }_{k=0}^n\left[\genfrac{}{}{0ex}{}{n}{k}\right]x^k$ 是最初的定义。

顺带一提，上升幂的系数还有一些特征，例如 $\left[\genfrac{}{}{0ex}{}{n}{1}\right]=(n-1)!$。

再顺带一提，根据上升幂和无符号的系数，读者可以相应地定义下降幂和有符号的系数，但此部分与后文无关，故这里略去。

1.3 威尔逊定理 (Wilson’s theorem) 及其扩展 (Gauss’s generalization)

威尔逊定理主要证明的是：当整数 $n>1$ 时，$n$ 是质数当且仅当 $(n-1)!\equiv -1(\mathrm{m}\mathrm{o}\mathrm{d}n)$。 虽然后文只会使用到其中的一个方向，但是出于科普的考虑，本小节将证明其充要性。

• 若 $n$ 能表示成两个不相等且大于 1 的整数之积，即 $n=xy$ $(1<x<y)$，则 $x,y<n$，于是 $xy|(n-1)!$，这意味着 $(n-1)!\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}n)$。
• 若 $n$ 是一个大于 2 的整数的平方，即 $n=x^2$ $(x>2)$，则 $2x<x^2=n$，于是 $2x^2|(n-1)!$，这意味着 $(n-1)!\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}n)$。
• 若 $n$ 是一个奇质数，则对于 $1\le u<n$ 且 $u^2\not\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}n)$ 的整数 $u$，总能找到唯一的整数 $v$ 满足 $1\le v<n$ 且 $uv\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}n)$，这意味着它们在 $((n-1)!\mathrm{m}\mathrm{o}\mathrm{d}n)$ 里能够一一抵消，另外由于模 $n$ 意义下存在原根（或者说模 $n$ 意义的缩系是一个循环群），则 $u^2\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}n)$ 最多只会有两个解，实际上这两个解恰好是 $1$ 和 $(n-1)$，故可以得到 $(n-1)!\equiv -1(\mathrm{m}\mathrm{o}\mathrm{d}n)$。（顺带一提，对于模意义下非零多项式方程，其解数可以参考拉格朗日定理 (Lagrange’s theorem)。）
• 剩下两个情况为：若 $n=4$，则 $(n-1)!\equiv 2\not\equiv -1(\mathrm{m}\mathrm{o}\mathrm{d}n)$；若 $n=2$，则 $(n-1)!\equiv 1\equiv -1(\mathrm{m}\mathrm{o}\mathrm{d}n)$。

---

高斯在此基础上进行了扩展，其结论是：当且仅当模 $n$ 意义下有原根时（即 $n=2,4,p^e,2p^e$，其中 $p$ 是奇质数，$e>0$ 时），${\prod }_{1\le u<n,\gcd (u,n)=1}u$ 在模 $n$ 意义下与 $-1$ 同余，否则该乘积与 $1$ 同余。 原证明略冗长，故这里略去。

顺带一提，若读者对模 $n$ 意义的缩系有一定了解，知道缩系可以表示成多个循环群的直积，则可以尝试这样进行充要证明：缩系中阶不超过 $2$ 的元素构成的子群可以表示成多个大小为 $2$ 的循环群的直积（设为 $k$ 个），其中每个直积因子中非单位元的元素会产生 $2^{k-1}$ 次贡献，因此当这样的直积因子有至少两个时，最终乘积必然是单位元（即 $1$），否则（若不存在至少两个这样的直积因子）此时各直积因子的大小互质，这意味着模 $n$ 意义下的原根存在。

注：这个证法可以证明更一般性的结论，即任意交换群的所有元素之积要么是单位元，要么存在唯一一个阶为 $2$ 的元素，且所有元素之积等于它。

1.4 伯努利公式 (Bernoulli’s formula)

伯努利公式也即等幂公式，主要是将 ${\sum }_{i=0}^{n-1}{i}^m$ 表示成关于 $n$ 的 $(m+1)$ 次多项式 $S_m(n)=\frac{1}{m+1}{\sum }_{i=1}^{m+1}\left(\genfrac{}{}{0ex}{}{m+1}{i}\right)B_{m+1-i}{n}^i$，其中 $B_i$ 是伯努利数列的第 $i$ 项。

考虑 $S_m(n)$ 看作关于 $m$ 的数列时对应的指数级生成函数 ${\sum }_{m\ge 0}\frac{S_m(n)}{m!}{x}^m$，可以得到

$$\begin{array}{rl}& \sum _{m\ge 0}\frac{S_m(n)}{m!}{x}^m\\ =& \sum _{i=0}^{n-1}\sum _{m\ge 0}\frac{(ix{)}^m}{m!}\\ =& \sum _{i=0}^{n-1}\exp (ix)\\ =& \frac{\exp (nx)-1}{x}/\frac{\exp (x)-1}{x}\text{,}\end{array}$$

其中 $\frac{\exp (nx)-1}{x}={\sum }_{m\ge 0}\frac{n^{m+1}}{(m+1)!}{x}^m$，而 $\frac{x}{\exp (x)-1}$ 则是伯努利数的指数型生成函数，即 $\frac{x}{\exp (x)-1}={\sum }_{m\ge 0}\frac{B_m}{m!}{x}^m$。

将伯努利数的生成函数代入即可求出 $S_m(n)$ 的公式，也可以由 $B_0=1$ 知道 $S_m(n)$ 是关于 $n$ 的 $(m+1)$ 次多项式。

注：不使用 $\exp (x)-1$ 而使用 $\frac{\exp (x)-1}{x}$ 作为分母是因为前者常数项为 $0$，而其作为形式幂级数的乘法逆元是不存在的。

1.5 拉格朗日插值法 (Lagrange interpolation)

对于一个已知次数为 $m$ 的多项式 $F(x)$，若已知对于 $(x_i,y_i)$ $(i=0,1,\dots ,m)$ 满足 $F(x_i)=y_i$，则 $F(x)={\sum }_{i=0}^m{y}_i{\ell }_i(x)$，其中 ${\ell }_i(x)={\prod }_{0\le j\le m,j\ne i}\frac{x-x_j}{x_i-x_j}$ 是满足 ${\ell }_i(x_i)=1$ 且 ${\ell }_i(x_j)=0$ $(j\ne i)$ 的插值基函数。 它的存在性和唯一性是显然的，故这里略去证明。

---

当 $x_i=i$ ($i=0,1,\dots ,m$) 时，上述式子可以化简为

$$F(x)=\sum _{i=0}^{m}F(i)\frac{(-1{)}^{m-i}{\prod }_{j=0}^m(x-j)}{i!(m-i)!(x-i)}=\sum _{i=0}^{m}F(i)(-1{)}^{m-i}(\genfrac{}{}{0ex}{}{x}{i}\left)\right(\genfrac{}{}{0ex}{}{x-i-1}{m-i})\text{.}$$

这意味着此时的插值基函数 ${\ell }_i(x)$ 在 $x$ 取整数时必然也是整数，那么即使 $F(x)$ 的系数不能在模 $p^e$ 意义下表示（即存在 $F(x)$ 的某个系数，满足其最简分数表示里分母是 $p$ 的倍数），只要所有相关的 $x$ 和 $F(x)$ 是整数，那么我们就可以在预处理的基础上，使用拉格朗日插值法在 $\mathcal{O}(m)$ 次大数乘法的时间内求出一个 $F(x)$。

令 $k!=p^{d_k}{r}_k$ $(d_k\ge 0,\gcd (p,r_k)=1)$，这里需要预处理 $(F(i)\mathrm{m}\mathrm{o}\mathrm{d}{p}^e)$, $d_i$ 和 $(r_i^{-1}\mathrm{m}\mathrm{o}\mathrm{d}{p}^e)$ $(i=0,1,\dots ,m)$。此处

• $F(i)$ 的计算代价需要根据其定义确定，暂时待定。
• 计算 $d_i$ 的总代价是 $\mathcal{O}\left({\nu }_p(m!)+m\right)=\mathcal{O}(m)$ 次大数除小数。
• 若先计算所有 $r_i$，再进行一次辗转相除得到 $r_m^{-1}$，然后依次得到所有 $r_i^{-1}$，那么只需要 $\mathcal{O}(m)$ 次大数乘法和一次求逆的代价。
• 利用辗转相除法进行一次求逆的代价是 $\mathcal{O}(e{\log }_{2}p)$ 次大数乘法和大数除法。

2. 计算方法

经过 Min_25 改进的计算方法已经变得相对来说浅显易懂，适合非数学专业的读者阅读思考。与 Andrew Granville 的方法相比，在此基础上增加的优化技巧也变得相对来说容易接受。

2.1 问题转化

对于整数 $u\ge 0$ 和 $0\le v<p$，定义

$$((up+v)!{)}_p=\prod _{0\le x\le up+v,\gcd (x,p)=1}x=\left(\prod _{i=0}^{u-1}\prod _{j=1}^{p-1}(ip+j)\right)\prod _{j=1}^v(up+j)\text{,}$$

则有

$$n!=p^{{\nu }_p(n!)}\prod _{i\ge 0}{\left(\left(\lfloor \frac{n}{p^{i+1}}\rfloor p+n_i\right)!\right)}_p\text{,}$$

其中每个 $((up+v)!{)}_p$ 都是与 $p$ 互质的值，因此计算 $\left(\genfrac{}{}{0ex}{}{n}{m}\right)\mathrm{m}\mathrm{o}\mathrm{d}{p}^e$ 可以化为计算 $\mathcal{O}({\log }_{p}n)$ 次 $\left(((up+v)!{)}_p\mathrm{m}\mathrm{o}\mathrm{d}{p}^{e-{\nu }_p\left(\left(\genfrac{}{}{0ex}{}{n}{m}\right)\right)}\right)$、相应次大数乘法，以及一次求逆。对于 ${\nu }_p\left(\left(\genfrac{}{}{0ex}{}{n}{m}\right)\right)\ge e$ 的情况，可以直接得出 $\left(\genfrac{}{}{0ex}{}{n}{m}\right)\mathrm{m}\mathrm{o}\mathrm{d}{p}^e=0$。

---

注意到

$$\prod _{j=1}^v(up+j)=\frac{(up{)}^{\overline{v+1}}}{up}=\sum _{k=0}^v\left[\genfrac{}{}{0ex}{}{v+1}{k+1}\right](up{)}^k\text{,}$$

而 $((up{)}^k\mathrm{m}\mathrm{o}\mathrm{d}{p}^e)$ 在 $k\ge p$ 时一定为 $0$，故计算一段连续的与 $p$ 互质的数之积可以在 O ( min ⁡ { p , e } ) \mathcal{O}(\min\{p, e\}) O(min{p,e}) 次大数乘法运算内完成。

同样地，我们也可以得到

∏ i = 0 u − 1 ∏ j = 1 p − 1 ( i   p + j ) ≡ ∏ i = 0 u − 1 ( ∑ k = 0 min ⁡ { p , e } − 1 [ p k + 1 ] ( i   p ) k ) ( m o d p e ) . \prod_{i = 0}^{u - 1}{\prod_{j = 1}^{p - 1}{(i\,p + j)}} \equiv \prod_{i = 0}^{u - 1}{\left(\sum_{k = 0}^{\min\{p, e\} - 1}{{p \brack k + 1} (i\,p)^k}\right)} \pmod{p^e}\text{.} i=0∏u−1​j=1∏p−1​(ip+j)≡i=0∏u−1​⎝⎛​k=0∑min{p,e}−1​[k+1p​](ip)k⎠⎞​(modpe).

2.2 优化分析

注意到上一小节末尾的式子是最难计算的部分，暴力计算它需要 O ( u min ⁡ { p , e } ) \mathcal{O}(u \min\{p, e\}) O(umin{p,e}) 次大数乘法和大数加法，而 $u$ 本身便可以是一个大数字，因此我们需要优化这部分的计算量。

---

注意到 $\gcd ([\genfrac{}{}{0ex}{}{p}{1}],p)=\gcd ((p-1)!,p)=1$，故 $\left({\left[\genfrac{}{}{0ex}{}{p}{1}\right]}^{-1}\mathrm{m}\mathrm{o}\mathrm{d}{p}^e\right)$ 存在，则

∏ i = 0 u − 1 ( ∑ k = 0 min ⁡ { p , e } − 1 [ p k + 1 ] ( i   p ) k ) ≡ [ p 1 ] u ∏ i = 0 u − 1 ( 1 + ∑ k = 1 min ⁡ { p , e } − 1 [ p k + 1 ] [ p 1 ] ( i   p ) k ) ( m o d p e ) . \prod_{i = 0}^{u - 1}{\left(\sum_{k = 0}^{\min\{p, e\} - 1}{{p \brack k + 1} (i\,p)^k}\right)} \equiv {p \brack 1}^{u} \prod_{i = 0}^{u - 1}{\left(1 + \sum_{k = 1}^{\min\{p, e\} - 1}{\frac{{p \brack k + 1}}{{p \brack 1}} (i\,p)^k}\right)} \pmod{p^e}\text{.} i=0∏u−1​⎝⎛​k=0∑min{p,e}−1​[k+1p​](ip)k⎠⎞​≡[1p​]ui=0∏u−1​⎝⎛​1+k=1∑min{p,e}−1​[1p​][k+1p​]​(ip)k⎠⎞​(modpe).

---

令 A ( x ) = ∑ k = 1 min ⁡ { p , e } − 1 a k x k A(x) = \sum_{k = 1}^{\min\{p, e\} - 1}{a_k x^k} A(x)=∑k=1min{p,e}−1​ak​xk，其中 $a_k=\left({\left[\genfrac{}{}{0ex}{}{p}{1}\right]}^{-1}\left[\genfrac{}{}{0ex}{}{p}{k+1}\right]\right)\mathrm{m}\mathrm{o}\mathrm{d}{p}^e$ 为整数。注意，这里 $A(x)$ 是定义在有理数域上的多项式，而有理数域上多项式的系数是可以表示成最简分数的有理数，但是分母可能是 $p$ 的倍数。

考虑形式幂级数 $B_u(x)={\prod }_{i=0}^{u-1}(1+A(ix))$，那么有

B u ( x ) = ∏ i = 0 u − 1 ( 1 + ∑ k = 1 min ⁡ { p , e } − 1 a k i k x k ) = 1 + ∑ k ≥ 1 t u , k x k , B_u(x) = \prod_{i = 0}^{u - 1}{\left(1 + \sum_{k = 1}^{\min\{p, e\} - 1}{a_k i^k x^k}\right)} = 1 + \sum_{k \geq 1}{t_{u, k} x^k}\text{,} Bu​(x)=i=0∏u−1​⎝⎛​1+k=1∑min{p,e}−1​ak​ikxk⎠⎞​=1+k≥1∑​tu,k​xk,

这里 $t_{u,k}$ 一定是整数，因为 $a_k{i}^k$ 是整数，这和 $x$ 的取值范围无关。

注意到我们希望计算的是 $(B_u(p)\mathrm{m}\mathrm{o}\mathrm{d}{p}^e)$，而当 $k\ge e$ 时，$p^e|p^k$ 且 $(t_{u,k}\mathrm{m}\mathrm{o}\mathrm{d}{p}^e)$ 一定可以表示出来，那么

$$B_u(p)\equiv 1+\sum _{k=1}^{e-1}{t}_{u,k}{p}^k(\mathrm{m}\mathrm{o}\mathrm{d}{p}^e)\text{,}$$

这启发我们寻找 $t_{u,k}$ 的表达式，从而计算出所求的式子。

---

令 $C(x)=\ln (1+A(x))={\sum }_{k\ge 1}\frac{(-1{)}^{k-1}}{k}{A}^k(x)={\sum }_{k\ge 1}{c}_k{x}^k$（注意 $c_k$ 是有理数，可能不是整数），于是有

$$\ln B_u(x)=\sum _{i=0}^{u-1}\ln (1+A(ix))=\sum _{i=0}^{u-1}\sum _{k\ge 1}{c}_k(ix{)}^k=\sum _{k\ge 1}{c}_k{s}_{u,k}{x}^k\text{,}$$

这里 $s_{u,k}={\sum }_{i=0}^{u-1}{i}^k$ 是一个关于 $u$ 的 $(k+1)$ 次多项式。

由

$$B_u^{\prime }(x)={\exp }^{\prime }(\ln B_u(x))=\exp (\ln B_u(x)){\ln }^{\prime }{B}_u(x)=B_u(x){\ln }^{\prime }{B}_u(x)$$

可知

$$\sum _{k\ge 1}k{t}_{u,k}{x}^{k-1}=\left(\sum _{k\ge 1}{t}_{u,k}{x}^k\right)\left(\sum _{k\ge 1}k{c}_k{s}_{u,k}{x}^{k-1}\right)\text{,}$$

取其中 $x^{k-1}$ 的系数可知 $t_{u,k}=\frac{1}{k}{\sum }_{i=1}^k{t}_{u,k-i}i{c}_i{s}_{u,i}$。

---

由 $s_{u,k}$ 是关于 $u$ 的多项式，我们可以归纳得到 $t_{u,k}$ 也是关于 $u$ 的多项式，从而确定 $B_u(p)$ 也是关于 $u$ 的多项式。

具体来说，令 ${\mathrm{d}\mathrm{e}\mathrm{g}}_u(f)$ 表示多项式 $f$ 关于 $u$ 的次数，则 ${\mathrm{d}\mathrm{e}\mathrm{g}}_u(t_{u,k})={\max }_{i=1}^k({\mathrm{d}\mathrm{e}\mathrm{g}}_u(t_{u,k-i})+{\mathrm{d}\mathrm{e}\mathrm{g}}_u(s_{u,i}))$，于是可以归纳证明 ${\mathrm{d}\mathrm{e}\mathrm{g}}_u(t_{u,k})=k{\mathrm{d}\mathrm{e}\mathrm{g}}_u(s_{u,1})=2k$，而 ${\mathrm{d}\mathrm{e}\mathrm{g}}_u(1+{\sum }_{k=1}^{e-1}{t}_{u,k}{p}^k)={\max }_{k=1}^{e-1}{\mathrm{d}\mathrm{e}\mathrm{g}}_u(t_{u,k})=2e-2$。

虽然这个多项式的系数可能不是整数，但当 $u$ 取非负整数时，多项式的值都是整数，故可以使用拉格朗日插值法计算一个 $B_u(p)$ 在模 $p^e$ 意义下的值。

---

此外，注意到当 $u\ge p^{e-1}$ 时，$B_u(p)$ 里会出现多个 ${\prod }_{1\le x<p^e,\gcd (x,p^e)=1}x$，而根据威尔逊定理的扩展，这个部分在模 $p^e$ 意义下只会是 $\pm 1$。不妨记这个值为 $\delta$，那么 $B_u(p)\equiv {\delta }^{\lfloor u/p^{e-1}\rfloor }{B}_{u\mathrm{m}\mathrm{o}\mathrm{d}{p}^{e-1}}(p)(\mathrm{m}\mathrm{o}\mathrm{d}{p}^e)$。为了便于计算，不妨将这里的 $p^{e-1}$ 换为 $2p^{e-1}$，则前面的正负号计算可以省去。

将插值法用到的 $u$ 限制在 $[0,2p^{e-1})$ 的好处是，计算 ${\nu }_p(u-i)$ $(i=0,1,\dots ,2e-2)$ 的代价可以从 $\mathcal{O}({\log }_{p}n+e)$ 次大数除小数变成 $\mathcal{O}(e)$ 次大数除小数。

---

还有一处需要分析的是计算 ${\left[\genfrac{}{}{0ex}{}{p}{1}\right]}^u$，注意到多次计算 $((up+v)!{)}_p$ 时的这部分是可以合并的。合并之后可以看出，$\left(\genfrac{}{}{0ex}{}{n}{m}\right)$ 中所需计算的部分是 ${\left[\genfrac{}{}{0ex}{}{p}{1}\right]}^{{\nu }_p\left(\left(\genfrac{}{}{0ex}{}{n}{m}\right)\right)}$，其幂次的大小是 $\mathcal{O}({\log }_{p}n)$ 的，而且只有当 ${\nu }_p\left(\left(\genfrac{}{}{0ex}{}{n}{m}\right)\right)<e$ 时才会产生计算，于是暴力计算即可。

2.3 复杂度计算

最后，我们来算一下大数运算的复杂度，这里略去小数运算是因为后者的复杂度远小于前者的复杂度。

要说明的一点是，我们默认存储的进制是 $p$ 的幂次，大数对 $p^e$ 取模可以 $\mathcal{O}(1)$ 实现，故将一次模乘的代价视为一次大数乘法的代价。（但模意义下求逆可能需要多次大数除法。）如果要考虑大数取模的复杂度，可以直接视所有单次大数乘除和一次大数取模同阶。

---

本文中所有涉及到的部分可以做如下分析：

• $n$ 和 $m$ 视为同阶，大小均为 $\mathcal{O}(n)$，而它们转化为 $p$ 进制的代价是 $\mathcal{O}({\log }_{p}n)$ 次大数除小数。
• 为了高效计算一段连续的与 $p$ 互质的数之积，需要预处理 $\left[\genfrac{}{}{0ex}{}{i}{j}\right]$ ( i = 1 , 2 , … , p ; j = 1 , 2 , … , min ⁡ { i , e } ) (i = 1, 2, \ldots, p; j = 1, 2, \ldots, \min\{i, e\}) (i=1,2,…,p;j=1,2,…,min{i,e})，代价是 O ( p min ⁡ { p , e } ) \mathcal{O}(p \min\{p, e\}) O(pmin{p,e}) 次大数乘法和大数加法。
• 预处理 $B_i(p)$ $(i=0,1,\dots ,2e-2)$ 的代价是 O ( e min ⁡ { p , e } ) \mathcal{O}(e \min\{p, e\}) O(emin{p,e}) 次大数乘法和大数加法。
• 预处理插值系数的代价是 $\mathcal{O}(e+e{\log }_{2}p)$ 次大数乘法、$\mathcal{O}(e{\log }_{2}p)$ 次大数除法和 $\mathcal{O}(e)$ 次大数除小数。
• $B_u(p)$ 单次求值的代价是 $\mathcal{O}(e)$ 次大数加法、大数乘法和大数除小数。
• 计算 $\mathcal{O}({\log }_{p}n)$ 个涉及到的 $((up+v)!{)}_p$ 需要相应次 $B_u(p)$ 求值，以及额外的 $\mathcal{O}({\log }_{p}n)$ 次大数加法和大数乘法。

故总代价为

• 大数加法、大数除小数：比乘除次数少，故略去；
• 大数乘法： O ( p min ⁡ { p , e } + e min ⁡ { p , e } + e + e log ⁡ 2 p + e log ⁡ p n ) = O ( p   e + e log ⁡ p n ) \mathcal{O}(p \min\{p, e\} + e \min\{p, e\} + e + e \log_2{p} + e \log_p{n}) = \mathcal{O}(p\,e + e \log_p{n}) O(pmin{p,e}+emin{p,e}+e+elog2​p+elogp​n)=O(pe+elogp​n)；
• 大数除法：$\mathcal{O}(e{\log }_{2}p)$，或者视为常数次求逆操作。

3. 一些改进

本节主要介绍一些可能没什么用或者严重超纲的改进方法，供有兴趣的读者研究思考。（本节基本上是 Min_25 没有讨论的内容。）

3.1 减少预处理

当 $p$ 远小于 $e$ 时，可以省去无符号第一类斯特林数的预处理，但复杂度不变。

---

这一点比较显然，留给读者分类讨论。

3.2 降低多项式次数

当 $p$ 很大而导致 $B_u(x)$ 的前 $e$ 项系数（有理数）的最简分母大多数都与 $p$ 互质时，模 $p^e$ 意义下多项式的次数可能降低，虽然也不会降低太多就是了。

---

举例来讲，对于 $B_u(p)$ 中 $t_{u,k}{p}^k$ 的部分，若 $t_{u,k}={\sum }_{i=0}^{2k}{\mathrm{c}\mathrm{o}\mathrm{e}\mathrm{f}\mathrm{f}}_i{u}^i$ 且 ${\mathrm{c}\mathrm{o}\mathrm{e}\mathrm{f}\mathrm{f}}_i$ 都能在模 $p^{e-k}$ 意义下表示，那么根据欧拉定理 (Euler’s theorem) 的扩展形式，总存在 $i_0\le e-k$，使得对于任意 $u$， $u^{i_0+\delta }\equiv u^{i_0+(\delta \mathrm{m}\mathrm{o}\mathrm{d}\phi (p^{e-k}))}(\mathrm{m}\mathrm{o}\mathrm{d}{p}^{e-k})$，于是 $t_{u,k}$ 的次数就能直接降低到 min ⁡ { 2 k , i 0 + φ ( p e − k ) − 1 } \min\{2 k, i_0 + \varphi(p^{e - k}) - 1\} min{2k,i0​+φ(pe−k)−1}，从而降低整体的次数。

通过上述对局部的分析，我们发现只要系数能在模意义下表示，次数就能降低，甚至 $p$ 不用太大，例如当 $p>e$ 时，所有的系数一定都能表示出来。

注：笔者对系数能表示的证明有点长，故这里略去，推荐有兴趣的读者先证明 $C(x)$ 中的 $k{c}_k$ 是整数。

顺带一提，我们可以利用差分的技巧，求出 $B_u(p)$ 的多阶差分（显然差分值都是整数），从而确定最优的次数。

注：差分可以得出 $B_u(p)\equiv {\sum }_{k=0}^{{\mathrm{d}\mathrm{e}\mathrm{g}}_u(B_u(p))}{\mathrm{c}\mathrm{o}\mathrm{e}\mathrm{f}\mathrm{f}}_k\left(\genfrac{}{}{0ex}{}{u}{k}\right)$，其中 ${\mathrm{c}\mathrm{o}\mathrm{e}\mathrm{f}\mathrm{f}}_k$ 是整数。然而，从这个结论中我们只能知道，当 $p>2e-2$ 时 $t_{u,k}$ 一定是整数，若想把这个条件进一步扩大（例如 $p>e$），读者还是需要思考别的方法。

3.3 优化求逆过程

注意到只有求逆可能需要大数除法，而当 $p^e$ 很大时，求逆部分还有优化空间，甚至避免大数除法，虽然求逆已经不是主要运算量。

---

注意到求逆本身可以看作是模 $p^e$ 意义下的方程 $f(x)=ax-1$，当 $\gcd (a,p)=1$ 时，$f(x)\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}{p}^e)$ 有唯一解系 $(a^{-1}\mathrm{m}\mathrm{o}\mathrm{d}{p}^e)$。

不妨设 $f(x)\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}{p}^k)$ 的解是 $r_k$ $(k=1,2,\dots ,e)$，我们将 $(f(r_{2k})\mathrm{m}\mathrm{o}\mathrm{d}{p}^{2k})$ 在 $x=r_k$ 处泰勒展开，有

$$0\equiv f(r_{2k})\equiv f(r_k)+f^{\prime }(r_k)(r_{2k}-r_k)(\mathrm{m}\mathrm{o}\mathrm{d}{p}^{2k})\text{,}$$

故 $r_{2k}\equiv r_k-\frac{f(r_k)}{f^{\prime }(r_k)}(\mathrm{m}\mathrm{o}\mathrm{d}{p}^{2k})$。

由于 $p^k|f(r_k)$，则上式中的 $\frac{1}{f^{\prime }(r_k)}$ 只需要知道在模 $p^k$ 意义下的值即可，而 $a^{-1}\equiv r_k(\mathrm{m}\mathrm{o}\mathrm{d}{p}^k)$，于是有 $r_{2k}\equiv r_k-(a{r}_k-1)r_k\equiv r_k(2-a{r}_k)(\mathrm{m}\mathrm{o}\mathrm{d}{p}^{2k})$。

我们可以通过小数字的辗转相除求出 $r_1$，然后通过倍增（或者说牛顿迭代法）得到 $r_e$，这样做只需要 $\mathcal{O}({\log }_{2}p)$ 次小数字乘除和 $\mathcal{O}({\log }_{2}e)$ 次大数乘法，从而避免了大数除法。事实上，如果大数乘法使用快速傅里叶变换 (Fast Fourier transform) 实现，则实际的大数乘法运算量等价于 $\mathcal{O}(1)$ 次 $p^e$ 规模的大数乘法。

3.4 卷积加速预处理

当 $p$ 远大于 $e$ 和 ${\log }_{p}n$ 时，预处理的斯特林数里有很多之后用不上的值，若在使用快速傅里叶变换 (Fast Fourier transform) 实现大数乘法基础上再实现系数为大数的多项式乘法，则复杂度可以进一步优化，如果不是那么难写的话。

---

设模 $p^e$ 意义下的数字 $u={\sum }_{i=0}^{e-1}{u}_i{p}^i$ ($0\le u_i<p$)，而模 $p^e$ 意义下关于 $x$ 的次数为 ${\mathrm{d}\mathrm{e}\mathrm{g}}_x(F)$ 的多项式 $F(x)={\sum }_{i=0}^{{\mathrm{d}\mathrm{e}\mathrm{g}}_x(F)}{f}_i{x}^i={\sum }_{i=0}^{{\mathrm{d}\mathrm{e}\mathrm{g}}_x(F)}{\sum }_{j=0}^{e-1}{f}_{i,j}{p}^j{x}^i$。

注意到 $F(x)$ 可以看作是关于 $p$ 和 $x$ 的二维多项式，那么只需要（在复数域上）实现小数字 $f_{i,j}$ 的二维离散傅里叶变换 (Discrete Fourier transform)，那么就可以在 $\mathcal{O}({\mathrm{d}\mathrm{e}\mathrm{g}}_x(F)e{\log }_2({\mathrm{d}\mathrm{e}\mathrm{g}}_x(F)e))$ 的时间内完成多项式乘法，这里笔者为了简写，决定将一次乘法的代价视作 $\mathcal{O}({\mathrm{d}\mathrm{e}\mathrm{g}}_x(F){\log }_2({\mathrm{d}\mathrm{e}\mathrm{g}}_x(F)))$ 次大数加法和 $\mathcal{O}({\mathrm{d}\mathrm{e}\mathrm{g}}_x(F))$ 次大数乘法。

令 $F_v(x)={\prod }_{j=1}^v(x+j)$，则有 $F_{2v}(x)=F_v(x)F_v(x+v)$，注意到当 $v<p$ 时，$x$ 和 $(x+v)$ 不会同时是 $p$ 的倍数，于是我们无法从 $F_v(x)$ 的前 $e$ 项系数得到 $F_{2v}(x)$ 的前 $e$ 项系数，我们需要进一步的观察。

不妨设定阈值 $\alpha =\lfloor \sqrt{v}\rfloor$，那么有 $F_v(x)=\left({\prod }_{i=0}^{\alpha -1}{F}_{\alpha }(x+\alpha i)\right){\prod }_{j={\alpha }^2+1}^v(x+j)$，而求 $F_v(x)$ 在某个 $x$ 处的点值便可以转化为求 $F_{\alpha }(x)$ 在 $\alpha$ 个 $x$ 处的点值。由于这 $\alpha$ 个位置是等差数列，而且 $F_{\alpha }(x)$ 的次数是 $\alpha$，读者不难想到使用拉格朗日插值法快速求解这些点值。

举例来说，假设已知 $F_{\alpha }(\alpha i)$ $(i=0,1,\dots ,\alpha )$，现在要求 $F_{\alpha }(up+\alpha k)$ $(k=0,1,\dots ,\alpha -1)$，那么有

$$\begin{array}{rl}& F_{\alpha }(up+\alpha k)\\ =& \sum _{i=0}^{\alpha }{F}_{\alpha }(\alpha i)\prod _{0\le j\le \alpha ,j\ne i}\frac{up+\alpha k-\alpha j}{\alpha i-\alpha j}\\ =& \left(\sum _{0\le i\le \alpha ,i\ne k}\frac{(-1{)}^{\alpha -i}{F}_{\alpha }(\alpha i)}{i!(\alpha -i)!}\frac{\alpha }{up+\alpha (k-i)}\right)\left(\prod _{0\le j\le \alpha }\frac{up+\alpha (k-j)}{\alpha }\right)\\ & +\frac{(-1{)}^{\alpha -k}{F}_{\alpha }(\alpha k)}{k!(\alpha -k)!}\prod _{0\le j\le \alpha ,j\ne k}\frac{up+\alpha (k-j)}{\alpha }\text{,}\end{array}$$

这里明显出现了一个卷积式，用多项式乘法加速即可。

注：由于 $\alpha \le v<p$，故上式中的分母都是有逆元的，而且可以沿用求阶乘逆元的方法将上式中所需的求逆操作优化到常数次（甚至总数也是常数次）。

同理，我们可以利用拉格朗日插值法辅助倍增求出 $F_{\alpha }(\alpha i)$ $(i=0,1,\dots ,\alpha )$。例如，根据 $F_d(\alpha i)$ $(i=0,1,\dots ,d)$，我们可以计算出 $F_d(\alpha i+d)$, $F_d(\alpha (i+d))$ 和 $F_d(\alpha (i+d)+d)$，从而得到 $F_{2d}(\alpha k)$ $(k=0,1,\dots ,2d)$，类似地可以得到 $F_{2d+1}$ 的点值。（如果对这部分有困惑，可以参考 fjzzq2002 的翻译文章。）

这样做所需的等价大数加法和大数乘法运算量是

$$\begin{array}{rlrl}{T}_{add}(\alpha )& =T_{add}(\alpha /2)+\mathcal{O}(\alpha {\log }_2\alpha )& \Rightarrow & T_{add}(\alpha )=\mathcal{O}(\alpha {\log }_2\alpha )\\ T_{mul}(\alpha )& =T_{mul}(\alpha /2)+\mathcal{O}(\alpha )& \Rightarrow & T_{mul}(\alpha )=\mathcal{O}(\alpha )\text{.}\end{array}$$

通过这样的改进（使用快速傅里叶变换和倍增法等等），我们可以将整体所需的运算量改进到 $\mathcal{O}((e+{\log }_{p}n)\sqrt{p}{\log }_{2}p)$ 次大数加法、$\mathcal{O}((e+{\log }_{p}n)\sqrt{p})$ 次大数乘法，附加常数次求逆。

注：本小节讨论的内容是基于 p ≥ max ⁡ { e , log ⁡ p n } p \geq \max\{e, \log_p{n}\} p≥max{e,logp​n} 分析的，而得到的结果只有在 $p$ 远大于 $e$ 和 ${\log }_{p}n$ 时才比第 2 节中的结果更优。

4. 一些讨论

根据本文得出的结论，$((2n-1)!!\mathrm{m}\mathrm{o}\mathrm{d}{2}^e)$ 将会是一个关于 $n$ 的次数不超过 $(2e-2)$ 的多项式，同理 $({\sum }_{i=1}^n(2i-1)!!\mathrm{m}\mathrm{o}\mathrm{d}{2}^e)$ 将会是关于 $n$ 的次数不超过 $(2e-1)$ 的多项式。

---

待补充