HAOI2018染色——容斥

Posted 2021-02-14 ylsoi

题目大意

loj

思路

设(f_i)表示至少出现了i种颜色的方案数

[ egin{aligned} f_i&={m choose i} imes frac{(s imes i)!}{(s!)^{i}} imes {nchoose s imes i} imes (m-i)^{n-s imes i}f_i&={m choose i} imes frac{n!}{(s!)^{i} imes (n-s imes i)!} imes (m-i)^{n-s imes i}\end{aligned} ]

设(g_i)表示恰好出现了i种颜色的方案数，不难得到(g_i)的容斥式子。

[ egin{aligned} g_i=&f_i-sum_{j=i+1}^{m}g_{j} imes {jchoose i}g_i=&f_i-frac{sum_{j=i+1}^{m} imes frac{g_j imes j!}{(j-i)!}}{i!} end{aligned} ]

上面是一个经典的分治FFT形式的式子，使用分治FFT可以做到(nlog^2n)。

其实可以更加优化，上面的式子同样也是一个经典的二项式反演。

[ egin{aligned} f_i&=sum_{j=i}^{m}{jchoose i} imes g_{j}g_i&=sum_{j=i}^{m}(-1)^{j-i} imes {jchoose i} imes f_{j} end{aligned} ]

直接用NTT一遍求出g的值即可。

/*=======================================
 * Author : ylsoi
 * Time : 2019.3.22
 * Problem : loj2527_ntt
 * E-mail : [email protected]
 * ====================================*/
#include<bits/stdc++.h>

#define REP(i,a,b) for(int i=a,i##_end_=b;i<=i##_end_;++i)
#define DREP(i,a,b) for(int i=a,i##_end_=b;i>=i##_end_;--i)
#define debug(x) cout<<#x<<"="<<x<<" "
#define fi first
#define se second
#define mk make_pair
#define pb push_back
typedef long long ll;

using namespace std;

void File(){
    freopen("loj2527_ntt.in","r",stdin);
    freopen("loj2527_ntt.out","w",stdout);
}

template<typename T>void read(T &_){
    _=0; T f=1; char c=getchar();
    for(;!isdigit(c);c=getchar())if(c=='-')f=-1;
    for(;isdigit(c);c=getchar())_=(_<<1)+(_<<3)+(c^'0');
    _*=f;
}

string proc(){
    ifstream f("/proc/self/status");
    return string(istreambuf_iterator<char>(f),istreambuf_iterator<char>());
}

const int maxn=1e7+10;
const int maxm=1e5+10;
const int mod=1004535809;
int n,m,S,w[maxm],ans;
int fac[maxn],ifac[maxn];

int qpow(int x,int y=mod-2){
    int ret=1; x%=mod;
    while(y){
        if(y&1)ret=1ll*ret*x%mod;
        x=1ll*x*x%mod;
        y>>=1;
    }
    return ret;
}

void inc(int &x,int y){
    x+=y;
    if(x>=mod)x-=mod;
    else if(x<0)x+=mod;
}

int C(int x,int y){
    if(x<0 || y<0 || x<y)return 0;
    return 1ll*fac[x]*ifac[y]%mod*ifac[x-y]%mod;
}

namespace Poly{
    int om[maxm<<2],dn[maxm<<2],lim,cnt;
    void dft(int *A,int ty){
        if(ty==-1)reverse(A+1,A+lim);
        REP(i,0,lim-1)if(i<dn[i])swap(A[i],A[dn[i]]);
        for(int len=1;len<lim;len<<=1){
            int w=om[len<<1];
            for(int L=0;L<lim;L+=len<<1){
                int wk=1;
                REP(i,L,L+len-1){
                    int u=A[i],v=1ll*A[i+len]*wk%mod;
                    A[i]=(u+v)%mod;
                    A[i+len]=(u-v)%mod;
                    wk=1ll*wk*w%mod;
                }
            }
        }
        if(ty==-1){
            int inv=qpow(lim);
            REP(i,0,lim-1)A[i]=(1ll*A[i]*inv%mod+mod)%mod;
        }
    }
}

void init(){
    read(n),read(m),read(S);
    REP(i,0,m)read(w[i]);
    fac[0]=1;
    int lim=max(n,m);
    REP(i,1,lim)fac[i]=1ll*fac[i-1]*i%mod;
    ifac[lim]=qpow(fac[lim]);
    DREP(i,lim-1,0)ifac[i]=1ll*ifac[i+1]*(i+1)%mod;
}

int f[maxm<<2],g[maxm<<2];

using namespace Poly;

