【jzoj3214】【SDOI2013】【方程】【中国剩余定理】【组合数取模】

题目大意

给定方程

X1+X 2+…+Xn=m

我们对第 1.. n1 个变量 进行一些限制 ：

X1≤A1

X2≤A2

…

Xn1 ≤An1

我们对第 n1+1.. n1+1.. n1+ n2 个变量 进行一些限制 ：

X_(n1+1)≥A_(n1+1)

X_(n1+2)≥A_(n1+2)

…

X_(n1+n2) ≥A_(n1+n2)

求：在满足这些限制的前提下， 该方程正整数解的个数。

答案可能很大，请输出对 p取模 后的答案 ，也即 答案除以 p的余数。

解题思路

观察可知，没有限制可以用隔板原理组合数解决，后面的限制可以强制分配满足限制。前面的限制可以用容斥原理解决。

关键是这题组合数n，m很大，模数不是质数。这样我们就先把模数分解质因数，对每种质因数分别求模后的组合数，再用中国剩余定理组合答案。

由于组合数要用逆元，对于单独质因子p^x，我们来求阶乘，可以考虑把质因子抽出来考虑，其他的再考虑。考虑将p的倍数先抽出来，前一部分有循环节，后一部分可以递归求解。

$19!\%9=(1*2*3*4*5*6*7*8*9*10*11*12*13*14*15*16*17*18*19) \%9$

$=(1*2*4*5*7*8*10*11*13*14*16*17*19)*3^6*(1*2*3*4*5*6) \%9$

这样阶乘就由两部分组成，不含因子的可以逆元，含因子的直接加减。这样成功求出了对小模数取模的组合数。

我们用中国剩余定理合并，我们要求一个数满足所有对小模数取模得到取模后的数。对于每一个小模数，求一个数对其取模得到想要的答案，对其他小模数取模等于零，我们考虑那个数是想要的数乘以一，但是这个一是由所有除当前模数的模数积乘以它的逆元构成的，由于前者包含的因子包含了其他模数所以取模后为零。将所有这些数加起来就是所要求的答案。

code

#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
#define LF double
#define LL long long
#define min(a,b) ((a<b)?a:b)
#define max(a,b) ((a>b)?a:b)
#define fo(i,j,k) for(int i=j;i<=k;i++)
#define fd(i,j,k) for(int i=j;i>=k;i--)
using namespace std;
int const MxN=1e7;
int N,M,Mo,a[10],Pri[10],Phi[10],Cnt[10],PC[10];LL F[10][10007+10];
int Pow(LL X,int Y,int Mo){
    LL Ans=1;
    while(Y){
        if(Y&1)Ans=Ans*X%Mo;
        X=X*X%Mo;
        Y>>=1;
    }
    return Ans;
}
int FF(int i,int P){
    F[i][P]=1;
    fo(j,1,P)if(j%Pri[i])F[i][P]=F[i][P]*j%PC[i];
    return F[i][P];
}
int NP;
int Fact(int N,LL &A,LL &B,int i){
    A=Pow((F[i][PC[i]])?F[i][PC[i]]:FF(i,PC[i]),(N/PC[i]),PC[i]);B=N/Pri[i];
    NP=N%PC[i];A=A*((F[i][NP])?F[i][NP]:FF(i,NP))%PC[i];
    LL AA=1,BB=0;if(N/Pri[i])Fact(N/Pri[i],AA,BB,i);
    A=A*AA%PC[i];B+=BB;
}
LL A,B,AA,BB;
int CC(int N,int M,int i){
    if(N<M)return 0;
    Fact(N,A,B,i);
    Fact(M,AA,BB,i);
    A=A*Pow(AA,Phi[i]-1,PC[i])%PC[i];B-=BB;
    Fact(N-M,AA,BB,i);
    A=A*Pow(AA,Phi[i]-1,PC[i])%PC[i];B-=BB;
    return (B<Cnt[i])?A*Pow(Pri[i],B,PC[i])%PC[i]:0;
}
int Ans;
int C(int N,int M){
    Ans=0;
    fo(i,1,Pri[0])Ans=(Ans+1ll*Pow(Mo/PC[i],Phi[i],Mo)*CC(N,M,i))%Mo;
    return Ans;
}
int main(){
    freopen("d.in","r",stdin);
    freopen("d.out","w",stdout);
    int T;scanf("%d%d",&T,&Mo);
    int MM=Mo,Mx=sqrt(MM);
    fo(i,2,Mx)if(MM%i==0){
        Pri[++Pri[0]]=i;Cnt[Pri[0]]++;Phi[Pri[0]]=i-1;MM/=i;
        while(MM%i==0)Cnt[Pri[0]]++,Phi[Pri[0]]*=i,MM/=i;
        PC[Pri[0]]=pow(Pri[Pri[0]],Cnt[Pri[0]]);
    }
    if(MM!=1)Pri[++Pri[0]]=MM,Cnt[Pri[0]]++,Phi[Pri[0]]=MM-1,PC[Pri[0]]=pow(Pri[Pri[0]],Cnt[Pri[0]]);
    int N1,N2,X,S,tag,Ans;
    fo(cas,1,T){
        scanf("%d%d%d%d",&N,&N1,&N2,&M);
        fo(i,1,N1)scanf("%d",&a[i]);
        fo(i,1,N2){scanf("%d",&X);M-=max(X-1,0);}
        Mx=(1<<N1)-1;Ans=0;
        fo(s,0,Mx){
            S=s,MM=M,tag=1;
            fo(i,1,N1){
                if(S&1)MM-=a[i],tag=-tag;
                S>>=1;
            }
            Ans=(Ans+tag*C(MM-1,N-1))%Mo;
        }
        printf("%d\n",(Ans+Mo)%Mo);
    }
    return 0;
}