「SDOI2013」方程 - 计数+容斥原理+扩展Lucas定理

题目大意

    给定方程

$$X_1+X_2+\cdots+X_n=M$$

    我们对第$1\cdots N_1$个变量进行一些限制：

$$X_i\le A_i$$

    我们对第$N_1+1\cdots N_1+N_2$个变量进行一些限制：

$$X_i\ge A_i$$

求：在满足这些限制的前提下，该方程正整数解的个数。
答案可能很大，请输出对$p$取模后的答案，也即答案除以$p$的余数。

---

题目分析

考虑若没有限制，通过隔板法不难分析出答案为$C_{m-1}^{n-1}$（将$m$个物品间隙中选出$n-1$个位置放隔板）
限制考虑题目的两个限制，发现第二个大于等于的限制很好搞，只需要先事先放上$A_i-1$个，答案即为$C_{m-1-\sum_{i=N_1+1}^{N_1+N_2}(A_i-1)}^{n-1}$。

考虑如何处理第一个限制，其限制取反即为大于等于的限制，故对其进行容斥，考虑有多少个违背限制。

题目很简单，难点在于，$p$不一定是素数，且经过检验，$p$分解后有次数不为$1$的项，因此必须使用扩展Lucas定理，注意记忆化。

---

代码

#include<algorithm>
#include<iostream>
#include<iomanip>
#include<cstring>
#include<cstdlib>
#include<climits>
#include<vector>
#include<cstdio>
#include<cmath>
#include<queue>
#include<map>
using namespace std;

typedef long long LL;

inline const LL Get_Int() {
	LL num=0,bj=1;
	char x=getchar();
	while(x<'0'||x>'9') {
		if(x=='-')bj=-1;
		x=getchar();
	}
	while(x>='0'&&x<='9') {
		num=num*10+x-'0';
		x=getchar();
	}
	return num*bj;
}

void Exgcd(LL a,LL b,LL &gcd,LL &x,LL &y) {
	if(!b)gcd=a,x=1,y=0;
	else Exgcd(b,a%b,gcd,y,x),y-=x*(a/b);
}
 
LL inv(LL a,LL mod) {
	LL gcd,x,y;
	Exgcd(a,mod,gcd,x,y);
	return (x+mod)%mod;
}
 
LL Quick_Pow(LL a,LL b,LL mod) {
	LL ans=1;
	for(; b; b>>=1,a=a*a%mod)if(b&1)ans=ans*a%mod;
	return ans;
}

int cnt=0;
LL fac[25][10505];
map<LL,int> M;

LL Fac(LL n,LL p,LL num) {
	if(n==0)return 1;
	int id=M[num];
	if(fac[id][0]==-1) {
		fac[id][0]=fac[id][1]=1;
		for(LL i=2; i<=num; i++)
			if(i%p)fac[id][i]=fac[id][i-1]*i%num;
			else fac[id][i]=fac[id][i-1];
	}
	LL sum=Quick_Pow(fac[id][num],n/num,num)*fac[id][n%num];
	return sum*Fac(n/p,p,num)%num;
}
 
LL Cal(LL n,LL p) {
	LL sum=0;
	for(LL i=n; i; i/=p)sum+=i/p;
	return sum;
}
 
LL C(LL n,LL m,LL p,LL num) {
	if(n<m)return 0;
	if(!M[num])M[num]=++cnt;
	LL x=Fac(n,p,num),y=Fac(m,p,num),z=Fac(n-m,p,num);
	LL sum=Cal(n,p)-Cal(m,p)-Cal(n-m,p);
	return x*inv(y,num)%num*inv(z,num)%num*Quick_Pow(p,sum,num)%num;
}

vector<LL> A,B;

LL ExLucas(LL n,LL m,LL mod) {
	LL x=mod,ans=0;
	if(A.empty()) {
		for(LL i=2; i<=sqrt(mod); i++)
			if(x%i==0) {
				LL num=1;
				while(x%i==0)x/=i,num*=i;
				A.push_back(i);
				B.push_back(num);
				LL Mi=mod/num;
				ans=(ans+Mi*inv(Mi,num)%mod*C(n,m,i,num)%mod)%mod;
			}
		if(x>1) {
			A.push_back(x);
			B.push_back(x);
			ans=(ans+mod/x*inv(mod/x,x)%mod*C(n,m,x,x)%mod)%mod;
		}
	} else {
		for(int i=0; i<A.size(); i++) {
			LL x=A[i],num=B[i];
			LL Mi=mod/num;
			ans=(ans+Mi*inv(Mi,num)%mod*C(n,m,x,num)%mod)%mod;
		}
	}
	return ans;
}

LL n,m,n1,n2,p,a[25];

#define C(n,m) ExLucas(n,m,p)

int main() {
	memset(fac,-1,sizeof(fac));
	int t=Get_Int();
	p=Get_Int();
	while(t--) {
		n=Get_Int();
		n1=Get_Int();
		n2=Get_Int();
		m=Get_Int();
		for(int i=1; i<=n1+n2; i++)a[i]=Get_Int();
		for(int i=n1+1; i<=n1+n2; i++)m-=a[i]-1;
		LL ans=C(m-1,n-1);
		for(int S=1; S<(1<<n1); S++) {
			LL cnt=0,tmpm=m;
			for(int i=0; i<n1; i++)if(S&(1<<i))cnt++,tmpm-=a[i+1];
			cnt=(cnt&1)?-1:1;
			ans=(ans+cnt*C(tmpm-1,n-1)+p)%p;
		}
		printf("%lld\n",ans);
	}
	return 0;
}