扩展Lucas定理
Fading
2019-07-05 10:52:08
### 扩展 Lucas 定理 exLucas
其他题解讲的都很不清楚，而且很奇怪。。。
所以我来发一篇良心题解，求大家点赞~
### 问题
求
$$C_n^m \bmod{p}$$
其中$p\leq 10^6$
### 算法
设
$$p=p_1^{\alpha_1}p_2^{\alpha_2}...$$
求出
$$\left\{\begin{aligned}C_n^m \bmod\ {p_1^{\alpha_1}} \\C_n^m \bmod\ {p_2^{\alpha_2}}\\ ......\end{aligned} \right.$$
最后用中国剩余定理合并即可。
假设我们现在要求
$$C_n^m\bmod{P^k}([P\in {\texttt{prime}}])$$
$$\therefore\ =\frac {n!}{m!(n-m)!}\bmod{P^k}$$
我们无法求得$m!,(n-m)!$的逆元，理由很简单，不一定有解。
怎么办呢？发现$a$对$p$有逆元的充要条件为
$$\gcd(a,p)=1$$
所以我们换一个式子：
$$=\frac {\frac {n!}{P^{x}}}{\frac {m!}{P^{y}}\frac {(n-m)!}{P^{z}}}P^{x-y-z}\bmod{P^k}$$
其中
$x$为$n!$中$P$因子的个数（包含多少$P$因子）。
其他都是同理的。
所以如果我们对于一个$n$，可以求出$\frac {n!}{P^x}$，那不就完事了吗？（这样我们就可以求逆元了）
问题等价于求
$$\frac {n!}{P^x}\bmod{P^k}$$
我们对$n!$进行变形：
$$n!=1\cdot 2\cdot 3\cdot ...\cdot n$$
$$=(P\cdot 2P\cdot 3P...)(1\cdot 2\cdot ...)$$
左边是$1\sim n$中所有$\leq n$的$P$的倍数。
右边反之。
因为$1\sim n$中有$\lfloor \frac nP\rfloor$个$P$的倍数，
$$\therefore\ =P^{\lfloor \frac nP\rfloor}(1 \cdot 2 \cdot 3...)(1\cdot 2\cdot ...)$$
$$=P^{\lfloor \frac nP\rfloor}(\lfloor \frac nP\rfloor)!\prod_{i=1,i\not\equiv 0\pmod {P}}^ni$$
显然后面这个$\bmod\ P$是有循环节的。
$$=P^{\lfloor \frac nP\rfloor}(\lfloor \frac nP\rfloor)!(\prod_{i=1,i\not\equiv 0\pmod {P}}^{P^k}i)^{\lfloor \frac n{P^k}\rfloor}(\prod_{i=P^k\lfloor \frac n{P^k}\rfloor,i\not\equiv 0\pmod {P}}^ni)$$
其中
$$(\prod_{i=1,i\not\equiv 0\pmod {P}}^{P^k}i)^{\lfloor \frac n{P^k}\rfloor}$$
是循环节$1\sim P$中所有无$P$因子数的乘积。
$$(\prod_{i=P^k\lfloor \frac n{P^k}\rfloor,i\not\equiv 0\pmod {P}}^ni)$$
是余项的乘积。
比如
$$22!\equiv(1\cdot 2\cdot 4\cdot 5\cdot 7\cdot 8)(10\cdot 11\cdot 13\cdot 14\cdot 16\cdot 17)(19\cdot 20\cdot 22)(3\cdot 6\cdot 9\cdot 12\cdot 15\cdot 18\cdot 21)\pmod {3^2}$$
$$\equiv(1\cdot 2\cdot 4\cdot 5\cdot 7\cdot 8)^2(19\cdot 20\cdot 22)3^7(1\cdot 2\cdot 3\cdot 4\cdot 5\cdot 6\cdot 7)\pmod {3^2}$$
$$\equiv3^77!(1\cdot 2\cdot 4\cdot 5\cdot 7\cdot 8)^2(19\cdot 20\cdot 22)\pmod {3^2}$$
发现正好是一样的。$\lfloor \frac {22}{3^2}\rfloor=2 $
理解了吧...
$$=P^{\lfloor \frac nP\rfloor}(\lfloor \frac nP\rfloor)!(\prod_{i=1,i\not\equiv 0\pmod {P}}^{P^k}i)^{\lfloor \frac n{P^k}\rfloor}(\prod_{i=P^k\lfloor \frac n{P^k}\rfloor,i\not\equiv 0\pmod {P}}^ni)$$
发现前面这个$P^{\lfloor \frac nP\rfloor}$是要除掉的。
然而$(\lfloor \frac nP\rfloor)!$里可能还包含$P$。
所以，我们定义函数：
$$f(n)=\frac {n!}{P^x}$$
$$\therefore\ f(n)=f(\lfloor \frac nP\rfloor)(\prod_{i=1,i\not\equiv 0\pmod {P}}^{P^k}i)^{\lfloor \frac n{P^k}\rfloor}(\prod_{i=P^k\lfloor \frac n{P^k}\rfloor,i\not\equiv 0\pmod {P}}^ni)$$
这样就好了。时间复杂度是$O(\log_{P}n)$。
边界$f(0)=1$
看看原式
$$=\frac {\frac {n!}{P^{x}}}{\frac {m!}{P^{y}}\frac {(n-m)!}{P^{z}}}P^{x-y-z}\bmod{P^k}$$
$$=\frac {f(n)}{f(m)f(n-m)}P^{x-y-z}\bmod{P^k}$$
由于$f(m)$一定与$P^k$互质，所以我们可以直接用$exgcd$求逆元啦。
最后我们还有一个问题，$P^{x-y-z}$咋求?
其实很好求。
比如说，我们要求$f(n)=\frac {n!}{P^x}$中的$x$
设$g(n)=x$(上述的$x$)
再看看阶乘的式子
$$n!=P^{\lfloor \frac nP\rfloor}(\lfloor \frac nP\rfloor)!(\prod_{i=1,i\not\equiv 0\pmod {P}}^{P^k}i)^{\lfloor \frac n{P^k}\rfloor}(\prod_{i=P^k\lfloor \frac n{P^k}\rfloor,i\not\equiv 0\pmod {P}}^ni)$$
很显然我们要的就是$P^{\lfloor \frac nP\rfloor}$
而且$(\lfloor \frac nP\rfloor)!$可能还有$P$因子。
所以我们可以得到以下递推式
$$g(n)=\lfloor \frac nP\rfloor+g(\lfloor \frac nP\rfloor)$$
时间复杂度是$O(\log_{P}n)$。
边界$g(n)=0(n<P)$
所以答案就是
$$=\frac {\frac {n!}{P^{x}}}{\frac {m!}{P^{y}}\frac {(n-m)!}{P^{z}}}P^{g(n)-g(m)-g(n-m)}\bmod{P^k}$$
最后用中国剩余定理合并答案即可得到$C_n^m$。
代码:
```cpp
#include<bits/stdc++.h>
#define ll long long
using namespace std;
#ifndef Fading
inline char gc(){
static char now[1<<16],*S,*T;
if (T==S){T=(S=now)+fread(now,1,1<<16,stdin);if (T==S) return EOF;}
return *S++;
}
#endif
#ifdef Fading
#define gc getchar
#endif
void exgcd(ll a,ll b,ll &x,ll &y){
if (!b) return (void)(x=1,y=0);
exgcd(b,a%b,x,y);
ll tmp=x;x=y;y=tmp-a/b*y;
}
ll gcd(ll a,ll b){
if (b==0) return a;
return gcd(b,a%b);
}
inline ll INV(ll a,ll p){
ll x,y;
exgcd(a,p,x,y);
return (x+p)%p;
}
inline ll lcm(ll a,ll b){
return a/gcd(a,b)*b;
}
inline ll mabs(ll x){
return (x>0?x:-x);
}
inline ll fast_mul(ll a,ll b,ll p){
ll t=0;a%=p;b%=p;
while (b){
if (b&1LL) t=(t+a)%p;
b>>=1LL;a=(a+a)%p;
}
return t;
}
inline ll fast_pow(ll a,ll b,ll p){
ll t=1;a%=p;
while (b){
if (b&1LL) t=(t*a)%p;
b>>=1LL;a=(a*a)%p;
}
return t;
}
inline ll read(){
ll x=0,f=1;char ch=gc();
while (!isdigit(ch)) {if (ch=='-') f=-1;ch=gc();}
while (isdigit(ch)) x=x*10+ch-'0',ch=gc();
return x*f;
}
inline ll F(ll n,ll P,ll PK){
if (n==0) return 1;
ll rou=1;//循环节
ll rem=1;//余项
for (ll i=1;i<=PK;i++){
if (i%P) rou=rou*i%PK;
}
rou=fast_pow(rou,n/PK,PK);
for (ll i=PK*(n/PK);i<=n;i++){
if (i%P) rem=rem*(i%PK)%PK;
}
return F(n/P,P,PK)*rou%PK*rem%PK;
}
inline ll G(ll n,ll P){
if (n<P) return 0;
return G(n/P,P)+(n/P);
}
inline ll C_PK(ll n,ll m,ll P,ll PK){
ll fz=F(n,P,PK),fm1=INV(F(m,P,PK),PK),fm2=INV(F(n-m,P,PK),PK);
ll mi=fast_pow(P,G(n,P)-G(m,P)-G(n-m,P),PK);
return fz*fm1%PK*fm2%PK*mi%PK;
}
ll A[1001],B[1001];
//x=B(mod A)
inline ll exLucas(ll n,ll m,ll P){
ll ljc=P,tot=0;
for (ll tmp=2;tmp*tmp<=P;tmp++){
if (!(ljc%tmp)){
ll PK=1;
while (!(ljc%tmp)){
PK*=tmp;ljc/=tmp;
}
A[++tot]=PK;B[tot]=C_PK(n,m,tmp,PK);
}
}
if (ljc!=1){
A[++tot]=ljc;B[tot]=C_PK(n,m,ljc,ljc);
}
ll ans=0;
for (ll i=1;i<=tot;i++){
ll M=P/A[i],T=INV(M,A[i]);
ans=(ans+B[i]*M%P*T%P)%P;
}
return ans;
}
signed main(){
ll n=read(),m=read(),P=read();
printf("%lld\n",exLucas(n,m,P));
return 0;
}
```