bzoj 3129

非常好的一道数学题，考察了大量数论和组合数学的知识

在做本题之前强烈建议先完成下列两个背景知识：

①：

bzoj 2142礼物

因为本题的一部分数据需要利用到拓展卢卡斯定理，而礼物是拓展卢卡斯定理的裸题，先做礼物是一个比较好的选择

有困难戳这里https://blog.csdn.net/lleozhang/article/details/82884768

②：

CF451E

本题的核心思想和CF451E完全相同，CF451E稍简单一些，所以先理解这里的思想再做本题会发现难度降了不少

有困难戳这里https://blog.csdn.net/lleozhang/article/details/83590652

接下来的讨论在你掌握了上两个背景知识的基础上进行：

首先我们注意到一点：他要求方程的解为正整数，那么这样不利于操作，所以我们把m减掉一个n，相当于预设每个x至少为1，然后仅需要求新方程的解非负即可

不要忘记修改上下界

那么这样就可以发现：前面1-n1的要求可以与CF451E完全一致，那么可以使用完全相同的方法来处理。

至于下界的问题：我们可以实现预设所有解均满足下界，那么我们可以再将m减掉所有的下界，这样其他所有x都可以随便选了。

这样就可以用于CF451E完全相同的方法解决掉这道题

可是不要忘记：本题的模数很有特点：

对于40%的数据，模数为10007，这是个质数，直接卢卡斯定理求组合数即可

对于另30%的数据，模数为一个大合数，但是可以质因子分解成几个不同的质数的乘积，那么对每个质数跑卢卡斯，最后中国剩余定理合并即可

对于最后30%的数据，模数为一个大合数，而且质因子分解之后有同一质数的幂次，所以只能使用拓展卢卡斯合并了

// luogu-judger-enable-o2
 
