需要求 n,m≤lim 的一些 Cnm 模 p 的结果,但是递推预处理组合数时间复杂度无法接受时,该怎么办呢?
0x0 普通方法
最基础的方法是预处理出阶乘和阶乘的逆元,直接套公式计算:
inline void init()
{
fra[0]=1;
for(int i=1;i<=lim;i++)
{
fra[i]=1ll*fra[i-1]*i%p;
}
inv[lim]=qpow(fra[lim],p-2);
for(int i=lim;i>=1;i--)
{
inv[i-1]=1ll*inv[i]*i%p;
}
}
inline int C(int n,int m)
{
return n<m?0:1ll*fra[n]*fra[n-m]%p*fra[m]%p;
}
这样做的话单次时间复杂度是 O(1) 的,预处理时间复杂度则是 O(lim)。
1x0 特殊手段
普通方法已经能满足大多数情况,但是有两种特殊情况除外。
1x1 lim≥p,但 p 是质数的情况(Lucas)
当 lim≥p 时,所有 i≥p 的 i! 都不与 p 互质,没有模 p 意义下的逆元,普通方法就行不通了。
这时可以借助 Lucas 定理:
Cnm≡Cnmodpmmodp⋅C⌊pn⌋⌊pm⌋(modp)
来实现递归求解:
inline void init()
{
int lim2=min(lim,p-1);
fra[0]=1;
for(int i=1;i<=lim2;i++)
{
fra[i]=1ll*fra[i-1]*i%p;
}
inv[lim2]=qpow(fra[lim2],p-2);
for(int i=lim2;i>=1;i--)
{
inv[i-1]=1ll*inv[i]*i%p;
}
}
inline int C(int n,int m)
{
return n<m?0:1ll*fra[n]*inv[m]%p*inv[n-m]%p;
}
int lucas(int n,int m)
{
return m==0?1:1ll*C(n%p,m%p)*lucas(n/p,m/p)%p;
}
这样做的话单次时间复杂度是 O(logpm) 的,预处理时间复杂度则是 O(min(lim,p))。
定理证明:
【引理 1】 对于任意一个满足 1≤x<p 的正整数 x,均有 Cpx≡0(modp),证明如下:
Cpx≡x!(p−x)!p!(modp)≡p⋅x!(p−x)!(p−1)!(modp)≡0(modp)
【引理 2】 (a+b)p≡ap+bp(modp),证明如下:
根据二项式定理:
(a+b)p≡i=0∑pCpi⋅ai⋅bp−i(modp)≡ap+bp+i=1∑p−1Cpi⋅ai⋅bp−i(modp)≡ap+bp(modp)(由引理 1 得)
【定理 1(Lucas 定理)】 Cnm≡Cnmodpmmodp⋅C⌊pn⌋⌊pm⌋(modp),证明如下:
设 x 为一个满足 1≤x<p 的正整数,由二项式定理得:(1+x)n≡i=0∑nCnixi。
令 pn=⌊pn⌋,pm=⌊pm⌋,rn=nmodp,rm=mmodp,那么有:
(1+x)n≡(1+x)p⋅pn+rn(modp)≡[(1+x)p]pn⋅(1+x)rn(modp)≡(1+xp)pn⋅(1+x)rn(modp)≡i=0∑pnj=0∑rnCpni⋅Crnj⋅xip+j(modp)(根据二项式定理)≡i=0∑nCpn⌊pi⌋⋅Crnimodp⋅xi(modp)
所以:
i=0∑nCnixi≡i=0∑nCpn⌊pi⌋⋅Crnimodp⋅xi(modp)
Cnm≡Cpnpm⋅Crnrm(modp)
证毕。
练习:P3773 [CTSC2017]吉夫特
1x2 p 不是质数的情况(exLucas)
未完待续。