您的位置:首页 > 其它

方程

2015-12-02 20:06 330 查看

题目大意

现有N个未知数X1..N,已知这N个未知数之和为M。

现有n1个不等式,第i个不等式为Xi<=Fi。

还有n2个不等式,第i个不等式为Xi+n1>=Gi。

求方程的正整数解个数,结果模P。

若P=p1c1∗p2c2∗...∗ptopctopP=p1^{c1}*p2^{c2}*...*ptop^{ctop}

则最大的pici<=105pi^{ci}<=10^5

N,M<=10^9,n1,n2<=8。

隔板问题

我们可以将模型转换,即有N个箱子,将M个球分到N个箱子里,并要求满足一些约束,且每个箱子至少有一个球。可知,在没有约束时,答案为CN−1M−1C_{M-1}^{N-1}

对于第二种约束Xi+n1>=Gi,我们可以在对应箱子放Gi-1个球,则无论如何分配,一定能满足这个约束。

因此先令M=M-∑n2i=1Gi−1\sum_{i=1}^{n2}Gi-1

如何对付第一种约束?

容斥原理

我们发现,我们可以使用容斥原理来解决。

不满足第i个约束,即不满足Xi<=Fi,那么就是Xi>Fi,我们可以在对应箱子放Fi个球,则无论如何分配一定不满足这个约束。即可解决此题。

由于N,M即大,但可以利用P的最大pcp^c不超过10510^5来分解质因数并快速阶乘,做组合数模。推荐一道经典的组合数模题reward,我的blog有题解。

参考程序

[code]#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
;
    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);
    }
}
内容来自用户分享和网络整理,不保证内容的准确性,如有侵权内容,可联系管理员处理 点击这里给我发消息
标签: