[HAOI2018]染色

[HAOI2018]染色

標籤： 容斥原理 NTT 多項式求逆

---

題目連接

提供一種與網上大多數博客中不太同樣的作法。

Solution

這題很明顯是容斥原理。
若是咱們直接計算不考慮重複，那麼能夠寫出式子
\[\sum_{i=0}^{min(n/s,m) } C_m^i C_n^{is} {(is)! \over (s!)^i} w_i (m-i)^{n-is} \]

想法是枚舉剛好出現了i種顏色出現s次，而後其餘的隨便取，可是這樣很明顯是會重複的。
由於在剩下的隨便取的過程當中，可能會又出現一些顏色剛好出現了s次。

咱們考慮給出現種數i加一個容斥係數\(f_i\)，顯然，假若有j種顏色出現了s次，那麼按照上面的式子，j的貢獻會被i計算\(C_j^i\)次。
咱們只須要對於每個\(w_i\),都有\[w_i=\sum_{j=0}^i C_i^j f_j\]
那麼這樣算就不會重了。
最終的式子是\[\sum_{i=0}^{min(n/s,m) } C_m^i C_n^{is} {(is)! \over (s!)^i} f_i (m-i)^{n-is} \]

如今惟一的難點就是如何算\(f\)，直接遞推是\(O(n^2)\)的。

注意到是一個卷積的形式，
\[\sum_{j=0} {f_j \over j!} {1 \over (i-j)!} = { w_i\over i!}\]

令\(F(x)=\sum_{i=0}^{\infty} {f_i\over i!}x^i,G(x)=\sum_{i=0}^{\infty} {1\over i!}x^i,W(x)=\sum_{i=0}^{\infty} {w_i \over i!}x^i\)

顯然有\(F(x)G(x)=W(x)\),即\(F(x)={W(x) \over G(x)}\)

因此只須要分治FFT或者多項式求逆便可。

Code

#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<cmath>
#include<queue>
#include<stack>
#include<set>
#include<map>
using namespace std;
#define ll long long
#define RG register
#define REP(i,a,b) for(RG int i=(a),_end_=(b);i<=_end_;i++)
#define DREP(i,a,b) for(RG int i=(a),_end_=(b);i>=_end_;i--)
#define EREP(i,a) for(int i=start[(a)];i;i=e[i].next)
inline int read()
{
    int sum=0,p=1;char ch=getchar();
    while(!(('0'<=ch && ch<='9') || ch=='-'))ch=getchar();
    if(ch=='-')p=-1,ch=getchar();
    while('0'<=ch && ch<='9')sum=sum*10+ch-48,ch=getchar();
    return sum*p;
}

const int maxn=3e5+20;
const int mod=1004535809;
const int maxp=1e7;

int n,m,s,g[maxn];
int jc[maxp+20],inv[maxp+20],jcn[maxp+20];

int rev[maxn];
int st[maxn];

inline int power(int a,int b)
{
    int ans=1;
    while(b)
    {
        if(b&1)ans=(ll)ans*a%mod;
        b>>=1;
        a=(ll)a*a%mod;
    }
    return ans;
}

inline void NTT(int *p,int n,int op)
{
    int l=0;
    while((1<<l)<n)l++;
    REP(i,0,n-1)rev[i]=(rev[i>>1]>>1)|((i&1)<<l-1);
    REP(i,1,n-1)if(i<rev[i])swap(p[i],p[rev[i]]);
    for(int i=1;i<n;i<<=1)
    {
        int w=power(3,(mod-1)/(i<<1));
        st[0]=1;REP(j,1,i-1)st[j]=(ll)st[j-1]*w%mod;
        for(int j=0;j<n;j+=i<<1)
        {
            for(int k=0;k<i;k++)
            {
                int x=p[j+k],y=(ll)p[i+j+k]*st[k]%mod;
                p[j+k]=x+y;if(p[j+k]>=mod)p[j+k]-=mod;
                p[i+j+k]=x-y;if(p[i+j+k]<0)p[i+j+k]+=mod;
            }
        }
    }
    if(op==-1)
    {
        int inv=power(n,mod-2)%mod;
        REP(i,0,n-1)p[i]=(ll)p[i]*inv%mod;
        reverse(p+1,p+n);
    }
}

inline void init()
{
    n=read();m=read();s=read();
    REP(i,0,min(m,n/s))g[i]=read();
    jc[0]=jc[1]=jcn[0]=jcn[1]=inv[1]=1;
    REP(i,2,max(n,m))jc[i]=(ll)i*jc[i-1]%mod,inv[i]=(ll)(mod-mod/i)*inv[mod%i]%mod,jcn[i]=(ll)inv[i]*jcn[i-1]%mod;
    
}

int f[maxn],h[maxn],H[maxn];
int A[maxn],B[maxn];

void Inv(int *p,int *q,int n)
{
    if(n==1)return q[0]=p[0],void();
    Inv(p,q,n>>1);
    REP(i,0,(n<<1)-1)A[i]=B[i]=0;
    REP(i,0,n-1)A[i]=q[i],B[i]=p[i];
    NTT(A,n<<1,1);NTT(B,n<<1,1);
    REP(i,0,(n<<1)-1)A[i]=(ll)A[i]*A[i]%mod*B[i]%mod;
    NTT(A,n<<1,-1);
    REP(i,0,n-1)q[i]=((ll)2*q[i]-A[i]+mod)%mod;
}

inline void get_f()
{
    int N=1,lim=min(m,n/s);
    while(N<=lim)N<<=1;
    REP(i,0,lim)g[i]=(ll)g[i]*jcn[i]%mod;
    REP(i,0,lim)h[i]=jcn[i];
    Inv(h,H,N);
    N<<=1;
    NTT(H,N,1);NTT(g,N,1);
    REP(i,0,N-1)f[i]=(ll)H[i]*g[i]%mod;
    NTT(f,N,-1);
    REP(i,0,lim)f[i]=(ll)f[i]*jc[i]%mod;
}

inline int C(int n,int m){ if(n<m)return 0;return (ll)jc[n]*jcn[m]%mod*jcn[n-m]%mod;}

inline void doing()
{
    int ans=0,S=1;
    REP(i,0,min(m,n/s))
    {
        ans=(ans+(ll)C(m,i)*C(n,i*s)%mod*jc[i*s]%mod*S%mod*f[i]%mod*power(m-i,n-i*s))%mod;
        S=(ll)S*jcn[s]%mod;
    }
    printf("%d\n",ans);
}

int main()
{
#ifndef ONLINE_JUDGE
    freopen("paint.in","r",stdin);
    freopen("paint.out","w",stdout);
#endif
    init();
    get_f();
    doing();
    return 0;
}