bzoj3129: [Sdoi2013]方程（容斥+exlucas）

传送门
直接$2^{n1}$枚举前面那些强制不满足条件进行容斥。
然后发现答案是一个组合数，然而模数可能是合数。。。
于是要用$exlucas$
代码：

#include<bits/stdc++.h>
#define ri register int
#define fi first
#define se second
using namespace std;
const int rlen=1<<18|1;
inline char gc(){
    static char buf[rlen],*ib,*ob;
    (ib==ob)&&(ob=(ib=buf)+fread(buf,1,rlen,stdin));
    return ib==ob?-1:*ib++;
}
inline int read(){
    int ans=0;
    char ch=gc();
    while(!isdigit(ch))ch=gc();
    while(isdigit(ch))ans=((ans<<2)+ans<<1)+(ch^48),ch=gc();
    return ans;
}
typedef pair<int,int> pii;
typedef long long ll;
int mod;
inline int add(const int&a,const int&b){return a+b>=mod?a+b-mod:a+b;}
inline int dec(const int&a,const int&b){return a>=b?a-b:a-b+mod;}
inline int mul(const int&a,const int&b){return (ll)a*b%mod;}
inline void Add(int&a,const int&b){a=a+b>=mod?a+b-mod:a+b;}
inline void Dec(int&a,const int&b){a=a>=b?a-b:a-b+mod;}
inline void Mul(int&a,const int&b){a=(ll)a*b%mod;}
inline int ksm(int a,int p){int ret=1;for(;p;p>>=1,a=mul(a,a))if(p&1)Mul(ret,a);return ret;}
inline int ksm(int a,ll p,int Mod){int ret=1;for(;p;p>>=1,a=(ll)a*a%Mod)if(p&1)ret=(ll)ret*a%Mod;return ret;}
const int N=25,ldxup=1e6,M=1e6+5;
int n,n1,n2,m,a[N],b[N],cnt[1<<8],fc[M];
vector<pii>md;
inline void init(){
    int t=mod;
    for(ri mt,i=2,up=sqrt(mod);i<=up;++i){
        if(t!=t/i*i)continue;
        mt=1;
        while(t==t/i*i)t/=i,mt*=i;
        md.push_back(pii(i,mt));
    }
    if(t^1)md.push_back(pii(t,t));
}
inline void exgcd(int a,int b,int&x,int&y){
    if(!b){x=1,y=0;return;}
    exgcd(b,a-a/b*b,x,y);
    int t=x;
    x=y,y=t-a/b*y;
}
inline int inv(int a,int p){
    int x,y;
    return exgcd(a,p,x,y),(x%p+p)%p;
}
inline int fac(ll n,int p,int pk){
    if(!n)return 1;
    int mt=1;
    for(ri l=1,r=p-1;l<=pk;l+=p,r+=p)for(ri i=l;i<=r;++i)mt=(ll)mt*i%pk;
    mt=ksm(mt,n/pk,pk);
    for(ri i=2,up=n%pk;i<=up;++i)if(i!=i/p*p)mt=(ll)mt*i%pk;
    return (ll)mt*fac(n/p,p,pk)%pk;
}
inline int C(ll n,ll m,int p,int pk){
    if(n<m)return 0;
    ll cnt=0;
    for(ll up=n;up>=p;up/=p)cnt+=up/p;
    for(ll up=m;up>=p;up/=p)cnt-=up/p;
    for(ll up=n-m;up>=p;up/=p)cnt-=up/p;
    return (ll)fac(n,p,pk)*inv(fac(m,p,pk),pk)%pk*inv(fac(n-m,p,pk),pk)%pk*ksm(p,cnt,pk)%pk;
}
inline int exlucas(ll n,ll m){
    if(n<m)return 0;
    int p=1,r=0;
    for(ri np,nr,t,a,b,i=md.size()-1;~i;--i){
        a=md[i].fi,b=md[i].se;
        t=C(n,m,a,b);
        np=p*b;
        nr=((ll)r*b%np*inv(b,p)%np+(ll)t*p%np*inv(p,b)%np)%np;
        p=np,r=nr;
    }
    return r;
}
inline int calc(vector<int>&a){
    ll sum=0;
    for(ri i=0,up=a.size();i<up;++i){
        sum+=a[i];
        if(sum>m)return 0;
    }
    return exlucas(m-sum-1,n-1);
}
int main(){
    ri tt=read();
    mod=read();
    init();
    while(tt--){
        n=read(),n1=read(),n2=read(),m=read();
        for(ri i=1;i<=n1;++i)a[i]=read();
        for(ri i=1;i<=n2;++i)b[i]=read();
        vector<int>t;
        int ans=0;
        for(ri sta=0,up=1<<n1;sta<up;++sta){
            cnt[sta]=sta?cnt[sta^(sta&-sta)]+1:0;
            t.clear();
            for(ri i=1;i<=n2;++i)t.push_back(b[i]-1);
            for(ri i=1;i<=n1;++i)if((sta>>(i-1))&1)t.push_back(a[i]);
            cnt[sta]&1?Dec(ans,calc(t)):Add(ans,calc(t));
        }
        cout<<ans<<'\n';
    }
    return 0;
}