SDOI2013 方程

传送门
看到这个题想起来熟悉的插板(正整数解，不能为空)。
首先看后一种限制，这个好办，强制先往这些盒子里放一些球，之后再插板。
前一种限制怎么办……？ 那我们考虑反着做，用容斥原理计算出不合法的情况，减掉就可以了。
想起来还是很简单的，但是这玩意实现起来真是恶心的要死…… 模数可以分解，于是我全部预处理了阶乘，不然是真的过不去……
看一下暴力的代码。

// luogu-judger-enable-o2
#include<bits/stdc++.h>
#define rep(i,a,n) for(ll i = a;i <= n;i++)
#define per(i,n,a) for(ll i = n;i >= a;i--)
#define enter putchar('
')

using namespace std;
typedef long long ll;
const int M = 200005;
const int INF = 2147483647;
const ll mod1 = 262203414;
const ll mod2 = 437367875;
const double eps = 1e-8;

ll read()
{
   ll ans = 0,op = 1;char ch = getchar();
   while(ch < '0' || ch > '9') {if(ch == '-') op = -1;ch = getchar();}
   while(ch >= '0' && ch <= '9') ans = ans * 10 + ch - '0',ch = getchar();
   return ans * op;
}

ll T,p,n,n1,n2,m,A[20],sum,ans,fac[M],inv[M];
ll bas[10] = {0,2,3,11,397,10007};

ll mul(ll a,ll b,ll t)
{
   ll res = a * b - (ll)((long double)a / t * b + eps) * t;
   return (res % t + t) % t;
}
ll exgcd(ll a,ll b,ll &x,ll &y)
{
   if(!b){x = 1,y = 0;return a;}
   ll d = exgcd(b,a%b,y,x);
   y -= a / b * x;
   return d;
}
ll Inv(ll a,ll b)
{
   ll x,y;
   exgcd(a,b,x,y);
   return (x % b + b) % b;
}
ll CRT(ll b,ll t) {return mul(mul(b,Inv(p/t,t),p),p/t,p);}
ll qpow(ll a,ll b,ll t)
{
   ll p = 1;
   while(b)
   {
      if(b & 1) p = mul(p,a,t);
      a = mul(a,a,t),b >>= 1;
   }
   return p;
}
void init(ll p,ll pk)
{
   fac[0] = 1;
   rep(i,1,pk) (i%p) ? fac[i] = mul(fac[i-1],i,pk) : fac[i] = fac[i-1];
}
ll Fac(ll x,ll pi,ll pk)
{
   if(!x) return 1;
   ll res = fac[pk];
   res = mul(qpow(res,x/pk,pk),fac[x%pk],pk);
   return mul(res,Fac(x/pi,pi,pk),pk);
}

ll C(ll n,ll m,ll pi,ll pk)
{
   init(pi,pk);
   ll d = Fac(n,pi,pk),d1 = Fac(m,pi,pk),d2 = Fac(n-m,pi,pk);
   ll k = 0;
   for(ll i = n;i;i /= pi) k += i / pi;
   for(ll i = m;i;i /= pi) k -= i / pi;
   for(ll i = n-m;i;i /= pi) k -= i / pi;
   return mul(mul(d,Inv(d1,pk),pk),mul(qpow(pi,k,pk),Inv(d2,pk),pk),pk);
}

ll exlucas(ll n,ll m)
{
   if(n < m) return 0;
   ll res = CRT(C(n,m,p,p),p);
   return res;
}

ll exlucas1(ll n,ll m)
{
   if(n < m) return 0;
   ll res = 0;
   rep(i,1,5) res += CRT(C(n,m,bas[i],bas[i]),bas[i]),res %= p;
   return res;
}

ll exlucas2(ll n,ll m)
{
   if(n < m) return 0;
   ll res = 0;
   res += CRT(C(n,m,5,125),125),res %= p;
   res += CRT(C(n,m,7,343),343),res %= p;
   res += CRT(C(n,m,101,10201),10201),res %= p;
   return res;
}

int get(int x)
{
   sum = 0;
   int cur = 1,tot = 0;
   while(x)
   {
      if(x&1) sum += A[cur],tot++;
      x >>= 1,cur++;
   }
   return tot;
}

int main()
{
   T = read(),p = read();
   while(T--)
   {
      n = read(),n1 = read(),n2 = read(),m = read();
      rep(i,1,n1+n2) {A[i] = read();if(i > n1) m -= (A[i]-1);}
      if(p == 10007)
      {
     ans = exlucas(m-1,n-1);
     rep(i,1,(1<<n1)-1)
     {
        int d = get(i);
        if(!((d-1)&1)) ans += (p - exlucas(m-1-sum,n-1)),ans %= p;
        else ans += exlucas(m-1-sum,n-1),ans %= p;
     }
     printf("%lld
",ans);
      }
      if(p == 262203414)
      {
     ans = exlucas1(m-1,n-1);
     rep(i,1,(1<<n1)-1)
     {
        int d = get(i);
        if(!((d-1)&1)) ans += (p - exlucas1(m-1-sum,n-1)),ans %= p;
        else ans += exlucas1(m-1-sum,n-1),ans %= p;
     }
     printf("%lld
",ans);
      }
      if(p == 437367875)
      {
     ans = exlucas2(m-1,n-1);
     rep(i,1,(1<<n1)-1)
     {
        int d = get(i);
        if(!((d-1)&1)) ans += (p - exlucas2(m-1-sum,n-1)),ans %= p;
        else ans += exlucas2(m-1-sum,n-1),ans %= p;
     }
     printf("%lld
",ans);
      }
   }
   return 0;
}