【BZOJ3129】[SDOI2013]方程（容斥，拓展卢卡斯定理）

【BZOJ3129】[SDOI2013]方程（容斥，拓展卢卡斯定理）

题面

BZOJ
洛谷

题解

因为答案是正整数，所先给每个位置都放一个就行了，然后\(A\)都要减一。
大于的限制和没有的区别不大，提前给他\(A_i\)个就好了。
假如没有小于的限制的话，那么就是经典的隔板法直接算答案。
如果提前给完之后，还剩\(M\)个球，要放进\(n\)个盒子，答案就是\(\displaystyle{M+n-1\choose n-1}\)
然而有一个小于的限制很烦人。发现数量很少，那么直接爆枚子集，强制一些不合法，然后容斥计算答案即可。
模数蛋疼无比，需要\(exLucas\)。
\(exLucas\)的时候因为模数相同，所以可以提前预处理模每个质因子情况下的阶乘，这样会快很多。

#include<iostream>
#include<cstdio>
using namespace std;
inline int read()
{
	int x=0;bool t=false;char ch=getchar();
	while((ch<'0'||ch>'9')&&ch!='-')ch=getchar();
	if(ch=='-')t=true,ch=getchar();
	while(ch<='9'&&ch>='0')x=x*10+ch-48,ch=getchar();
	return t?-x:x;
}
int fpow(int a,int b)
{
	int s=1;
	while(b){if(b&1)s*=a;a*=a;b>>=1;}
	return s;
}
int fpow(int a,int b,int MOD)
{
	int s=1;
	while(b){if(b&1)s=1ll*s*a%MOD;a=1ll*a*a%MOD;b>>=1;}
	return s;
}
void exgcd(int a,int b,int &x,int &y)
{
	if(!b){x=1;y=0;return;}
	exgcd(b,a%b,y,x);y-=a/b*x;
}
int inv(int n,int m)
{
	int x,y;exgcd(n,m,x,y);
	x=(x%m+m)%m;return x;
}
int fac[50],pw[50],tot=0;
int jc[50][11000];
int JC(int n,int p,int MOD,int id,int &z)
{
	if(!n){z=0;return 1;}
	int ret=JC(n/p,p,MOD,id,z);z+=n/p;
	if(n>MOD)ret=1ll*ret*fpow(jc[id][MOD],n/MOD,MOD)%MOD,n%=MOD;
	ret=1ll*ret*jc[id][n]%MOD;
	return ret;
}
int exLucas(int n,int m,int P)
{
	if(n<0||m<0||n<m)return 0;
	int M[50],V[50];
	for(int i=1;i<=tot;++i)
	{
		int MOD=fpow(fac[i],pw[i]),zero=0,a=1,b=1,z;
		a=JC(n,fac[i],MOD,i,z);zero+=z;
		b=JC(m,fac[i],MOD,i,z);zero-=z;
		b=1ll*b*JC(n-m,fac[i],MOD,i,z)%MOD;zero-=z;
		V[i]=1ll*a*inv(b,MOD)%MOD*fpow(fac[i],zero,MOD)%MOD;
		M[i]=MOD;
	}
	for(int i=2;i<=tot;++i)
	{
		int x,y;exgcd(M[1],M[i],x,y);
		x=(1ll*x*(V[i]-V[1])%M[i]+M[i])%M[i];
		V[1]=(1ll*M[1]*x+V[1])%(M[1]*M[i]);
		M[1]*=M[i];
	}
	return V[1];
}
int n,n1,n2,m;
int A[20];
int main()
{
	int T=read(),P=read(),MOD=P;
	for(int i=2;i*i<=P;++i)
		if(P%i==0)
		{
			fac[++tot]=i;pw[tot]=0;
			while(P%i==0)++pw[tot],P/=i;
		}
	if(P>1)++tot,fac[tot]=P,pw[tot]=1;
	for(int i=1;i<=tot;++i)
	{
		int MOD=fpow(fac[i],pw[i]);
		jc[i][0]=1;
		for(int j=1;j<=MOD;++j)jc[i][j]=1ll*jc[i][j-1]*((j%fac[i])?j:1)%MOD;
	}
	while(T--)
	{
		n=read();n1=read();n2=read();m=read()-n;
		for(int i=1;i<=n1+n2;++i)A[i]=read()-1;
		for(int i=1;i<=n2;++i)m-=A[n1+i];
		int S=1<<n1,ans=0;
		for(int i=0;i<S;++i)
		{
			int N=m,opt=1;
			for(int j=0;j<n1;++j)
				if(i&(1<<j))opt=MOD-opt,N-=A[j+1]+1;
			ans=(ans+1ll*opt*exLucas(N+n-1,n-1,MOD))%MOD;
		}
		printf("%d\n",ans);
	}
	return 0;
}