FastMod代码来源。证明和代码分开，可以根据自己的需要跳转。
Lucas定理
$${n \choose m} \mod p = {\lfloor \frac{n}{p} \rfloor \choose \lfloor \frac{m}{p} \rfloor} {\langle n \rangle _p \choose \langle m \rangle _p} \mod p$$
上式当$p$为素数时成立。
此时，若$p$为一个较小的素数，而$n,m$为一个较大的数（指不能直接通过预处理阶乘及其逆元的大数，比如$1e18$），那么我们就可以递归地求解该结果。
| |
根据上述代码，我们只要预处理出$p$以内的阶乘及其逆元即可，时间复杂度$O(p+\log{m})$。
证明
其实在上面的递归过程已经提示了一件事：$p$进制分解。
考虑$n,m$在$p$进制的表达式
$$n=\sum_{i=0}^{k}{a_kp^k},m=\sum_{i=0}^{k}{b_kp^k}$$
那么Lucas定理实际就是：
$${n \choose m}\equiv{a_0 \choose b_0}{a_1 \choose b_1}\dots{a_k \choose b_k}\pmod{p}$$
现在先来看一个简单的式子：
$$(a+b)^p=\sum_{i=0}^{p}{p \choose i}a^ib^{p-i}\equiv{a^p+b^p}\pmod{p}$$
因为其余的部分都有因子$p$
那么我们可以得到：
$$(1+x)^n\equiv(1+x)^{\sum a_ip^i}\equiv\prod\left(1+x^{p^i}\right)^{a_i}\pmod{p}$$
现在考虑$x^m$的系数，左边是${n \choose m}$，右边是$\prod{a_i \choose b_i}$，已经是Lucas的表达式啦。
拓展问题 - exLucas定理
如果上式中$p$为合数，同时$n,m$依旧很大，怎么求解呢？
由中国剩余定理可以知道，
如果合数$p$的唯一分解为$\prod_{i=1}^{s}{p_i}^{\alpha_i}$，则有： $$ x \equiv {n \choose m} \pmod{p} \Rightarrow \begin{cases} x & \equiv & {n \choose m} \pmod{p_1^{\alpha_1}} \newline x & \equiv & {n \choose m} \pmod{p_2^{\alpha_2}} \newline & \vdots & \newline x & \equiv & {n \choose m} \pmod{p_s^{\alpha_s}} \newline \end{cases} $$
所以只要求出$p$的唯一分解，再分别解出各个同余式，通过CRT合并答案就能完成全部解答。
现在的问题就转变为，如何求出${n \choose m} \mod {p_i}^{\alpha_i}$。
我们知道${n \choose m}=\frac{n!}{m!(n-m)!}$，如果$m!$与$p^{\alpha}$互素，那么只要拓展欧几里得求出阶乘逆元就能直接算出结果。但是，先不论$n,m$的大小是否支持我们直接求出阶乘。它们是否互素呢？显然是否定的。
不过，如果能分离出不互素的部分和互素的部分，我们可以就求逆元得到答案了。
注意到，$p^{\alpha}$只有一种质因子$p$，那么阶乘中与之不互素的部分显然都是与$p$不互素的。这种不互素的数的个数也可以轻松算出，即${pot_p(n!)}={\sum_{i=1}^{\infty}{\lfloor \frac{n}{p^i} \rfloor}}$。然后在模$p^{\alpha}$意义下所有数在$Z_{p^{\alpha}}$中，（注意$p^{\alpha}$比较小）也就是说大于$p^{\alpha}$的数通过取模必然在$[0,p-1]$中，且很容易知道这样的数是循环出现的。
所以现在的思路就是，预处理$p^{\alpha}$内所有与$p$互素的数的乘积（令为$base$），那么$n!$中必然有$\lfloor \frac{n}{p^{\alpha}} \rfloor$个这样的乘积，剩余部分为$n%p^{\alpha}$内与$p$互素的数的乘积。
$$ fac(n)=base^\frac{n}{p^{\alpha}} fac(n \bmod p^{\alpha}) $$
总结一下就是，现在我们把阶乘在$O({p^{\alpha}}+\log_{p}{n})$时间上，分成了互素数乘积fac(n)和不互素的p^pot两部分。
所以，上代码了。
| |
那么，现在我们就能得到： $${n \choose m} \equiv \frac{fac(n)}{fac(m) \times fac(n-m)} \times p^{pot_p(n)-pot_p(m)-pot_p(n-m)} \pmod{p^{\alpha}}$$
之后就是简单的CRT啦，这个就不单独上代码了。
直接给出全部代码吧：
因为最近沉迷压行，可能看起来有点不舒服，
但是我舒服就行但是可以复制到本地去格式化，所以我就不改啦。
| |
一些优化（LibreOJ #181. 二项式系数）
上面给的代码只能过一些数据比较水的单组测试样例，因为是板子题所以没有优化。
要过一些其它的题比如LibreOJ #181. 二项式系数，这种数据比较毒瘤的是肯定不行的。
以下优化基于模数确定，多组输入的情况。
预处理部分
在求分离阶乘的时候，我们其实可以通过前缀和预处理出，前k个数中与$p^{\alpha}$互素的数的乘积$sum[k]$，以及$base$的$\lambda$次幂$base_{sum}[\lambda]$。
由于$base$与$p^{\alpha}$互素，通过欧拉定理可知$base^{\varphi(p^{\alpha})}\equiv 1\pmod{p^{\alpha}}$。
$\varphi(p^k)=p^{k-1}\varphi(p)=p^{k-1}(p-1)$
所以预处理时，$sum$预处理到$p^{\alpha}$，$base_{sum}$预处理到$p^{k-1}(p-1)$即可。
再一个是CRT函数。
设$m_1,m_2,\dots,m_k$为$k$个两两互素的正整数，$m=\prod{m_i}$。
令$m=m_iM_i$，则同余方程组 $$\begin{cases}x &\equiv & b_1 \pmod{m_1} \newline x &\equiv & b_2 \pmod{m_2} \newline & \vdots & \newline x &\equiv & b_k \pmod{m_k} \end{cases}$$
有唯一解 $$x\equiv {\sum_{i=1}^{k}{b_iM_iM_i^{-1}}}\pmod{m},M_iM_i^{-1}\equiv1\pmod{m_i}$$
由于需要多次调用CRT，我们可以直接预处理$M_iM_i^{-1}$这一部分。
再预处理完这些之后，我们就来改写fac函数，使之同时求出pot。
| |
然后是exLucas函数
| |
复杂度分析
我们来算算复杂度
- 唯一分解加上预处理 $$\log{P}\sum(\varphi(p^k)+p^k+\log{p^k})$$
- 查询 $$T\log{P}(\log_p{(n)}+\log_p{(m)}+\log_p{(n-m)}+\log(pot))$$
总共大概$1e8$和一个不太大（实际很大）的常数，然后本地跑了3.6秒，一直TLE。。。
感觉太玄学了就去看了大佬提交的代码，然后就看到了一个很玄学的优化（取模）。
大佬的玄学优化（取模）
| |
| |
个人理解：除法的效率很低，用位运算代替除法提高效率的优化（？）
这个是我自己一直在用的，比较常规。
| |
我搜索了一下c++中各种运算的效率
来源：CSDN博客，不保证正误，但是至少知道除法很慢就对了。
移位 > 赋值 > 大小比较 > 加法 > 减法 > 乘法 > 取模 > 除法。
由于模数确定，我们处理一下再用位运算右移代替除法会比正常取模快。
AC代码
LibreOJ 提交链接 加了快读，以及压行，所以有点丑，但应该不影响阅读。
| |