题解：P4720 【模板】扩展卢卡斯定理/exLucas
TH911
·
2025-07-09 22:29:48
·
题解
题目传送门
算法介绍
exLucas 算法可以求 \dbinom{n}{m}\bmod P ，P 不一定为质数 。（但是，exLucas 实际上和 Lucas 没有多大关系……）
由唯一分解定理，可以对 P 进行质因数分解：
P=p_1^{c_1}p_2^{c_2}p_3^{c_3}\cdots p_n^{c^n}
算法流程
CRT 求解
我们如果能够求出 a_1,a_2,\cdots,a_n 使得：
\begin{cases}
\dbinom nm&\equiv a_1\pmod{p_1^{c_1}}\\
\dbinom nm&\equiv a_2\pmod{p_1^{c_2}}\\
&\cdots\\
\dbinom nm&\equiv a_n\pmod{p_n^{c_n}}\\
\end{cases}
那么我们就可以通过 CRT 求解出 \dbinom nm\bmod P 。因为在 CRT 中，恰好也有 P=p_1^{c_1}p_2^{c_2}\cdots p_n^{c_n} 。
求解模质数幂下余数
展开定义式：
\dbinom nm=\frac{n!}{m!(n-m)!}
因此：
\dbinom nm\equiv\frac{n!}{m!(n-m)!}\pmod{p_i^{c_i}}
因为 p_i^{c_i},m!,(n-m)! 不一定互质，因此乘法逆元不一定存在 。
考虑先提出 分子和分母中所有的 p_i 次幂，随后便可以使用逆元求解。
记 x 的质因数分解中质数 p 的幂次为 \nu_p(x) ，剩余积为 (x)_p ，即：
x=p^{\nu_p(x)}\cdot(x)_p
则有：
\frac{n!}{m!(n-m)!}\equiv\frac{(n!)_{p_i}}{(m!)_{p_i}((n-m)!)_{p_i}}\cdot p_i^{\nu_{p_i}(n!)-\nu_{p_i}(m!)-\nu_{p_i}((n-m)!)}\pmod{p_i^{c_i}}
那么，现在只需要考虑对于 x,p ，如何高效地 求出 \nu_{p}(x!)\bmod p_i^{c_i},(x!)_p\bmod p_i^{c_i} 。
考虑：
x!=1\times2\times\cdots\times p\times\cdots\times2p\times\cdots\times\left\lfloor\dfrac{x}{p}\right\rfloor p\times\cdots\times x
容易发现，在 p,2p,3p,\cdots,\left\lfloor\dfrac xp \right\rfloor p 中 p 的个数为 \left\lfloor\dfrac xp\right\rfloor+\nu_p\left(\left\lfloor\dfrac xp \right\rfloor!\right) 。其中 \nu_p\left(\left\lfloor\dfrac xp \right\rfloor!\right) 是 1,2,3,\cdots,\left\lfloor\dfrac xp\right\rfloor 中的 p 的个数。
将递推式展开：
$$
\begin{aligned}
\nu_p\left(x!\right)&=\left\lfloor\dfrac xp\right\rfloor+\nu_p\left(\left\lfloor\dfrac xp \right\rfloor!\right)\\
&=\left\lfloor\dfrac xp\right\rfloor+\left\lfloor\dfrac x{p^2}\right\rfloor+\nu_p\left(\left\lfloor\dfrac xp^2 \right\rfloor!\right)\\
&=\left\lfloor\dfrac xp\right\rfloor+\left\lfloor\dfrac x{p^2}\right\rfloor+\left\lfloor\dfrac x{p^3}\right\rfloor+\cdots
\end{aligned}
$$
因此可以 $\mathcal O(\log x)$ 计算 $\nu_p(x!)$：
```cpp
int v(ll n,ll p){
int ans=0;
do{
n/=p;
ans+=n;
}while(n);
return ans;//这里取不取模其实都无所谓,因为幂次不会太多,想取也可以.
}
```
***
现在考虑如何计算 $(x!)_p\bmod p^c$；显然不能利用定义式 $(x!)_p=\dfrac{x!}{\nu_p(x!)}$，而需要其他方法（因为无法得知 $x!$）。
不难进行递推：
$$
\begin{aligned}
(x!)_p&\equiv\prod_{i=1}^n(i)_p\\
&\equiv\left(\prod_{1\leq i\leq x,i\perp p}(i)_p\right)\left(\prod_{i=1}^{\left\lfloor\frac xp\right\rfloor}(p\cdot i)_p\right)\\
&\equiv\left(\prod_{1\leq i\leq x,i\perp p}(i)_p\right)\left(\left\lfloor\frac xp\right\rfloor!\right)_p\\
&\equiv\left(\prod_{1\leq i\leq p^c,i\perp p}i\right)^{\left\lfloor\frac{x}{p^c}\right\rfloor}\left(\prod_{1\leq i\leq x\bmod p^c,i\perp p}i\right)\left(\left\lfloor\frac xp\right\rfloor!\right)_p\\
\end{aligned}
$$
由 Wilson 定理的推论（$m=2,4,p^c,2p^c$ 时 $k\equiv-1$，否则为 $1$）：
$$
\prod_{1\leq k\leq m,k\perp m}k\equiv\pm1\pmod m
$$
因此：
$$
\begin{aligned}
(x!)_p&\equiv(\pm1)^{\left\lfloor\frac{x}{p^c}\right\rfloor}\left(\prod_{1\leq i\leq x\bmod p^c,i\perp p}i\right)\left(\left\lfloor\frac xp\right\rfloor!\right)_p
\end{aligned}
$$
可以发现，每次计算 $(x!)_p$ 时，$p^c$ 是固定的，因此可以预处理：
$$
f_j\equiv\prod_{1\leq i\leq j,i\perp p}i\pmod p^c
$$
因此，得到最终推导式：
$$
\begin{aligned}
(x!)_p&\equiv(\pm1)^{\left\lfloor\frac{x}{p^c}\right\rfloor}f_{x\bmod p^c}\left(\left\lfloor\frac xp\right\rfloor!\right)_p
\end{aligned}
$$
可以递归或迭代处理，时间复杂度 $\mathcal O(p^c+\log n x)$。
```cpp
int Wilson(int n,int p,int pc){
int ans=1;
vector<int>f(pc)
f[0]=1;
for(int i=1;i<pc;i++){
f[i]=i%p?f[i-1]*i%pc:f[i-1];
}
bool flag=p!=2||pc<=4;
while(n>1){
if((n/pc)&flag){
ans=pc-ans;
}
ans=ans*f[n%pc]%pc;
n/=p;
}
return ans;
}
```
***
# AC 代码
**一定要开 `__int128`！会爆 `long long`！**
```cpp
//#include<bits/stdc++.h>
#include<algorithm>
#include<iostream>
#include<cstring>
#include<iomanip>
#include<cstdio>
#include<string>
#include<vector>
#include<cmath>
#include<ctime>
#include<deque>
#include<queue>
#include<stack>
#include<list>
using namespace std;
typedef long long ll;
#define ll __int128
#define int __int128
constexpr const int PMAX=1e6;
void exgcd(ll a,ll &x,ll b,ll &y){
if(!b){
x=1;
y=0;
return;
}
ll tmp;
exgcd(b,tmp,a%b,x);
y=tmp-a/b*x;
}
ll inverse(ll a,ll p){
ll x,y;
exgcd(a,x,p,y);
x%=p;
if(x<0){
x+=p;
}
return x;
}
int qpow(int base,int n){
int ans=1;
while(n){
if(n&1){
ans=1ll*ans*base;
}
base=1ll*base*base;
n>>=1;
}
return ans;
}
int CRT(vector<int>a,vector<int>p){
ll L=1;
for(int &i:p){
L*=i;
}
ll x=0;
for(int i=0;i<a.size();i++){
ll q=L/p[i];
x=(x+1ll*a[i]*q%L*inverse(q,p[i])%L)%L;
}
if(x<0){
x+=L;
}
return x;
}
int v(ll n,ll p){
int ans=0;
do{
n/=p;
ans+=n;
}while(n);
return ans;
}
int Wilson(int n,int p,int pc){
int ans=1;
vector<int>f(pc);
f[0]=1;
for(int i=1;i<pc;i++){
f[i]=i%p?f[i-1]*i%pc:f[i-1];
}
bool flag=p!=2||pc<=4;
while(n>1){
if((n/pc)&flag){
ans=pc-ans;
}
ans=ans*f[n%pc]%pc;
n/=p;
}
return ans;
}
void breakDown(int P,vector<int>&p,vector<int>&pc){
int pp=P;
for(int i=2;i*i<=pp;i++){
if(pp%i==0){
p.push_back(i);
pc.push_back(1);
while(pp%i==0){
pc.back()*=i;
pp/=i;
}
}
}
if(pp>1){
p.push_back(pp);
pc.push_back(pp);
}
}
long long exLucas(ll n,ll m,ll P){
if(n<m){
return 0;
}
vector<int>p,pc;
breakDown(P,p,pc);
vector<int>a(pc.size());
for(int i=0;i<p.size();i++){
int nWilson=Wilson(n,p[i],pc[i]);
int mWilson=Wilson(m,p[i],pc[i]);
int nmWilson=Wilson(n-m,p[i],pc[i]);
a[i]=nWilson*inverse(mWilson*nmWilson%pc[i],pc[i])*qpow(p[i],v(n,p[i])-v(m,p[i])-v(n-m,p[i]));
}
return CRT(a,pc);
}
main(){
/*freopen("test.in","r",stdin);
freopen("test.out","w",stdout);*/
ios::sync_with_stdio(false);
cin.tie(0);cout.tie(0);
long long n,m,P;
cin>>n>>m>>P;
cout<<exLucas(n,m,P)<<'\n';
cout.flush();
/*fclose(stdin);
fclose(stdout);*/
return 0;
}
/*
1000000000000000000 500000000000000000 998243
0
*/
```