bzoj5306 [Haoi2018]染色（容斥原理+ntt）

一个不会数数的老年咸鱼选手的学习经历x

首先我很快的得到了一个式子：$C_m^i*\frac{n!}{(s!)^i*(n-is)!}*(m-i)^{n-is}$
但是我说不出它的意义…感觉很重复x，需要容斥x
然后就傻掉了gg

其实可以按套路来，给每一项设一个容斥系数b[i],
令$nn=min(m,n/S)$
即最后答案为
$Ans=\sum\limits_{i=0}^{nn}b[i]*C_m^i*\frac{n!}{(s!)^i*(n-is)!}*(m-i)^{n-is}$
需要满足对于每一个恰为k种的方案的累计贡献系数为w[k]，即
$\sum\limits_{i=0}^kb[i]*C_k^i=w[k]$
从$i=0$开始我们可以不断解出b[i]，复杂度$O(n^2)$，可以通过50分的数据。

我们考虑进一步优化，求恰为i种的方案数f[i]。
我们有$f[i]=C_m^i*\frac{n!}{(s!)^i*(n-is)!}*g(m-i,n-iS)$
其中$g(i,j)$表示把i种颜色放到j个格子中，没有”恰好占了S个格子”的颜色的方案数。可以容斥得到$g[i][j]=\sum\limits_{k=0}^i(-1)^kC_i^k\frac{j!}{(s!)^k(j-sk)!}(i-k)^{j-ks}$
代入原式，用j-i代替j，得到

$$f(i)=\frac{m!}{i!(m-i)!}\frac{n!}{(S!)^i(n-iS)!}\sum\limits_{j=i}^{nn}(-1)^{j-i}\frac{(m-i)!}{(j-i)!(m-j)!}\frac{(n-iS)!}{(S!)^{j-i}(n-jS)!}(m-j)^{n-jS}\\=\frac{m!n!}{i!}\sum\limits_{j=i}^{nn}(-1)^{j-i}\frac{1}{(j-i)!(m-j)!}\frac{1}{(S!)^{j}(n-jS)!}(m-j)^{n-jS}$$
最终答案就是
$$Ans=\sum\limits_{i=0}^{nn}w_i*f(i)=\sum_\limits{i=0}^{nn}\frac{m!n!w_i}{i!}\sum\limits_{j=i}^{nn}(-1)^{j-i}\frac{1}{(j-i)!(m-j)!}\frac{1}{(S!)^{j}(n-jS)!}(m-j)^{n-jS}\\=m!n!\sum_{j=0}^{nn}\frac{(m-j)^{n-jS}}{(m-j)!(S!)^{j}(n-jS)!}\sum_{i=0}^{j}\frac{w_i}{i!}\frac{(-1)^{j-i}}{(j-i)!}$$
发现后半部分就是个卷积，ntt计算即可，前面都可以预处理。
于是$O(mlogm)$解决。

