【BZOJ3129】[SDOI2013]方程（容斥，拓展卢卡斯定理）
题面
题解
因为答案是正整数，所先给每个位置都放一个就行了，然后\(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;
}