[SDOI2013]方程

Description
给定方程
\[ X_1+X_2+...+X_n=M \]
我们对第\(1...n_1\)个变量进行一些限制:
\(X_1\leqslant A_1\)
\(X_2\leqslant A_2\)
\(X_{n_1}\leqslant A_{n_1}\)
我们对第\(n_1+1...n_1+n_2\)个变量进行一些限制:
\(X_{n_1+1}\geqslant A_{n_1+1}\)
\(X_{n_1+2}\geqslant A_{n_1+2}\)
\(X_{n_1+n_2}\geqslant A_{n_1+n_2}\)
求:在满足这些限制的前提下,该方程正整数解的个数。
答案可能很大,请输出对p取模后的答案,也即答案除以p的余数。

Input
输入含有多组数据,第一行两个正整数T,p。T表示这个测试点内的数据组数,p的含义见题目描述。
对于每组数据,第一行四个非负整数n,n1,n2,m。
第二行nl+n2个正整数,表示A1..n1+n2。请注意,如果n1+n2等于0,那么这一行会成为一个空行。

Output
共T行,每行一个正整数表示取模后的答案。

Sample Input
3 10007
3 1 1 6
3 3
3 0 0 5
3 1 1 3
3 3

Sample Output
3
6
0


如果不考虑前\(n_1\)个变量的限制,那么答案为
\[ \binom{M-\sum\limits_{i=n_1+1}^{n_1+n_2}A_i-1}{n-1} \]
现在多了\(n_1\)个限制,但\(n_1\)很小,我们可以暴力容斥,枚举哪些变量不满足限制即可

/*program from Wolfycz*/
#include<cmath>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#define inf 0x7f7f7f7f
#define Fi first
#define Se second
#define MK std::make_pair
typedef long long ll;
typedef unsigned int ui;
typedef unsigned long long ull;
typedef std::pair<int,int> pii;
inline char gc(){
    static char buf[1000000],*p1=buf,*p2=buf;
    return p1==p2&&(p2=(p1=buf)+fread(buf,1,1000000,stdin),p1==p2)?EOF:*p1++;
}
template<typename T>inline T frd(T x){
    int f=1; char ch=gc();
    for (;ch<'0'||ch>'9';ch=gc())   if (ch=='-')    f=-1;
    for (;ch>='0'&&ch<='9';ch=gc()) x=(x<<1)+(x<<3)+ch-'0';
    return x*f;
}
template<typename T>inline T read(T x){
    int f=1;char ch=getchar();
    for (;ch<'0'||ch>'9';ch=getchar())  if (ch=='-')    f=-1;
    for (;ch>='0'&&ch<='9';ch=getchar())    x=(x<<1)+(x<<3)+ch-'0';
    return x*f;
}
inline void print(int x){
    if (x<0)    putchar('-'),x=-x;
    if (x>9)    print(x/10);
    putchar(x%10+'0');
}
template<typename T>inline T min(T x,T y){return x<y?x:y;}
template<typename T>inline T max(T x,T y){return x>y?x:y;}
template<typename T>inline T swap(T &x,T &y){T t=x; x=y,y=t;}
const int N=1e4;
namespace Math{
    int P[10],V[10],C[10],f[10][N+10],tot,SP;
    int mlt(int a,int b,int p){
        int res=1;
        for (;b;b>>=1,a=1ll*a*a%p)  if (b&1)    res=1ll*res*a%p;
        return res;
    }
    void prepare(int p){
        SP=p;
        for (int i=2;i*i<=p;i++){
            if (p%i)    continue;
            P[++tot]=i,V[tot]=1;
            while (p%i==0)  p/=i,V[tot]*=i,C[tot]++;
            f[tot][0]=1;
            for (int j=1;j<=V[tot];j++) f[tot][j]=1ll*f[tot][j-1]*(j%i?j:1)%V[tot];
        }
        if (p!=1){
            f[++tot][0]=1;
            P[tot]=V[tot]=p,C[tot]=1;
            for (int j=1;j<=V[tot];j++) f[tot][j]=1ll*f[tot][j-1]*(j%p?j:1)%V[tot];
        }
    }
    int gcd(int a,int b){return !b?a:gcd(b,a%b);}
    void exgcd(int a,int b,int &x,int &y){
        if (!b){x=1,y=0;return;}
        exgcd(b,a%b,x,y);
        int t=x; x=y,y=t-a/b*y;
    }
    int Ex_GCD(int a,int b,int c){
        int d=gcd(a,b),x,y;
        if (c%d)    return -1;
        a/=d,b/=d,c/=d;
        exgcd(a,b,x,y);
        x=(1ll*x*c%b+b)%b;
        return x;
    }
    pii work(int n,int i){
        if (n<=1)   return MK(1,0);
        int res=1ll*mlt(f[i][V[i]],n/V[i],V[i])*f[i][n%V[i]]%V[i];
        pii tmp=work(n/P[i],i); int cnt=n/P[i];
        return MK(1ll*res*tmp.Fi%V[i],tmp.Se+cnt);
    }
    int calc(int n,int m,int i){
        int res=1,cnt=0; pii tmp;
        tmp=work(n,i); cnt+=tmp.Se;
        res=1ll*res*tmp.Fi;
        
        tmp=work(m,i); cnt-=tmp.Se;
        res=1ll*res*Ex_GCD(tmp.Fi,V[i],1)%V[i];
        
        tmp=work(n-m,i); cnt-=tmp.Se;
        res=1ll*res*Ex_GCD(tmp.Fi,V[i],1)%V[i];
        
        return cnt<C[i]?1ll*res*mlt(P[i],cnt,V[i])%V[i]:0;
    }
    int Ex_C(int n,int m){
        if (n<m)    return 0;
        if (tot==1) return Ex_GCD(1,V[1],calc(n,m,1));
        int Ans=0;
        for (int i=1;i<=tot;i++)    Ans=(1ll*Ex_GCD(SP/V[i],V[i],1)*(SP/V[i])%SP*calc(n,m,i)%SP+Ans)%SP;
        return Ans;
    }
}
int main(){
    int T=read(0),p=read(0);
    Math::prepare(p);
    for (;T;T--){
        static int A[10];
        int n=read(0),x=read(0),y=read(0),m=read(0),Ans=0;
        for (int i=0;i<x;i++)   A[i]=read(0);
        for (int i=0;i<y;i++)   m-=read(0)-1;
        for (int sta=0;sta<1<<x;sta++){
            int del=0,cnt=0;
            for (int i=0;i<x;i++)   if (sta>>i&1)   cnt++,del+=A[i];
            Ans=(Ans+(cnt&1?-1:1)*Math::Ex_C(m-del-1,n-1))%p;
        }
        printf("%d\n",(Ans+p)%p);
    }
    return 0;
}
posted @ 2019-03-31 09:55 Wolfycz 阅读(...) 评论(...) 编辑 收藏