[bzoj3129][SDOI2013]方程

题目大意

现有N个未知数X1..N，已知这N个未知数之和为M。
现有n1个不等式，第i个不等式为Xi<=Fi。
还有n2个不等式，第i个不等式为Xi+n1>=Gi。
求方程的正整数解个数，结果模P。
若$P=p1^{c1}*p2^{c2}*...*ptop^{ctop}$
则最大的$pi^{ci}<=10^5$
N,M<=10^9，n1,n2<=8。

隔板问题

我们可以将模型转换，即有N个箱子，将M个球分到N个箱子里，并要求满足一些约束，且每个箱子至少有一个球。可知，在没有约束时，答案为$C_{M-1}^{N-1}$
对于第二种约束Xi+n1>=Gi，我们可以在对应箱子放Gi-1个球，则无论如何分配，一定能满足这个约束。
因此先令M=M-$\sum_{i=1}^{n2}Gi-1$
如何对付第一种约束？

容斥原理

我们发现，我们可以使用容斥原理来解决。
不满足第i个约束，即不满足Xi<=Fi，那么就是Xi>Fi，我们可以在对应箱子放Fi个球，则无论如何分配一定不满足这个约束。即可解决此题。
由于N,M即大，但可以利用P的最大$p^c$不超过$10^5$来分解质因数并快速阶乘，做组合数模。推荐一道经典的组合数模题reward，我的blog有题解。

参考程序

#include<cstdio>
#include<iostream>
#define fo(i,a,b) for(i=a;i<=b;i++)
using namespace std;
typedef long long ll;
ll f[10],a[20],b[20],c[20],d[20],e[20],pri[32000+10],fac[100000+10];
bool bz[32000+10];
ll i,j,k,l,t,n,m,n1,n2,ca,p,pp,num,top,xx,yy,cnt;
ll quicksortmi(ll x,ll y,ll p){
    if (!y) return 1;
    if (y==1) return x%p;
    ll t=quicksortmi(x,y/2,p);
    t=t*t%p;
    if (y%2) t=t*(x%p)%p;
    return t;
}
void gcd(ll a,ll b){
    if (!b){
        xx=1;
        yy=0;
    }
    else{
        gcd(b,a%b);
        swap(xx,yy);
        yy-=xx*(a/b);
    }
}
ll getny(ll x,ll y){
    gcd(x,y);
    xx=(xx%y+y)%y;
    return xx;
}
ll calcfac(ll n,ll p,ll pp){
    if (n<pp) return fac[n];
    ll t=quicksortmi(fac[p-1],n/p,p);
    t=t*fac[n%p]%p;
    cnt+=n/pp;
    t=t*calcfac(n/pp,p,pp)%p;
    return t;
}
ll calc(ll x,ll y,ll p,ll pp){
    ll i;
    fac[0]=1;
    fo(i,1,p-1)
        if (i%pp==0) fac[i]=fac[i-1];
        else fac[i]=fac[i-1]*i%p;
    cnt=0;
    ll A=calcfac(y,p,pp);
    ll tot=cnt;
    cnt=0;
    ll B=calcfac(x,p,pp);
    B=B*calcfac(y-x,p,pp)%p;
    B=getny(B,p);
    return A*B%p*quicksortmi(pp,tot-cnt,p)%p;
}
ll comb(ll x,ll y,ll p){
    if (x>y) return 0;
    fo(i,1,top) a[i]=calc(x,y,d[i],e[i]);
    fo(i,1,top) b[i]=getny(c[i],d[i]);
    ll t=0;
    fo(i,1,top) t=(t+a[i]*b[i]%p*c[i]%p)%p;
    return t;
}
void dfs(ll x,ll m,ll cnt){
    if (x==n1+1){
        ll t=comb(n-1,m-1,p);
        if (cnt%2) num=((num-t)%p+p)%p;
        else num=(num+t)%p;
        return;
    }
    dfs(x+1,m,cnt);
    if (m-f[x]) dfs(x+1,m-f[x],cnt+1);
}
int main(){
    fo(i,2,32000){
        if (!bz[i]) pri[++k]=i;
        fo(j,1,k){
            if (pri[j]*i>32000) break;
            bz[i*pri[j]]=1;
            if (i%pri[j]==0) break;
        }
    }
    scanf("%lld%lld",&ca,&p);
    pp=p;
    fo(i,1,k){
        if (pp%pri[i]==0){
            d[++top]=1;e[top]=pri[i];
            while (pp%pri[i]==0){
                d[top]*=pri[i];
                pp/=pri[i];
            }
        }
    }
    fo(i,1,top) c[i]=p/d[i];
    while (ca--){
        scanf("%lld%lld%lld%lld",&n,&n1,&n2,&m);
        fo(i,1,n1) scanf("%lld",&f[i]);
        fo(i,1,n2){
            scanf("%lld",&k);
            if (k) m-=k-1;
        }
        num=0;
        dfs(1,m,0);
        printf("%lld\n",num);
    }
}