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(\begin{matrix}m\\i\end{matrix}\right)(m-i)^{n-i\times s}\frac{n!}{(s!)^i(n-s*i)!} \]

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

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

则有

\[F_i=\sum_{j=i}^m\left(\begin{matrix}j\\i\end{matrix}\right)G_j \]

二项式反演一下,有

\[G_i=\sum_{j=i}^n(-1)^{j-i}\left(\begin{matrix}j\\i\end{matrix}\right)F_j \]

把组合数拆开

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

可以得到

\[i!G_i=\sum_{j=i}^n\frac{(-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\n",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\n",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\n",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\n",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\n",ans);
    return 0;
}
posted @ 2019-03-19 17:55  dreagonm  阅读(122)  评论(0编辑  收藏  举报