#include <cstdio>
#define ll long long
int T,p;
int n,m,n1,n2;
int aa[15];
int bb[15];
ll a[15];
ll s[10];
ll p0[4]={0,5,7,101};
ll pu[4]={0,125,343,10201};
ll num[4]={0,3,3,2};
ll mode[5]={10007,2,3,11,397};
ll inv[20005][5];
ll mul[20005][5];
struct node
{
    ll mi;
    ll val;
};
void init()
{
    for(int i=0;i<=4;i++)
    {
        inv[0][i]=inv[1][i]=1;
        mul[0][i]=mul[1][i]=1;
        for(int j=2;j<mode[i];j++)
        {
            inv[j][i]=(mode[i]-mode[i]/j)*inv[mode[i]%j][i]%mode[i];
        }
        for(int j=2;j<mode[i];j++)
        {
            mul[j][i]=mul[j-1][i]*j%mode[i];
            inv[j][i]=inv[j-1][i]*inv[j][i]%mode[i];
        }
    }
}
void ex_gcd(ll a,ll b,ll &x,ll &y)
{
    if(b==0)
    {
        x=1;
        y=0;
        return;    
    }
    ex_gcd(b,a%b,x,y);
    ll t=x;
    x=y;
    y=t-(a/b)*x;
}
ll C(ll x,ll y,int num)
{
    if(x<y)
    {
        return 0;
    }else if(x==y)
    {
        return 1;
    }
    if(x<mode[num])
    {
        return mul[x][num]*inv[y][num]%mode[num]*inv[x-y][num]%mode[num]; 
    }else
    {
        return C(x%mode[num],y%mode[num],num)*C(x/mode[num],y/mode[num],num)%mode[num];
    }
}
ll china()
{
    ll M=262203414;
    ll ret=0;
    for(int i=1;i<=5;i++)
    {
        ll M0=M/mode[i-1];
        ll x,y;
        ex_gcd(M0,mode[i-1],x,y);
        x=(x%mode[i-1]+mode[i-1])%mode[i-1];
        ret+=x*a[i]%M*M0%M;
        ret%=M;
    }
    return ret;
}
ll solve(ll x,ll y)
{
    for(int i=1;i<=5;i++)
    {
        a[i]=C(x,y,i-1);
    }
    return china();
}
ll pow_mul(ll x,ll y,ll mod)
{
    ll ans=1;
    while(y)
    {
        if(y&1)
        {
            ans*=x;
            ans%=mod;
        }
        x*=x;
        x%=mod;
        y>>=1;
    }
    return ans;
}
ll get_inv(ll a,ll b)
{
    ll x,y;
    ex_gcd(a,b,x,y);
    return (x%b+b)%b;
}
node get_mul(ll x,ll num)
{
    if(x==0)
    {
        return (node){0,1};
    }
    ll ans=1;
    ll p1=x/p0[num],p2=x/pu[num];
    if(p2)
    {
        for(ll i=2;i<pu[num];i++)
        {
            if(i%p0[num])
            {
                ans*=i;
                ans%=pu[num];
            }
        }
        ans=pow_mul(ans,p2,pu[num]);
    }
    for(ll i=pu[num]*p2+1;i<=x;i++)
    {
        if(i%p0[num])
        {
            ans*=i;
            ans%=p;
        }
    }
    node re=get_mul(p1,num);
    return (node){re.mi+x,ans*re.val%p};
}
ll get_C(ll x,ll y,ll num)
{
    if(x<y)
    {
        return 0;
    }else if(x==y)
    {
        return 1;
    }
    node f1=get_mul(x,num);
    node f2=get_mul(y,num);
    node f3=get_mul(x-y,num);
    ll t1=pow_mul(p0[num],f1.mi-f2.mi-f3.mi,pu[num])*f1.val%pu[num];
    ll t2=get_inv(f2.val,pu[num]);
    ll t3=get_inv(f3.val,pu[num]);
    return t1*t2*t3%pu[num];
}
ll china_again()
{
    ll M=p;
    ll ret=0;
    for(int i=1;i<=3;i++)
    {
        ll M0=M/pu[i];
        ll x,y;
        ex_gcd(M0,pu[i],x,y);
        x=(x%pu[i]+pu[i])%pu[i];
        ret+=x*M0%M*a[i]%M;
        if(ret>=M)
        {
            ret-=M;
        }
    }
    return ret;
}
ll ex_lucas(ll x,ll y)
{
    for(int i=1;i<=3;i++)
    {
        a[i]=get_C(x,y,i);
    }
    return china_again();
}
int main()
{
    scanf("%d%d",&T,&p);
    init();
    while(T--)
    {
        scanf("%d%d%d%d",&n,&n1,&n2,&m);
        m-=n;
        for(int i=1;i<=n1;i++)
        {
            scanf("%d",&aa[i]);
            aa[i]--;
        }
        for(int i=1;i<=n2;i++)
        {
            scanf("%d",&bb[i]);
            bb[i]--;
            m-=bb[i];
        }
        if(m<0)
        {
            printf("0\n");
            continue;
        }else if(m==0)
        {
            printf("1\n");
            continue;
        }
        if(p==10007)
        {
            ll ans=0;
            for(int i=0;i<(1<<n1);i++)
            {
                int t1=m;
                int flag=1;
                for(int j=0;j<n1;j++)
                {
                    if((1<<j)&i)
                    {
                        t1-=(aa[j+1]+1);
                        flag=-flag;
                    }
                }
                if(t1<0)
                {
                    continue;
                }
                ans+=flag*C(t1+n-1,n-1,0);
                ans=(ans%p+p)%p;
            }
            printf("%lld\n",ans);
        }else if(p==262203414)
        {
            ll ans=0;
            for(int i=0;i<(1<<n1);i++)
            {
                int t1=m;
                int flag=1;
                for(int j=0;j<n1;j++)
                {
                    if((1<<j)&i)
                    {
                        t1-=(aa[j+1]+1);
                        flag=-flag;
                    }
                }
                if(t1<0)
                {
                    continue;
                }
                ans+=flag*solve(t1+n-1,n-1);
                ans=(ans%p+p)%p;
            }
            printf("%lld\n",ans);
        }else
        {
            ll ans=0;
            for(int i=0;i<(1<<n1);i++)
            {
                int t1=m;
                int flag=1;
                for(int j=0;j<n1;j++)
                {
                    if((1<<j)&i)
                    {
                        t1-=(aa[j+1]+1);
                        flag=-flag;
                    }
                }
                if(t1<0)
                {
                    continue;
                }
                ans+=flag*ex_lucas(t1+n-1,n-1);
                if(ans<0)
                {
                    ans+=p;
                }else if(ans>=p)
                {
                    ans-=p;
                }
            }
            printf("%lld\n",ans);
        }
    }
    return 0;
}