扩展卢卡斯定理学习笔记

IV.exLucas（扩展卢卡斯定理）

虽然是这个名字，但是它跟常规卢卡斯没有半毛钱关系

exLucas也是用来计算 $\dbinom nm\bmod p$ 的。不同于普通Lucas，这里的 $p$ 可以不为质数。

对于不为质数的模数，一个常规的想法是对其分解质因数，然后考虑其对于每个质数的次幂 $b_i^{k_i}$ 的模数，最后用CRT将所有东西怼一块即可。

于是我们现在考虑 $\dbinom nm\bmod b^k$。首先，一个显然的想法是，发现仅有 $b$ 的倍数在模 $b^k$ 意义下无逆元，于是我们求出 $n!,m!,(n-m)!$ 这三个数中有 $b$ 的多少次幂，$b$ 的幂次部分可以直接相减得到最终组合数中有多少个 $b$ 的次幂，而剩余部分可以直接求逆元处理。

于是问题变成两个：如何求一个阶乘中有多少个 $b$，以及除去 $b$ 后剩余了什么。

首先，第一个问题是常规问题，令 $f(x)$ 表示 $x!$ 中有多少个 $b$，则 $1\sim x$ 中所有 $b$ 的倍数显然是 $b$ 的唯一来源，对这些 $b$ 的倍数各除上一个 $b$ 后显然又会出现新的一组阶乘，故如此递归下去即可。边界为 $f(0)=0$。递推式为 $f(x)=f\left(\left\lfloor\dfrac xb\right\rfloor\right)+\left\lfloor\dfrac xb\right\rfloor$。

然后，同理，第二个问题也可以分成两部分：$b$ 的倍数以及非 $b$ 的倍数。$b$ 的倍数就除去 $b$ 继续递归即可，非 $b$ 的倍数，显然其具有 $b^k$ 的循环节，于是预处理出 $0\sim b^k$ 中所有数的阶乘（阶乘中要除去 $b$ 的倍数），则单个循环节内所有东西的积明显是 $(b^k)!$（当然，除去 $b$ 的倍数），循环节个数是 $\left\lfloor\dfrac x{b^k}\right\rfloor$，于是做一个快速幂即可。至于剩余的未完结的循环节，就直接查表就行了，反正我们已经预处理出来了。

时间复杂度 $O\Big(p+\log n\log p\Big)$，因为对于每一个质因数 $b$ 都要 $O(b^k)$ 地预处理阶乘，这部分复杂度共 $O(p)$；质因数数量是 $\log$ 级别的，然后对于每个质因数都要递归地处理 $n$，明显递归层数为 $\log n$。

IV.I.【模板】扩展卢卡斯

代码：

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
ll n,m;
int P,M,tot,p[30],a[30],b[30],k[30],fac[1001000],phi[30],res;
ll calcfacb(ll ip,int id){
	if(!ip)return 0;
	return calcfacb(ip/b[id],id)+ip/b[id];
}
int ksm(int x,ll y,int mod){
	int z=1;
	for(;y;y>>=1,x=1ll*x*x%mod)if(y&1)z=1ll*z*x%mod;
	return z;
}
int calcfacm(ll ip,int id){
	if(!ip)return 1;
	return 1ll*ksm(fac[p[id]],ip/p[id],p[id])*fac[ip%p[id]]%p[id]*calcfacm(ip/b[id],id)%p[id];
}
int C(ll x,ll y,int id){
	int A=calcfacm(x,id);
	A=1ll*A*ksm(calcfacm(y,id),phi[id]-1,p[id])%p[id]*ksm(calcfacm(x-y,id),phi[id]-1,p[id])%p[id];
	ll B=calcfacb(x,id)-calcfacb(y,id)-calcfacb(x-y,id);
	return 1ll*A*ksm(b[id],B,p[id])%p[id];
}
void newprime(int i){
	b[++tot]=i,p[tot]=1;while(!(P%i))k[tot]++,P/=i,p[tot]*=i;
	phi[tot]=p[tot]/i*(i-1);
	fac[0]=1;for(int j=1;j<=p[tot];j++)if(j%i)fac[j]=1ll*fac[j-1]*j%p[tot];else fac[j]=fac[j-1];
	a[tot]=C(n,m,tot);
}
int main(){
	scanf("%lld%lld%d",&n,&m,&P),M=P;
	for(int i=2;i*i<=P;i++)if(!(P%i))newprime(i);
	if(P!=1)newprime(P);
	for(int i=1;i<=tot;i++)(res+=1ll*a[i]*ksm(M/p[i],phi[i]-1,p[i])%M*(M/p[i])%M)%=M;
	printf("%d\n",res);
	return 0;
}

IV.II.[国家集训队]礼物

其实也是模板，但是又写了一道练手。

