bzoj千题计划267：bzoj3129: [Sdoi2013]方程

http://www.lydsy.com/JudgeOnline/problem.php?id=3129

如果没有Ai的限制，就是隔板法，C（m-1，n-1）

>=Ai 的限制：m减去Ai

<=Ai 的限制：容斥原理，总数- 至少有一个数>Ai + 至少有两个数>Ai - ……

计算组合数取模，模数虽然很大也不是质数，但是质因数分解后最大的才 10201，所以用扩展卢卡斯即可

注意在用扩展卢卡斯计算阶乘的时候，要预处理不包含当前质因子的阶乘，否则会TLE 3个点

#include<cstdio>
#include<iostream>

using namespace std;

typedef long long LL;

LL p;

int up[9],down[9];

int num;
int PI[10001],PK[10001];

LL fac[10202];

template<typename T>
void read(T &x)
{
    x=0; char c=getchar();
    while(!isdigit(c)) c=getchar();
    while(isdigit(c)) { x=x*10+c-'0'; c=getchar(); }
}

void pre()
{
    LL t=p;
    for(LL i=2;i*i<=p;++i)
        if(!(t%i))
        {
            PI[++num]=i;
            PK[num]=1;
            while(!(t%i)) t/=i,PK[num]*=i;
        }
    if(t>1) 
    {
        PI[++num]=t;
        PK[num]=t;
    }
}

LL Pow(LL a,LL b,LL mod)
{
    LL res=1;
    for(;b;b>>=1,a=a*a%mod)
        if(b&1) res=res*a%mod;
    return res;
}

void exgcd(LL a,LL b,LL &x,LL &y)
{
    if(!b) { x=1; y=0; return; }
    exgcd(b,a%b,y,x); y-=a/b*x;
}

LL get_inv(LL a,LL b)
{
    LL x,y;
    exgcd(a,b,x,y);
    x=(x%b+b)%b;
    return x;
}

LL get_fac(int n,LL pk,LL pi)
{
    if(!n) return 1;
    LL ans=1;
    if(n/pk) ans=Pow(fac[pk],n/pk,pk);
    ans=ans*fac[n%pk]%pk;
    return ans*get_fac(n/pi,pk,pi)%pk;
}

LL exlucas(int n,int m,LL pk,LL pi)
{
    fac[0]=1;
    for(int i=1;i<=pk;++i)
    {
        fac[i]=fac[i-1];
        if(i%pi) fac[i]=fac[i]*i%pk;
    }
    LL fn=get_fac(n,pk,pi);
    LL fm=get_fac(m,pk,pi);
    LL fnm=get_fac(n-m,pk,pi);
    LL k=0;
    for(int i=n;i;i/=pi) k+=i/pi;
    for(int i=m;i;i/=pi) k-=i/pi;
    for(int i=n-m;i;i/=pi) k-=i/pi;
    LL ans=fn*get_inv(fm,pk)%pk*get_inv(fnm,pk)%pk*Pow(pi,k,pk)%pk;
    return ans*(p/pk)%p*get_inv(p/pk,pk)%p;
}

LL get_C(int n,int m)
{
    if(n<m) return 0;
    LL ans=0;
    LL pk;
    for(int i=1;i<=num;++i) 
        ans=(ans+exlucas(n,m,PK[i],PI[i]))%p;
    return ans;
}

int main()
{
    freopen("equation.in","r",stdin);
    freopen("equation.out","w",stdout);
    int T; 
    read(T); read(p);
    pre();
    int n,n1,n2,m;
    int mm,t;
    LL ans=0;
    while(T--)
    {
        read(n); read(n1); read(n2); read(m);
        for(int i=1;i<=n1;++i) read(up[i]);
        for(int i=1;i<=n2;++i) read(down[i]);
        for(int i=1;i<=n2;++i) m-=down[i]-1;
        ans=0;
        for(int i=0;i<(1<<n1);++i)
        {
            mm=m;
            t=0;
            for(int j=1;j<=n1;++j)
                if(i&(1<<j-1)) mm-=up[j],++t;
            t=(t&1) ? -1 : 1;
            ans=(ans+t*get_C(mm-1,n-1)+p)%p;
        }
        cout<<ans<<'
';
    }
}

3129: [Sdoi2013]方程

Time Limit: 30 Sec Memory Limit: 256 MB
Submit: 646 Solved: 375
[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