#include <bits/stdc++.h>
using namespace std;
#define ll long long
#define inf 0x3f3f3f3f
#define N 10000010
#define mod 1004535809
#define G 3
inline char gc(){
    static char buf[1<<16],*S,*T;
    if(S==T){T=(S=buf)+fread(buf,1,1<<16,stdin);if(T==S) return EOF;}
    return *S++;
}
inline int read(){
    int x=0,f=1;char ch=gc();
    while(ch<'0'||ch>'9'){if(ch=='-') f=-1;ch=gc();}
    while(ch>='0'&&ch<='9') x=x*10+ch-'0',ch=gc();
    return x*f;
}
int n,m,S,nn,fac[N],inv[N],a[100010<<2],b[100010<<2],ans=0;
inline void dec(int &x,int y){x-=y;if(x<0) x+=mod;}
inline void inc(int &x,int y){x+=y;if(x>=mod) x-=mod;}
inline int ksm(int x,int k){
    int res=1;for(;k;k>>=1,x=(ll)x*x%mod) if(k&1) res=(ll)res*x%mod;return res;
}
namespace fft{
    int n,m,L=0,R[100010<<2];
    inline void ntt(int *a,int f){
        for(int i=0;i<n;++i) if(i<R[i]) swap(a[i],a[R[i]]);
        for(int i=1;i<n;i+=i){
            int wn=ksm(G,f==1?(mod-1)/(2*i):mod-1-(mod-1)/(2*i));
            for(int j=0;j<n;j+=i<<1){
                int w=1;
                for(int k=0;k<i;++k,w=(ll)w*wn%mod){
                    int x=a[j+k],y=(ll)a[j+k+i]*w%mod;
                    a[j+k]=(x+y)%mod;a[j+k+i]=x-y;if(a[j+k+i]<0) a[j+k+i]+=mod;
                }
            }
        }if(f!=-1) return;int invn=ksm(n,mod-2);
        for(int i=0;i<n;++i) a[i]=(ll)a[i]*invn%mod; 
    }
    inline void calc(){
        m=nn+nn;for(n=1;n<=m;n+=n) ++L;
        for(int i=1;i<n;++i) R[i]=R[i>>1]>>1|(i&1)<<L-1;
        ntt(a,1);ntt(b,1);
        for(int i=0;i<n;++i) a[i]=(ll)a[i]*b[i]%mod;
        ntt(a,-1);
    }
}
int main(){
//  freopen("a.in","r",stdin);
    n=read();m=read();S=read();nn=min(n/S,m);
    fac[0]=1;inv[0]=inv[1]=1;
    for(int i=2;i<=max(n,m);++i) inv[i]=(ll)inv[mod%i]*(mod-mod/i)%mod;
    for(int i=1;i<=max(n,m);++i) fac[i]=(ll)fac[i-1]*i%mod,inv[i]=(ll)inv[i-1]*inv[i]%mod;
    for(int i=0;i<=nn;++i) a[i]=(ll)read()*inv[i]%mod,b[i]=(i&1)?mod-inv[i]:inv[i];
    fft::calc();
    for(int i=0;i<=nn;++i)
        inc(ans,(ll)ksm(m-i,n-i*S)*inv[m-i]%mod*inv[n-i*S]%mod*ksm(inv[S],i)%mod*a[i]%mod);
    ans=(ll)ans*fac[n]%mod*fac[m]%mod;
    printf("%d\n",ans);
    return 0;
}

#include <bits/stdc++.h>
using namespace std;
#define ll long long
#define inf 0x3f3f3f3f
#define N 10000010
#define mod 1004535809
inline char gc(){
    static char buf[1<<16],*S,*T;
    if(S==T){T=(S=buf)+fread(buf,1,1<<16,stdin);if(T==S) return EOF;}
    return *S++;
}
inline int read(){
    int x=0,f=1;char ch=gc();
    while(ch<'0'||ch>'9'){if(ch=='-') f=-1;ch=gc();}
    while(ch>='0'&&ch<='9') x=x*10+ch-'0',ch=gc();
    return x*f;
}
int n,m,S,nn,w[100010],fac[N],inv[N];
inline void dec(int &x,int y){x-=y;if(x<0) x+=mod;}
inline void inc(int &x,int y){x+=y;if(x>=mod) x-=mod;}
inline int ksm(int x,int k){
    int res=1;for(;k;k>>=1,x=(ll)x*x%mod) if(k&1) res=(ll)res*x%mod;return res;
}
inline int C(int x,int y){
    return (ll)fac[x]*inv[y]%mod*inv[x-y]%mod;
}
namespace sol1{
    int b[100010],ans=0;
    inline void gao(bool flag){
        if(!flag){
            for(int i=0;i<=nn;++i){
                b[i]=w[i];
                for(int j=0;j<i;++j)
                    dec(b[i],(ll)C(i,j)*b[j]%mod);
            }
        }else for(int i=0;i<=nn;++i) b[i]=(i&1)?mod-w[0]:w[0];
        for(int i=0;i<=nn;++i){
            inc(ans,(ll)b[i]*C(m,i)%mod*fac[n]%mod*ksm(inv[S],i)%mod*inv[n-i*S]%mod*ksm(m-i,n-i*S)%mod);
        }printf("%d\n",ans);
    }
}
int main(){
//  freopen("a.in","r",stdin);
    n=read();m=read();S=read();nn=min(n/S,m);bool flag=1;
    fac[0]=1;inv[0]=inv[1]=1;
    for(int i=2;i<=max(n,m);++i) inv[i]=(ll)inv[mod%i]*(mod-mod/i)%mod;
    for(int i=1;i<=max(n,m);++i) fac[i]=(ll)fac[i-1]*i%mod,inv[i]=(ll)inv[i-1]*inv[i]%mod;
    for(int i=0;i<=m;++i){
        w[i]=read();if(i&&w[i]) flag=0;
    }if(nn<=5000||flag){sol1::gao(flag);return 0;}
}