bzoj 3129: [Sdoi2013]方程（容斥原理+组合数学+数论）

来源：互联网发布：中科恒业医学软件编辑：程序博客网时间：2024/05/17 20:46

3129: [Sdoi2013]方程

Time Limit: 30 Sec Memory Limit: 256 MB
Submit: 386 Solved: 235
[Submit][Status][Discuss]

Description

给定方程
X1+X2+. +Xn=M
我们对第l..N1个变量进行一些限制：
Xl < = A
X2 < = A2
Xn1 < = An1
我们对第n1 + 1..n1+n2个变量进行一些限制：
Xn1+l > = An1+1
Xn1+2 > = An1+2
Xnl+n2 > = Anl+n2
求：在满足这些限制的前提下，该方程正整数解的个数。
答案可能很大，请输出对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

【样例说明】
对于第一组数据，三组解为(1，3，2)，(1,4,1),(2,3,1)
对于第二组数据，六组解为(1，1，3)，(1,2,2),(1,3,1),(2,1,2),(2,2,1),(3,1,1)

HINT

n < = 10^9 , n1 < = 8 , n2 < = 8 , m < = 10^9 ,p<=437367875

对于l00%的测试数据： T < = 5，1 < = A1..n1_n2 < = m，n1+n2 < = n

Source

[Submit][Status][Discuss]

题解：对于an1+1~an2这些数，给m减去ai-1就没有影响了，因为限制是Xi>=ai，而插板法的限制是Xi>=1，所以直接减去ai-1就能消除影响
对于a1~an1这些数，只有8个，可以容斥，所有情况-1个不满足的+2个不满足的-...
这里不满足就是强制某些数>ai，也就是给m减去ai，然后用组合数求解。

组合数取模用lucas+中国剩余定理来做就可以了，参见bzoj 2142

#include<iostream>#include<cstdio>#include<cstring>#include<algorithm>#include<cmath>#define N 200003#define LL long long #define pa pair<LL,LL>using namespace std;LL n,m,p;LL prime[N],cnt[N],mod[N],a[N],a1[N],ans1,ans;int t,n1,n2,tot,q[N],len;void get_prime(){    LL x=p;    for (LL i=2;i*i<=x;i++)    if (x%i==0)    {        prime[++tot]=i; cnt[tot]=0; mod[tot]=1;        while (x%i==0)         {            x/=i;            cnt[tot]++;            mod[tot]*=i;         }    }    if (x>1)    {        prime[++tot]=x;        cnt[tot]=1;         mod[tot]=x;    }}int quickpow(LL num,LL x,LL m){    LL ans=1; LL base=num%m;    while(x)    {        if (x&1)  ans=(ans*base)%m;        x>>=1;        base=(base*base)%m;    }    return ans%m;}void exgcd(LL a,LL b,LL &x,LL &y){    if (b==0)    {        x=1; y=0;        return;    }    exgcd(b,a%b,x,y);    LL t=y;    y=x-(a/b)*y;    x=t;}LL inverse(LL a,LL b){    LL x,y;    exgcd(a,b,x,y);    return (x%b+b)%b;}pa solve(int k,LL n){    if (n==0)  return pa(0,1);    LL x=n/prime[k]; LL y=n/mod[k];    LL ans=1;    if (y)    {        for (LL i=2;i<mod[k];i++)         if (i%prime[k]!=0)          ans=(ans*i)%mod[k];        ans=quickpow(ans,y,mod[k]);    }    for (LL i=y*mod[k]+1;i<=n;i++)     if (i%prime[k]!=0)      ans=(ans*i)%mod[k];    pa t=solve(k,x);    return pa(x+t.first,ans*t.second%mod[k]);}LL china(){    LL ans=0,x,y;    for (int i=1;i<=tot;i++)    {        LL r=p/mod[i];        exgcd(r,mod[i],x,y);        ans=(ans+r*x*a[i])%p;    }    return (ans%p+p)%p;}LL calc(int k,LL n,LL m){    if (n<m)  return 0;    pa a=solve(k,n); pa b=solve(k,m); pa c=solve(k,n-m);    LL t1=a.second*quickpow(prime[k],a.first-b.first-c.first,mod[k])%mod[k];    LL t2=inverse(b.second,mod[k])%mod[k];    LL t3=inverse(c.second,mod[k])%mod[k];    return t1*t2*t3%mod[k];}LL work(LL n,LL m){    for (int i=1;i<=tot;i++)     a[i]=calc(i,n,m);    return china();}void dfs(LL x,LL m,LL cnt){    if (x==n1+1){        LL t=work(m-1,n-1);        if (cnt%2) ans=((ans-t)%p+p)%p;        else ans=(ans+t)%p;        return;    }    dfs(x+1,m,cnt);    if (m-a1[x]) dfs(x+1,m-a1[x],cnt+1);}int main(){    scanf("%d%lld",&t,&p);    get_prime();    for (int T=1;T<=t;T++)    {        scanf("%lld%d%d%lld",&n,&n1,&n2,&m);        for (int i=1;i<=n1;i++)  scanf("%lld",&a1[i]);        for (int j=1;j<=n2;j++)         {            LL x; scanf("%lld",&x);            m-=(x-1);         }        ans=0;        if (m-1<n-1)         {            printf("0\n");            continue;         }        dfs(1,m,0);        printf("%lld\n",(ans%p+p)%p);    }}

0 0