SDOI2013 方程

P3301 [SDOI2013] 方程

题目大意

给定方程$x_1+x_2+\cdots +x_n=m$。

对于$1\le i\le n_1$，需要满足$x_i\le a_i$。对于$n_1+1\le i\le n_1+n_2$，需要满足$x_i\ge a_i$。

求该方程的正整数解的个数。输出答案对$p$取模后的结果。

有$T$组数据。

$1\le T\le 5,1\le n,m\le 10^9,0\le n_1,n_2\le \min (8,m)$，$1\le p\le 437367875$

题解

前置知识：扩展lucas定理

对于$x_i\ge a_i$，可以将$m$减去$a_i-1$，使得$x_i$需要满足的条件变为$x_i\ge 1$，相当于没有限制。

对于$x_i\le a_i$，因为$n_1$比较小，所以可以容斥解决。不满足要求的条件为$x_i\ge a_i+1$，这样就转化成了上面的情况。

每种情况用插板法可得结果为$C_{m-1}^{n-1}$，因为$p$不一定是质数，所以要用扩展lucas定理。

还要加上一些优化，就是要预处理出模数的每个因子的幂。

code

#include<bits/stdc++.h>
using namespace std;
int T,n1,n2,tot,s1=0;
long long mod,x,y,b[25],a[105],r[105],wt[105];
long long f[105][20005];
void exgcd(long long c,long long d){
	if(d==0){
		x=1;y=0;
		return;
	}
	exgcd(d,c%d);
	long long t=x;x=y;y=t-c/d*y;
}
long long mi(long long t,long long v){
	if(!v) return 1;
	long long re=mi(t,v/2);
	re=re*re%mod;
	if(v&1) re=re*t%mod;
	return re;
}
long long gt(long long v,long long p,long long q,int wh){
	if(!v) return 1;
	long long re=mi(f[wh][q],v/q)%q;
	re=re*f[wh][v%q]%q*gt(v/p,p,q,wh)%q;
	return re;
}
long long C(long long v1,long long v2,long long p,long long q,int wh){
	if(v1<v2) return 0;
	long long f1=gt(v1,p,q,wh),f2=gt(v2,p,q,wh),f3=gt(v1-v2,p,q,wh),vt=0;
	for(long long i=p;i<=v1;i*=p) vt+=v1/i;
	for(long long i=p;i<=v2;i*=p) vt-=v2/i;
	for(long long i=p;i<=v1-v2;i*=p) vt-=(v1-v2)/i;
	return mi(p,vt)%q*f1%q*(mi(f2,q-q/p-1)%q)%q*(mi(f3,q-q/p-1)%q)%q;
}
long long exlucas(long long n,long long m){
	long long v=mod;
	tot=0;
	for(int i,j=1;j<=s1;j++){
		i=wt[j];
		r[++tot]=1;
		while(v%i==0){
			r[tot]*=i;
			v/=i;
		}
		a[tot]=C(n,m,i,r[tot],tot);
	}
	v=mod;
	long long re=0;
	for(int i=1;i<=tot;i++){
		exgcd(v/r[i],r[i]);
		x=(x%r[i]+r[i])%r[i];
		re=(re+v/r[i]*a[i]*x%v)%v;
	}
	return re;
}
long long solve(long long n,long long m){
	long long re=0;
	for(int i=0;i<1<<n1;i++){
		long long tmp=1,wm=m;
		for(int j=1;j<=n1;j++){
			if((i>>j-1)&1){
				tmp=mod-tmp;
				wm-=b[j];
			}
		}
		re=(re+tmp*exlucas(wm-1,n-1))%mod;
	}
	return re;
}
void init(){
	long long v=mod;
	for(int i=2;i*i<=v;i++){
		if(v%i==0){
			wt[++s1]=i;
			while(v%i==0) v/=i;
		}
	}
	if(v>1) wt[++s1]=v;
	v=mod;
	for(int i=1;i<=s1;i++){
		long long tmp=1;
		while(v%wt[i]==0){
			v/=wt[i];
			tmp*=wt[i];
		}
		f[i][0]=1;
		for(int j=1;j<=tmp;j++){
			f[i][j]=f[i][j-1];
			if(j%wt[i]) f[i][j]=f[i][j]*j%tmp;//预处理出模数的每个因子的幂
		}
	}
}
int main()
{
	long long n,m;
	scanf("%d%lld",&T,&mod);
	init();
	while(T--){
		scanf("%lld%d%d%lld",&n,&n1,&n2,&m);
		for(int i=1;i<=n1+n2;i++) scanf("%lld",&b[i]);
		for(int i=n1+1;i<=n1+n2;i++) m-=b[i]-1;
		printf("%lld\n",solve(n,m));
	}
	return 0;
}