扩展卢卡斯定理学习笔记

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;
}

原文地址:https://www.cnblogs.com/Troverld/p/14620910.html