感觉这次的代码比上一题要清爽很多？

代码：

#include<bits/stdc++.h>
using namespace std;
int p,n,m,a[10];
int ksm(int x,int y,int mod){
	int z=1;
	for(;y;y>>=1,x=1ll*x*x%mod)if(y&1)z=1ll*z*x%mod;
	return z;
}
int B,K,P,phi,fac[100100];
int M,A;
int calcB(int x){int ret=0;while(x)ret+=(x=x/B);return ret;}
int calcfac(int x){if(!x)return 1;return 1ll*ksm(fac[P-1],x/P,P)*fac[x%P]%P*calcfac(x/B)%P;}
int C(int x,int y){return 1ll*calcfac(x)*ksm(calcfac(y),phi-1,P)%P*ksm(calcfac(x-y),phi-1,P)%P*ksm(B,calcB(x)-calcB(y)-calcB(x-y),P)%P;}
void func(int ip){
	B=ip,P=1,K=0;while(!(p%B))p/=B,K++,P*=B;
	phi=P/B*(B-1);
	fac[0]=1;for(int i=1;i<P;i++)if(i%B)fac[i]=1ll*fac[i-1]*i%P;else fac[i]=fac[i-1];
	int S=1;for(int i=1,j=n;i<=m;i++)S=1ll*S*C(j,a[i])%P,j-=a[i];
//	printf("%d %d %d %d:%d\n",B,K,P,phi,S);
	(A+=1ll*S*(M/P)%M*ksm((M/P),phi-1,P)%M)%=M;
}
int main(){
	scanf("%d%d%d",&p,&n,&m),M=p;for(int i=1;i<=m;i++)scanf("%d",&a[i]);
	for(int i=1,j=n;i<=m;i++){j-=a[i];if(j<0){puts("Impossible");return 0;}}
	for(int i=2;i*i<=p;i++)if(!(p%i))func(i);
	if(p!=1)func(p);
	printf("%d\n",A);
	return 0;
}

IV.III.[SDOI2013]方程

稍微不那么板子一点了。

首先，对于 $X\geq A$ 一类的限制，明显可以提前将 $A-1$ 配给 $x$ 使得其限制同其它数一样，都是 $\geq1$。

然后，对于 $X\leq A$ 一类限制，反正 $n_1$ 很小，直接容斥就行了。

于是问题转换为 $n$ 个正整数和为 $m$ 的方案数。依照小学数学·隔板法，其方案数为 $\dbinom{m-1}{n-1}$。用exLucas随便算即可。

时间复杂度 $O(2^{n_1}\times\log P\times\log^2n)$。

代码：

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
int T,p,n,n1,n2,a[10];
ll m; 
int ksm(int x,int y,int mod){
	int z=1;
	for(;y;y>>=1,x=1ll*x*x%mod)if(y&1)z=1ll*z*x%mod;
	return z;
}
int B,K,P,phi,fac[100100];
int M,A;
int calcB(int x){int ret=0;while(x)ret+=(x=x/B);return ret;}
int calcfac(int x){if(!x)return 1;return 1ll*ksm(fac[P-1],x/P,P)*fac[x%P]%P*calcfac(x/B)%P;}
int C(int x,int y){if(x<0||x<y)return 0;return 1ll*calcfac(x)*ksm(calcfac(y),phi-1,P)%P*ksm(calcfac(x-y),phi-1,P)%P*ksm(B,calcB(x)-calcB(y)-calcB(x-y),P)%P;}
void func(int ip){
	B=ip,P=1,K=0;while(!(p%B))p/=B,K++,P*=B;
	phi=P/B*(B-1);
	fac[0]=1;for(int i=1;i<P;i++)if(i%B)fac[i]=1ll*fac[i-1]*i%P;else fac[i]=fac[i-1];
	int S=0;for(int i=0;i<(1<<n1);i++){
		ll newm=m;
		for(int j=0;j<n1;j++)if(i&(1<<j))newm-=a[j];
		if(newm<0)continue;
		if(__builtin_popcount(i)&1)(S+=P-C(newm-1,n-1))%=P;else (S+=C(newm-1,n-1))%=P;
	}
	(A+=1ll*S*(M/P)%M*ksm((M/P),phi-1,P)%M)%=M;
}
int main(){
	scanf("%d%d",&T,&M);
	while(T--){
		p=M,scanf("%d%d%d%lld",&n,&n1,&n2,&m),A=0;
		for(int i=0;i<n1;i++)scanf("%d",&a[i]);
		for(int j=0,x;j<n2;j++)scanf("%d",&x),m-=(x-1);
		if(m<0){puts("0");continue;}
		for(int i=2;i*i<=p;i++)if(!(p%i))func(i);
		if(p!=1)func(p);
		printf("%d\n",A);
	}
	return 0;
}