[SDOI2013]方程
Description
给定方程
\[X_1+X_2+...+X_n=M
\]
我们对第\(1...n_1\)个变量进行一些限制：
\(X_1\leqslant A_1\)
\(X_2\leqslant A_2\)
\(X_{n_1}\leqslant A_{n_1}\)
我们对第\(n_1+1...n_1+n_2\)个变量进行一些限制：
\(X_{n_1+1}\geqslant A_{n_1+1}\)
\(X_{n_1+2}\geqslant A_{n_1+2}\)
\(X_{n_1+n_2}\geqslant A_{n_1+n_2}\)
求：在满足这些限制的前提下，该方程正整数解的个数。
答案可能很大，请输出对p取模后的答案，也即答案除以p的余数。
Input
输入含有多组数据，第一行两个正整数T，p。T表示这个测试点内的数据组数，p的含义见题目描述。
对于每组数据，第一行四个非负整数n，n1，n2，m。
第二行nl+n2个正整数，表示A1..n1+n2。请注意，如果n1+n2等于0，那么这一行会成为一个空行。
Output
共T行，每行一个正整数表示取模后的答案。
Sample Input
3 10007
3 1 1 6
3 3
3 0 0 5
3 1 1 3
3 3
Sample Output
3
6
0
如果不考虑前\(n_1\)个变量的限制，那么答案为
\[\binom{M-\sum\limits_{i=n_1+1}^{n_1+n_2}A_i-1}{n-1}
\]
现在多了\(n_1\)个限制，但\(n_1\)很小，我们可以暴力容斥，枚举哪些变量不满足限制即可
/*program from Wolfycz*/
#include<cmath>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#define inf 0x7f7f7f7f
#define Fi first
#define Se second
#define MK std::make_pair
typedef long long ll;
typedef unsigned int ui;
typedef unsigned long long ull;
typedef std::pair<int,int> pii;
inline char gc(){
static char buf[1000000],*p1=buf,*p2=buf;
return p1==p2&&(p2=(p1=buf)+fread(buf,1,1000000,stdin),p1==p2)?EOF:*p1++;
}
template<typename T>inline T frd(T x){
int f=1; char ch=gc();
for (;ch<'0'||ch>'9';ch=gc()) if (ch=='-') f=-1;
for (;ch>='0'&&ch<='9';ch=gc()) x=(x<<1)+(x<<3)+ch-'0';
return x*f;
}
template<typename T>inline T read(T x){
int f=1;char ch=getchar();
for (;ch<'0'||ch>'9';ch=getchar()) if (ch=='-') f=-1;
for (;ch>='0'&&ch<='9';ch=getchar()) x=(x<<1)+(x<<3)+ch-'0';
return x*f;
}
inline void print(int x){
if (x<0) putchar('-'),x=-x;
if (x>9) print(x/10);
putchar(x%10+'0');
}
template<typename T>inline T min(T x,T y){return x<y?x:y;}
template<typename T>inline T max(T x,T y){return x>y?x:y;}
template<typename T>inline T swap(T &x,T &y){T t=x; x=y,y=t;}
const int N=1e4;
namespace Math{
int P[10],V[10],C[10],f[10][N+10],tot,SP;
int mlt(int a,int b,int p){
int res=1;
for (;b;b>>=1,a=1ll*a*a%p) if (b&1) res=1ll*res*a%p;
return res;
}
void prepare(int p){
SP=p;
for (int i=2;i*i<=p;i++){
if (p%i) continue;
P[++tot]=i,V[tot]=1;
while (p%i==0) p/=i,V[tot]*=i,C[tot]++;
f[tot][0]=1;
for (int j=1;j<=V[tot];j++) f[tot][j]=1ll*f[tot][j-1]*(j%i?j:1)%V[tot];
}
if (p!=1){
f[++tot][0]=1;
P[tot]=V[tot]=p,C[tot]=1;
for (int j=1;j<=V[tot];j++) f[tot][j]=1ll*f[tot][j-1]*(j%p?j:1)%V[tot];
}
}
int gcd(int a,int b){return !b?a:gcd(b,a%b);}
void exgcd(int a,int b,int &x,int &y){
if (!b){x=1,y=0;return;}
exgcd(b,a%b,x,y);
int t=x; x=y,y=t-a/b*y;
}
int Ex_GCD(int a,int b,int c){
int d=gcd(a,b),x,y;
if (c%d) return -1;
a/=d,b/=d,c/=d;
exgcd(a,b,x,y);
x=(1ll*x*c%b+b)%b;
return x;
}
pii work(int n,int i){
if (n<=1) return MK(1,0);
int res=1ll*mlt(f[i][V[i]],n/V[i],V[i])*f[i][n%V[i]]%V[i];
pii tmp=work(n/P[i],i); int cnt=n/P[i];
return MK(1ll*res*tmp.Fi%V[i],tmp.Se+cnt);
}
int calc(int n,int m,int i){
int res=1,cnt=0; pii tmp;
tmp=work(n,i); cnt+=tmp.Se;
res=1ll*res*tmp.Fi;
tmp=work(m,i); cnt-=tmp.Se;
res=1ll*res*Ex_GCD(tmp.Fi,V[i],1)%V[i];
tmp=work(n-m,i); cnt-=tmp.Se;
res=1ll*res*Ex_GCD(tmp.Fi,V[i],1)%V[i];
return cnt<C[i]?1ll*res*mlt(P[i],cnt,V[i])%V[i]:0;
}
int Ex_C(int n,int m){
if (n<m) return 0;
if (tot==1) return Ex_GCD(1,V[1],calc(n,m,1));
int Ans=0;
for (int i=1;i<=tot;i++) Ans=(1ll*Ex_GCD(SP/V[i],V[i],1)*(SP/V[i])%SP*calc(n,m,i)%SP+Ans)%SP;
return Ans;
}
}
int main(){
int T=read(0),p=read(0);
Math::prepare(p);
for (;T;T--){
static int A[10];
int n=read(0),x=read(0),y=read(0),m=read(0),Ans=0;
for (int i=0;i<x;i++) A[i]=read(0);
for (int i=0;i<y;i++) m-=read(0)-1;
for (int sta=0;sta<1<<x;sta++){
int del=0,cnt=0;
for (int i=0;i<x;i++) if (sta>>i&1) cnt++,del+=A[i];
Ans=(Ans+(cnt&1?-1:1)*Math::Ex_C(m-del-1,n-1))%p;
}
printf("%d\n",(Ans+p)%p);
}
return 0;
}
