0. 介绍
组合数取模，即给定 (n,m,p)，求
下面按数据范围讨论一些解法：
| (p) / (n,m) | (n,mle 200) | (n,mle 10^6) | (n,mle 10^{18}) |
|---|---|---|---|
| (ple 10^6) 且 (p) 为素数 | 杨辉三角暴力递推 | Lucas 定理 | Lucas 定理 |
| (p>10^9) 且 (p) 为素数 | 杨辉三角暴力递推 | 预处理阶乘及其逆元 | exLucas |
| (p=p_1p_2p_3cdots p_k) 且 (p_{1cdots k}le 2000) 为互不相同的素数 | 杨辉三角暴力递推 | Lucas 定理 CRT 合并 | Lucas 定理 CRT 合并 |
| (p) 无特殊限制 | 杨辉三角暴力递推 | exLucas | exLucas |
（(p) 的第二档本来想表达 (p) 远大于 (n) 的，但这样 (n,mle 10^{18}) 的档就不可做了 qwq）
（(p) 的第三档就是 (p) 为 square-free number）
（CRT 指中国剩余定理）
时空复杂度：
| 做法 | 预处理时间复杂度 | 询问时间复杂度 | 空间复杂度 |
|---|---|---|---|
| 杨辉三角暴力递推 | (O(n^2)) | (O(1)) | (O(n^2)) |
| 预处理阶乘及其逆元 | (O(n)) | (O(1)) | (O(n)) |
| Lucas 定理 | 预处理 (1sim p) 的组合数，取决于所用的算法 | (O(plog n)) | 预处理 (1sim p) 的组合数，取决于所用的算法 |
| Lucas 定理 CRT 合并 | 预处理 (1sim p) 的组合数，取决于所用的算法 | (O(Plog n+plog P)) | 预处理 (1sim p) 的组合数，取决于所用的算法 |
| exLucas | 无预处理，(O(1)) | (O(sum p^klog n(log n-k)+plog P)) | (O(1)) |
| 多组询问 exLucas（模数相同） | (O(sum p^k)) | (O(log p^kn)) | (O(1)) |
在 Lucas 定理 CRT 合并中， (P) 是模数，(p) 是最大质因子
1. 杨辉三角暴力递推
我们有
（考虑 (n) 里选 (m)，讨论选不选第一个即可证明）
(O(n^2)) 递推即可 .
2. 预处理阶乘及其逆元
众所周知
(O(n)) 预处理阶乘及其逆元即可
3. Lucas 定理
Lucas 定理：
形式 一
设 (n,m) 在 (p)（(p) 是质数）进制下表示为
[overline{n_kn_{k-1}n_{k-2}dots n_1}\,,overline{m_km_{k-1}m_{k-2}dots m_1} ]其中 (0le n_i<p)，(0le m_i<p) .
则有：
[dbinom{n}{m}equiv dbinom{n_k}{m_k}dbinom{n_{k-1}}{m_{k-1}}cdotsdbinom{n_1}{m_1}pmod p ]
形式二（简单）
[C^n_mequiv C^{nmod p}_{mmod p}cdot C^{lfloor n/p floor}_{lfloor m/p floor} ]
Lucas 定理的证明
设 (xinmathbb N^+)，且 (x<p)，则：
由于 (x<p) 且 (p) 是素数，所以 (x) 存在模 (p) 意义下的逆元，因此：
因为 (pcdot x^{-1}cdot C^{x-1}_{p-1}equiv 0pmod p) ，故 (C^x_pequiv0pmod p).
由二项式定理，
除了 (i=0,i=p) 的项数，别的模 (p) 均为 (0)，所以：
现在我们设：
从而：
再由二项式定理有
而同时又有：
注意到 (ip+j) 正好能取到 (0) 到 (m) 内所有整数一次，所以枚举 (k=ip+j)，得
结合二项式定理，得
对比系数，得：
命题获证 .
用形式一可以写出一个递推的代码：
int Lucas(int n,int m,int p)
{
for (int i=0;i<p;i++)
{
C[i][0]=1;
for (int j=1;j<=i;j++)
{
C[i][j]=C[i-1][j-1]+C[i-1][j];
if (C[i][j]>=p) C[i][j]-=p;
}
} // init
int ans=1;
while (n||m){ans=ans*C[n%p][m%p]; n/=p; m/=p;}
return ans;
}
而用形式二则可以写出一个递归的代码：
ll lucas(ll n,ll m){return (n<m)?0:(n<p)?C(n,m):C(n%p,m%p)*lucas(n/p,m/p)%p;}
4. Lucas 定理 CRT 合并
令
则
用 Lucas 定理求组合数，然后 CRT 合并即可
正确性显然
5. exLucas
先把 (p) 分解质因数，令
用同样的做法，列出同余方程组，求出 (dbinom nmmod p_i^{k_i}) 然后 CRT 合并 .
那么问题来了，(dbinom nmmod p_i^{k_i}) 怎么求？
我们知道，(dbinom nm=dfrac{n!}{m!(n-m)!})，由于逆元可能不存在，考虑提取出质因子 (p_i) 来
其中 (A) 为 (n!) 中质因子 (p_i) 的次数，(B,C) 同理 .
此时分母就存在逆元了，现在问题又转化成了求形如
先考虑如何计算 (n!mod p_k) .
举个例子（(22!mod 3^2)）
可以看出上式分为三个部分：
第一部分是 (p) 的幂，次数是小于等于 (n) 的 (p) 的倍数的个数（即 (leftlfloordfrac{n}{p}
ight
floor)）.
第二部分是一个阶乘 (7!)，即 (leftlfloordfrac{n}{p} ight floor!)，可以递归求
第三部分是 (n!) 中与 (p) 互质的部分的乘积，这一部分具有如下性质（加个 (p^k)）：
一般的，有
整块快速幂小块暴力即可求出 (n!mod p^k)
回到原问题，求
分母全部由上述第一部分和第二部分贡献（第三部分和 (p) 互质）。而递归计算第二部分的时候已经除去了第二部分中的因子 (p)，所以最终的答案就是上述第二部分递归返回的结果和第三部分的乘积（与第一部分无关） .
然后一层层回去即可求出组合数取模 .
typedef long long ll;
ll qpow(ll a,ll n,ll p)
{
ll ans=1;
while (n)
{
if (n&1) ans=ans*a%p;
a=a*a%p; n>>=1;
} return ans;
}
namespace extra_lucas_detail
{
ll fac(ll n,ll p,ll pk)
{
if (!n) return 1;
ll ans=1;
for (int i=1;i<pk;i++)
if (i%p) ans=ans*i%pk;
ans=qpow(ans,n/pk,pk);
for (int i=1; i<=n%pk;i++)
if (i%p) ans=ans*i%pk;
return ans*fac(n/p,p,pk)%pk;
}
ll exgcd(ll a,ll b,ll& x,ll& y)
{
if (!b){x=1; y=0; return a;}
ll d=exgcd(b,a%b,x,y);
int z=x; x=y; y=z-y*(a/b);
return d;
}
ll inv(ll a,ll P){ll x,y,d=exgcd(a,P,x,y); return (d==1)?x%P:-1;}
ll C(ll n,ll m,ll p,ll pk)
{
if (n<m) return 0;
ll A=fac(n,p,pk),B=fac(m,p,pk),C=fac(n-m,p,pk),cnt=0;
for (ll i=n;i;i/=p) cnt+=i/p;
for (ll i=m;i;i/=p) cnt-=i/p;
for (ll i=n-m;i;i/=p) cnt-=i/p;
return A*inv(B,pk)%pk*inv(C,pk)%pk*qpow(p,cnt,pk)%pk;
}
}
ll exlucas(ll n,ll m,ll p)
{
using namespace extra_lucas_detail;
ll ans=0,P=p;
for (int i=2;(p>1)&&(i*i<=p);i++)
{
ll pk=1;
while (!(p%i)){p/=i; pk*=i;}
if (pk>1)
{
ll now=C(n,m,i,pk);
ans=(ans+now*(P/pk)%P*inv(P/pk,pk)%P)%P;
}
}
if (p>1)
{
ll now=C(n,m,p,p);
ans=(ans+now*(P/p)%P*inv(P/p,p)%P)%P;
} return (ans%P+P)%P;
}