int main(){
//  File();
    init();
    REP(i,0,m)if(S*i<=n)f[i]=1ll*C(m,i)*fac[S*i]%mod*
        qpow(ifac[S],i)%mod*C(n,S*i)%mod*qpow(m-i,n-S*i)%mod*fac[i]%mod;
    REP(i,0,m)g[i]=1ll*(i&1 ? -1 : 1)*ifac[i]%mod;
    reverse(g,g+m+1);

    lim=1,cnt=0;
    while(lim<=m+m)lim<<=1,++cnt;
    if(!cnt)cnt=1;
    om[lim]=qpow(3,(mod-1)/lim);
    for(int i=lim>>1;i;i>>=1)
        om[i]=1ll*om[i<<1]*om[i<<1]%mod;
    REP(i,0,lim-1)dn[i]=dn[i>>1]>>1|((i&1)<<(cnt-1));

    dft(f,1),dft(g,1);
    REP(i,0,lim-1)g[i]=1ll*g[i]*f[i]%mod;
    dft(g,-1);

    REP(i,0,m)inc(ans,1ll*g[i+m]*w[i]%mod*ifac[i]%mod);
    printf("%d
",ans);

    return 0;
}

/*=======================================
 * Author : ylsoi
 * Time : 2019.3.21
 * Problem : loj2527
 * E-mail : [email protected]
 * ====================================*/
#include<bits/stdc++.h>

#define REP(i,a,b) for(int i=a,i##_end_=b;i<=i##_end_;++i)
#define DREP(i,a,b) for(int i=a,i##_end_=b;i>=i##_end_;--i)
#define debug(x) cout<<#x<<"="<<x<<" "
#define fi first
#define se second
#define mk make_pair
#define pb push_back
typedef long long ll;

using namespace std;

void File(){
    freopen("loj2527.in","r",stdin);
    freopen("loj2527.out","w",stdout);
}

template<typename T>void read(T &_){
    _=0; T f=1; char c=getchar();
    for(;!isdigit(c);c=getchar())if(c=='-')f=-1;
    for(;isdigit(c);c=getchar())_=(_<<1)+(_<<3)+(c^'0');
    _*=f;
}

string proc(){
    ifstream f("/proc/self/status");
    return string(istreambuf_iterator<char>(f),istreambuf_iterator<char>());
}

const int maxn=1e7+10;
const int maxm=1e5+10;
const int mod=1004535809;
int n,m,S,w[maxm],ans;
int fac[maxn],ifac[maxn];

int qpow(int x,int y=mod-2){
    int ret=1; x%=mod;
    while(y){
        if(y&1)ret=1ll*ret*x%mod;
        x=1ll*x*x%mod;
        y>>=1;
    }
    return ret;
}

int C(int x,int y){
    if(x<0 || y<0 || x<y)return 0;
    return 1ll*fac[x]*ifac[y]%mod*ifac[x-y]%mod;
}

void inc(int &x,int y){
    x+=y;
    if(x>=mod)x-=mod;
    else if(x<0)x+=mod;
}

namespace Poly{
    int om[maxm<<2],dn[maxm<<2],lim,cnt;
    void dft(int *A,int ty){
        if(ty==-1)reverse(A+1,A+lim);
        REP(i,0,lim-1)if(i<dn[i])swap(A[i],A[dn[i]]);
        for(int len=1;len<lim;len<<=1){
            int w=om[len<<1];
            for(int L=0;L<lim;L+=len<<1){
                int wk=1;
                REP(i,L,L+len-1){
                    int u=A[i],v=1ll*A[i+len]*wk%mod;
                    A[i]=(u+v)%mod;
                    A[i+len]=(u-v)%mod;
                    wk=1ll*wk*w%mod;
                }
            }
        }
        if(ty==-1){
            int inv=qpow(lim);
            REP(i,0,lim-1)A[i]=1ll*A[i]*inv%mod;
        }
    }
    void multi(int *A,int *B,int *C,int la,int lb){
        lim=1,cnt=0;
        while(lim<=la+lb)lim<<=1,++cnt;
        REP(i,0,lim-1){
            dn[i]=dn[i>>1]>>1|((i&1)<<(cnt-1));
            if(i>la)A[i]=0;
            if(i>lb)B[i]=0;
        }
        dft(A,1),dft(B,1);
        REP(i,0,lim-1)C[i]=1ll*A[i]*B[i]%mod;
        dft(C,-1);
    }
}

void math_init(){
    fac[0]=1;
    int lim=max(n,m);
    REP(i,1,lim)fac[i]=1ll*fac[i-1]*i%mod;
    ifac[lim]=qpow(fac[lim],mod-2);
    DREP(i,lim-1,0)ifac[i]=1ll*ifac[i+1]*(i+1)%mod;

    using namespace Poly;
    lim=1;
    while(lim<=m+m)lim<<=1;
    om[lim]=qpow(3,(mod-1)/lim);
    for(int i=lim>>1;i;i>>=1)
        om[i]=1ll*om[i<<1]*om[i<<1]%mod;
}

int f[maxm<<2],g[maxm<<2],ap[maxm<<2],bp[maxm<<2];

void divide(int l,int r){
    if(l==r){
        g[l]=(f[l]-1ll*g[l]*ifac[l]%mod)%mod;
        inc(ans,1ll*g[l]*w[l]%mod);
        return;
    }

    using namespace Poly;
    int mid=(l+r)>>1;

    divide(mid+1,r);

    REP(i,0,r-mid-1)ap[i]=1ll*g[i+mid+1]*fac[i+mid+1]%mod;
    REP(i,0,r-l-1)bp[i]=ifac[i+1];
    reverse(bp,bp+r-l);
    multi(ap,bp,ap,r-mid-1,r-l-1);
    REP(i,l,mid)inc(g[i],ap[i-mid-1+r-l]);

    divide(l,mid);
}

int main(){
    //File();

    read(n),read(m),read(S);
    REP(i,0,m)read(w[i]);

    math_init();

    REP(i,0,m)if(S*i<=n)
        f[i]=1ll*C(m,i)*fac[n]%mod*qpow(ifac[S],i)%mod*ifac[n-S*i]%mod*qpow(m-i,n-S*i)%mod;

    divide(0,m);

    printf("%d
",ans);

    return 0;
}