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

}