P4491 [HAOI2018]染色

反思

二项式反演的板子，但是还要加上一个NTT，就变得恶心人了

我NTT写挂调了好长时间。。。

注意

NTT应该这么写

int tmp=pow((opt)?invG:G,(MOD-1)/i);

而不这样写（因为单位元是(G^{frac{mod-1}{i}})）

int tmp=pow((opt)?invG:G,n/i);

这里应该这样写

for(int i=2;i<=n;i<<=1){
    ...
}

而不应该这样写（每次长度要倍增）

for(int i=2;i<=n;i++){
    ...
}

多项式还是写的少了

思路

({x_1*y_1,x_2*y_2,x_3*y_3,dots,x_k*y_k})的一个多重集合，设大小为n，它的排列数是

[frac{n!}{x_1!x_2!x_3!dots x_4!} ]

题目让求的是有恰好出现s次的颜色个数恰好是i个的方案数

考虑容斥一发

设(F_i)表示出现s次的颜色个数最少是i个的方案数

考虑钦定i个

则有

[F_i=left(egin{matrix}m\iend{matrix} ight)(m-i)^{n-i imes s}frac{n!}{(s!)^i(n-s*i)!} ]

从前到后分别是选哪i个，后面的(m-i)种颜色随便选(n-i imes s)个位置，剩下的是一个多重集合的排列

设(G_i)表示出现s次的颜色个数恰好是i个的方案

则有

[F_i=sum_{j=i}^mleft(egin{matrix}j\iend{matrix} ight)G_j ]

二项式反演一下，有

[G_i=sum_{j=i}^n(-1)^{j-i}left(egin{matrix}j\iend{matrix} ight)F_j ]

把组合数拆开

[G_i=sum_{j=i}^n(-1)^{j-i}frac{j!}{(j-i)!(i)!}F_j ]

可以得到

[i!G_i=sum_{j=i}^nfrac{(-1)^{j-i}F_jj!}{(j-i)!} ]

设(f_i=F_ii!)，(g_i=frac{(-1)^{j-i}}{(j-i)!})

则

[i!G_i=sum_{j=i}^nf_ig_{j-i} ]

是一个卷积的形式，按照力的做法，把(f)反向成(f')，卷积之后的第(lim-i)项就对应原来的第(i)项

上NTT就好了

代码

#include <cstdio>
#include <algorithm>
#include <cstring>
#define int long long 
using namespace std;
const int G = 3,MOD = 1004535809 , invG = 334845270;
int pow(int a,int b){
    int ans=1;
    while(b){
        if(b&1)
            ans=(ans*a)%MOD;
        a=(a*a)%MOD;
        b>>=1;
    }
    return ans;
}
void NTT(int *a,int n,int opt){
    int lim=0;
    while((1<<lim)<n)
        lim++;
    n=(1<<lim);
    for(int i=0;i<n;i++){
        int t=0;
        for(int j=0;j<lim;j++)
            if((i>>j)&1)
                t|=(1<<(lim-j-1));
        if(i<t)
            swap(a[i],a[t]);
    }
    for(int i=2;i<=n;i<<=1){
        int len=i/2;
        int tmp=pow((opt)?invG:G,(MOD-1)/i);
        for(int j=0;j<n;j+=i){
            int mid=1;
            for(int k=j;k<len+j;k++){
                int t=mid*a[k+len]%MOD;
                a[k+len]=((a[k]-t)%MOD+MOD)%MOD;
                a[k]=(a[k]+t)%MOD;
                mid=(mid*tmp)%MOD;
            }
        }
    }
    if(opt){
        int invn=pow(n,MOD-2);
        for(int i=0;i<n;i++)
            a[i]=(a[i]*invn)%MOD;
    }
}
int a[20001000],b[20001000],Gx[20001000],f[20001000],jc[20001000],jc_inv[20001000],n,m,s,w[20001000];
int C(int n,int m){
    if(n<m)
        return 0;
    return jc[n]*jc_inv[n-m]%MOD*jc_inv[m]%MOD;
}
signed main(){
    // printf("INVG=%lld
",pow(G,MOD-2));
    scanf("%lld %lld %lld",&n,&m,&s);
    for(int i=0;i<=m;i++)
        scanf("%lld",&w[i]);
    int lim=min(n/s,m),maxx=max(n,m);
    jc[0]=1;
    jc_inv[0]=1;
    for(int i=1;i<=maxx;i++)
        jc[i]=(jc[i-1]*i)%MOD;
    jc_inv[maxx]=pow(jc[maxx],MOD-2);
    for(int i=maxx-1;i>=1;i--)
        jc_inv[i]=(jc_inv[i+1]*(i+1))%MOD;
    for(int i=0;i<=min(n/s,m);i++){
        f[i]=(C(m,i)*jc[n]%MOD*pow(pow(jc[s],i),MOD-2)%MOD*(pow(jc[n-i*s],MOD-2))%MOD*pow(m-i,n-i*s)%MOD)%MOD;
        // printf("f[%lld]=%lld
",i,f[i]);
    }
    for(int i=0;i<=lim;i++){
        a[i]=f[i]*jc[i]%MOD;
        b[i]=(((i&1)?-1:1)*jc_inv[i]%MOD+MOD)%MOD;
        // b[i]=jc_inv[lim-i];
        // if((lim-i)&1)
        //     b[i]=MOD-b[i];
        // printf("A[%lld]=%lld B[%lld]=%lld
",i,a[i],i,b[i]);
    }
    for(int i=0,j=lim;i<j;i++,j--)
        swap(a[i],a[j]);
    int ttt=0;
    while((1<<ttt)<(2*lim+4))
        ttt++;
    int lent=(1<<ttt);
    NTT(a,lent,0);
    NTT(b,lent,0);
    for(int i=0;i<lent;i++)
        Gx[i]=(a[i]*b[i])%MOD;
    // for(int i=0;i<lent;i++)
    // printf("A[%lld]=%lld B[%lld]=%lld
",i,a[i],i,b[i]);
    NTT(Gx,lent,1);
    int ans=0;
    for(int i=0;i<=lim;i++)
        ans=(ans+Gx[lim-i]*jc_inv[i]%MOD*w[i]%MOD)%MOD;
    printf("%lld
",ans);
    return 0;
}