Luogu P3301 [SDOI2013]方程

Link
\(X_i\ge k\)的限制可以通过\(M-(k-1)\)来消掉。
\(X_i\le k\)的限制可以通过容斥转化为\(X_i\le k\)的限制然后跟上面一样消掉。
这样方案数就是一个经典的隔板问题，ExLucas+CRT就行了。

#include<cstdio>
typedef long long i64;
i64 n,n1,n2,m,p,a[21],f[1000007];
int read(){int x;scanf("%d",&x);return x;}
i64 pow(i64 a,i64 k,i64 p){i64 r=1;for(;k;k>>=1,a=a*a%p)if(k&1)r=a*r%p;return r;}
void exgcd(i64 a,i64 b,i64&x,i64&y){!b? (x=1,y=0):(exgcd(b,a%b,y,x),y-=x*(a/b));}
i64 inv(i64 a,i64 p){i64 x,y;exgcd(a,p,x,y);return (x+p)%p;}
i64 fac(i64 n,i64 pi,i64 pk){return !n? 1:(n<pi? f[n]:pow(f[pk-1],n/pk,pk)*f[n%pk]%pk*fac(n/pi,pi,pk)%pk);}    
i64 C(i64 n,i64 m,i64 pi,i64 pk)
{
    if(n<m) return 0;
    f[0]=1;
    for(int i=1;i<=pk;++i) f[i]=(i%pi)? f[i-1]*i%pk:f[i-1];
    i64 a=fac(n,pi,pk),b=fac(m,pi,pk),c=fac(n-m,pi,pk);int k=0;
    for(i64 i=n;i;i/=pi) k+=i/pi;
    for(i64 i=m;i;i/=pi) k-=i/pi;
    for(i64 i=n-m;i;i/=pi) k-=i/pi;
    return a*inv(b,pk)%pk*inv(c,pk)%pk*pow(pi,k,pk)%pk;
}
i64 CRT(i64 a,i64 pk){i64 x=p/pk;return a*x%p*inv(x,pk)%p;}
i64 solve(i64 n,i64 m,i64 pi,i64 pk)
{
    i64 ans=0;
    for(int i=0;i<1<<n1;++i)
    {
	int f=1;i64 sum=n;
        for(int j=0;j<n1;++j) if(i>>j&1) f=-f,sum-=a[j+1];
        ans=(ans+f*C(sum,m,pi,pk)+pk)%pk;
    }
    return ans;
}
i64 ExLucas(i64 n,i64 m,i64 q)
{
    if(n<m) return 0;
    i64 ans=0;
    for(i64 i=2,pk;i*i<=q;++i)
	if(!(q%i))
	{
	    for(pk=1;!(q%i);q/=i,pk*=i);
	    ans=(ans+CRT(solve(n,m,i,pk),pk))%p;
	}
    if(q^1) ans=(ans+CRT(solve(n,m,q,q),q))%p;
    return ans;
}
int main()
{
    int T=read();p=read();
    for(;T;--T)
    {
	n=read(),n1=read(),n2=read(),m=read();
	for(int i=1;i<=n1+n2;++i) a[i]=read();
	for(int i=n1+1;i<=n1+n2;++i) m-=a[i]-1;
	printf("%lld\n",ExLucas(m-1,n-1,p));
    }
}