题链：
http://www.lydsy.com/JudgeOnline/problem.php?id=3129
题解：
容斥，扩展Lucas，中国剩余定理
先看看不管限制，只需要每个位置都是正整数时的方案数的求法。
假设有 N 个未知数，加起来的和为 M。
转化一下问题变为："小球分配"
有 M 个相同的小球，放在 N 个盒子里，且每个盒子至少有一个的方案数。
那么方案数为 ${C}_{M-1}^{N-1}$
怎么理解这个式子呢？"插隔板法"。
使 M个小球放在一排，考虑可以在相邻小球的空隙间(共 M-1个空隙)插入一个隔板，总共插入 N-1个隔板。
把相邻的两个隔板(把首尾也看作另外2个隔板)中间的区域看做一个个的盒子，区域内的小球即为该盒子里的小球。
则任意一种插隔板的方法都对应一种把小球放入盒子的方法。
所以方案数为 ${C}_{M-1}^{N-1}$
对于题目给的两种限制，不要被吓到了。
其实大于等于(>=W)这一个限制很好处理，直接把 M-=(W-1)，即事先给这些位置分配 W-1个小球。
因为接下来按照上面的 "小球分配" 方式，这些限制一定可以满足。
然后对于小于等于(<=W) 就需要用到容斥了。
定义 f[S] 表示至少 S 集合里面的这些"小于等于"限制都不能满足，那么方案数的求法就是：
先强制给这些盒子分配 Wi 个小球，
使得接下来按照上面的 "小球分配" 方式求出方案数，这些盒子一定会超过限制。
所以按照容斥的做法:
答案 ANS = 至少 0 个盒子超出限制的方案数
- 至少 1 个盒子超出限制的方案数
+至少 2 个盒子超出限制的方案数
-...+...
但是求那个方案数 C(N,M) 就比较麻烦了，因为 N,M很大，且模数不保证为质数，
所以用到 扩展Lucas(+国剩) 来求。
代码：
#include<cstdio>
#include<cstring>
#include<iostream>
#define _ % P
#define filein(x) freopen(#x".in","r",stdin);
#define fileout(x) freopen(#x".out","w",stdout);
using namespace std;
int h[300],k[300],d[300],pi[10],ppi[10],fc[10][100500],pnt;
int T,N,N1,N2,M,ANS;
void init(){
ANS=0;
memset(h,0,sizeof(h));
memset(k,0,sizeof(k));
memset(d,0,sizeof(d));
}
void pre_prime(int P){
static bool np[100500];
static int prime[100500],cnt;
for(int i=2;i<=100000;i++){
if(!np[i]){
prime[++cnt]=i;
if(P%i==0){
pi[++pnt]=i; ppi[pnt]=1;
while(P%i==0) ppi[pnt]*=i,P/=i;
}
}
for(int j=1;j<=cnt&&i<=100000/prime[j];j++){
np[i*prime[j]]=1;
if(i%prime[j]==0) break;
}
}
for(int i=1;i<=pnt;i++){
fc[i][0]=1;
for(int j=1;j<=100000;j++)
fc[i][j]=(1ll*fc[i][j-1]*(j%pi[i]?j:1))%ppi[i];
}
}
int pow(int a,int b,int P){
int ret=1; a=(a)_;
while(b){
if(b&1) ret=(1ll*ret*a)_;
a=(1ll*a*a)_; b>>=1;
}
return ret;
}
void exEuclidean(int a,int b,int &g,long long &x,long long &y){
if(!b) g=a,x=1,y=0;
else exEuclidean(b,a%b,g,y,x),y-=(a/b)*x;
}
int inv(int a,int P){
int g; long long x,y;
exEuclidean(a,P,g,x,y);
x=((1ll*x)_+P)_;
return (int)x;
}
int fac(int n,int I,int Pi,int P){
if(!n) return 1; int ret=1,r=n%P;
//for(int i=2;i<=P;i++) if(i%Pi) ret=(1ll*ret*i)_;
ret=fc[I][P];
ret=pow(ret,n/P,P);
//for(int i=2;i<=r;i++) if(i%Pi) ret=(1ll*ret*i)_;
ret=(1ll*ret*fc[I][r])_;
return (1ll*ret*fac(n/Pi,I,Pi,P))_;
}
int C(int n,int m,int I,int Pi,int P){
int x,y,z,c=0;
x=fac(n,I,Pi,P); y=fac(m,I,Pi,P); z=fac(n-m,I,Pi,P);
for(int i=n;i;i/=Pi) c+=i/Pi;
for(int i=m;i;i/=Pi) c-=i/Pi;
for(int i=n-m;i;i/=Pi) c-=i/Pi;
return ((((1ll*x*inv(y,P))_)*inv(z,P))_*pow(Pi,c,P))_;
}
int exLucas(int n,int m,int P){
int ret=0,tmp; if(n<m) return 0;
for(int i=1;i<=pnt;i++){
tmp=C(n,m,i,pi[i],ppi[i]);
tmp=((1ll*tmp*(P/ppi[i]))_*inv(P/ppi[i],ppi[i]))_;
ret=(1ll*ret+tmp)_;
}
return ret;
}
int main()
{
int P;
scanf("%d%d",&T,&P);
pre_prime(P);
while(T--){
init();
scanf("%d%d%d%d",&N,&N1,&N2,&M);
for(int i=1;i<=N1;i++) scanf("%d",&d[1<<(i-1)]);
for(int i=1,x;i<=N2;i++) scanf("%d",&x),M-=(x-1);
for(int s=1,p;p=s&-s,s<(1<<N1);s++)
h[s]=h[s^p]+1,k[s]=k[s^p]+d[p];
for(int s=0,m,tmp;s<(1<<N1);s++){
m=M-k[s];
tmp=exLucas(m-1,N-1,P);
if(h[s]&1) tmp=(-1ll*tmp+P)_;
ANS=((1ll*ANS+tmp)_+P)_;
}
printf("%d
",ANS);
}
return 0;
}