【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的余数。

(⊙ ▽ ⊙)

利用容斥原理，可以将所有要满足的条件一转化为条件二。
而条件二可以利用隔板法，Ans=Cn−1m−1−∑xi。

---

但是，本道题的最大Trick是组合数取模：Lucas定理+中国剩余定理

两个定理

中国剩余定理

内容：
设a1,a2,a3...ak两两互质，且一个数X mod这些数，分别得到m1,m2,m3...,mk。
设Mj=∏ni=1mi(i!=j)，
则X=∏ni=1ai∗Mi∗M−1i。

Lucas定理

内容：
举个例子，比如说要算19! mod 32。
19!=1∗2∗3∗4∗...∗19
把所有3的倍数提取出来：
19!=(1∗2∗4∗5∗7∗8∗10∗11∗13∗14∗16∗17∗19)∗36∗(1∗2∗3∗4∗5∗6)
很显然，(1∗2∗3∗4∗5∗6)可以利用相同的办法算出，
(1∗2∗4∗5∗7∗8∗10∗11∗13∗14∗16∗17∗19)则可以发现，
每32属于一个循环节，于是可以使用快速幂优化。

---

为了可以使用中国剩余定理，
我们先对题目所给的模数p分解质因数，得∏pkii。
对于Cmn，要对p取模；
就让它先对每个pkii取模，这样就可以使用中国剩余定理。

---

问题又变成：Cmn对pk取模了，
Cmn=n!m!(n−m!)，
所以可以对n!,m!,(n−m)!分别使用Lucas定理：
n!⇒t1∗pu1
m!⇒t2∗pu2
(n−m)!⇒t3∗pu3；
最后，模出来的数即为pu1−u2−u3∗t1∗t2−1∗t3−1。

(￣～￣)

#include<stdio.h>
#define ll long long
using namespace std;
const ll maxn=110000;
ll t,n,m,n1,n2,mo,i,j,k;
ll a[maxn],ans,p[maxn],pp[maxn],P[maxn];
ll ch[maxn],fac[maxn];
ll qpower(ll a,ll b,ll mo){
    ll c=1;
    while (b){
        if (b&1) c=a*c%mo;
        a=a*a%mo;
        b>>=1;
    }
    return c;
}
ll exgcd(ll a,ll b,ll &x,ll &y){
    if (b==0){
        x=1;
        y=0;
        return a;
    }
    ll r=exgcd(b,a%b,x,y);
    ll t=x;
    x=y;
    y=t-a/b*y;
    return r;
}
ll ni(ll v,ll mo){
    ll x,y;
    exgcd(v,mo,x,y);
    return (x%mo+mo)%mo;
}
ll L;
void lucas(ll n,ll p,ll P,ll &tmp,ll &u){
    if (n==0) return ; 
    ll i,k=n/P,l=1;
    u+=n/p;
    if (k){
        l=qpower(L,k,P);
    }//循环节计算
    tmp=tmp*l%P;
    for (i=k*P+1;i<=n;i++) if (i%p) tmp=tmp*i%P;//计算剩余部分
    lucas(n/p,p,P,tmp,u);
}
ll count(ll m,ll n,ll p,ll P){
    ll t1=1,t2=1,t3=1,u1=0,u2=0,u3=0;
    L=1;
    for (int i=1;i<=P;i++) if (i%p) L=L*i%P;
    lucas(n,p,P,t1,u1);
    lucas(m,p,P,t2,u2);
    lucas(n-m,p,P,t3,u3);
    return qpower(p,u1-u2-u3,P)*t1%P*ni(t2,P)%P*ni(t3,P)%P;
}
ll china(){
    ll i,j=0;
    for (i=1;i<=p[0];i++){
        ll tmp=mo/P[i];
        j=(j+ch[i]*ni(tmp,P[i])*tmp)%mo;
    }
    return (j%mo+mo)%mo;
}
ll c(ll up,ll down){
    ll i,j,k;
    if (up>down) return 0;
    for (i=1;i<=p[0];i++) ch[i]=count(up,down,p[i],P[i]);
    return china();
}
void solve(ll v,ll l,ll num){
    ll i,j,k;
    if (l>n1){
        ans=(ans+c(n-1,v-1)*(num&1?-1:1))%mo;
        return;
    }
    if (v-a[l]>0) solve(v-(a[l]),l+1,num+1);
    solve(v,l+1,num);
}
void init(){
    ll i,j,k=mo;
    for (i=2;i*i<=k;i++){
        if (k%i==0){
            p[++p[0]]=i;
            P[p[0]]=1;
            while (k%i==0){
                k/=i;
                pp[p[0]]++;
                P[p[0]]*=i;
            }
        }
    }
    if (k>1){
        p[++p[0]]=k;
        pp[p[0]]=1;
        P[p[0]]=k;
    }
}
int main(){
    scanf("%lld%lld",&t,&mo);
    init();
    while (t--){
        scanf("%lld%lld%lld%lld",&n,&n1,&n2,&m);
        ans=0;
        for (i=1;i<=n1;i++) scanf("%d",&a[i]);
        for (i=1;i<=n2;i++) scanf("%d",&j),m-=j-1;
        solve(m,1,0);
        ans=(ans%mo+mo)%mo;
        printf("%lld\n",ans);
    }
}

(⊙v⊙)

1.使用Lucas定理时，可以预处理阶乘；
2.使用Lucas定理时，要另开变量存储循环节的快速幂。