[洛谷P4841]城市規劃

摘要： 我真的菜，這道題從冬令營day0開始想起想到現在才想出來，然後發現真的是一道多項式求逆板子題。我已經菜出一種境界了。 下面分享一下我做這道題的經歷（歡迎大家來嘲諷我）： 首先看到這道題我就想了一個神奇的$O(n^3)$的DP，貌似和$n!$暴力一樣只能拿20分，感覺特別不...

我真的菜，這道題從冬令營day0開始想起想到現在才想出來，然後發現真的是一道多項式求逆板子題。我已經菜出一種境界了。

下面分享一下我做這道題的經歷（歡迎大家來嘲諷我）：

首先看到這道題我就想了一個神奇的$O(n^3)$的DP，貌似和$n!$暴力一樣只能拿20分，感覺特別不爽。這個做法是這樣的：

考慮我們做Dinic的時候一開始對一個圖分的層（可能是我網路流做多了腦子做傻了），然後可以發現一張任意的無向簡單連通圖都可以以1為源點唯一地表示為一個分層的模式。所以我們DP一下這個分層的模式就可以了。我們發現轉移下一層的代價只與上一層的個數有關，所以我們定義狀態：$f(i,j)$表示i個節點，分層之後最後一層有j個的簡單無向聯通圖的個數。然後這個一個$2d/1d$動態規劃，我不知道怎麼優化。我一直覺得這個想法很精妙，所以一直陷在這個深坑裡跳不出來。這非常好地說明了我是傻逼。

我感覺正解狀態應該是一維的，但是我之前的想法怎麼也優化不到一維。最終我惱羞成怒，強行把它搞到一維。那狀態是什麼呢？肯定是$f(i)$表示i個點的無向簡單聯通圖的個數了。那麼我考慮轉移，考慮第i個點是不是割點，如果是割點，那麼列舉割掉i之後1所在聯通塊大小，然後隨便轉移一下就可以了，但是不是割點怎麼辦啊？我怎麼知道啊。

所以我又鴿了，我好菜了。我突然發現貌似可以把所有簡單無向圖的個數算出來，然後減去不聯通的個數。怎麼減呢？列舉1所在聯通塊有多少個點就可以了。於是我非常開心地寫出了一個$1d/1d$動態規劃：

$f(i)=2^{C_i^2}-\displaystyle\sum_{j=1}^{i-1}C_{i-1}^{j-1}f(j)2^{C_{i-j}^2}$

那麼我們把它拆開來：

$f(i)=2^{C_i^2}-(i-1)!\displaystyle\sum_{j=1}^{i-1}\frac{f(j)}{(j-1)!}\times\frac{2^{C_{i-j}^2}}{(i-j)!}$

感覺這個$(i-1)!$特別煩，那麼我們把它除掉：

$\frac{f(i)}{(i-1)!}=\frac{2^{C_i^2}}{(i-1)!}-\displaystyle\sum_{j=1}^{i-1}\frac{f(j)}{(j-1)!}\times\frac{2^{C_{i-j}^2}}{(i-j)!}$

設$A(x)=\frac{f(x)}{(x-1)!}$，$B(x)=\frac{2^{C_x^2}}{(x-1)!}$，$C(x)=\frac{2^{C_x^2}}{x!}$

那麼原式子相當於：

$A(x)=B(x)-A(x)C(x)$

變形一下：

$A(x)=\frac{B(x)}{C(x)+1}$

B和C我們都是知道的，那麼按照這個式子算一下A，然後用A就可以求出$f(n)$了。

程式碼：

\#include<bits/stdc++.h>
using namespace std;
const int MAX_N=1<<18|5,MOD=1004535809;
typedef long long ll;
int inv[MAX_N],fac[MAX_N],fac_inv[MAX_N];
inline int fpow(int x,ll n){
    int ret=1;
    for(;n;n>>=1,x=(ll)x*x%MOD)
        n&1?ret=(ll)ret*x%MOD:0;
    return ret;
}
namespace POLY{
    int w[MAX_N],inv[MAX_N];
    void init(){
        for(int i=0;i<=18;++i) w[1<<i]=fpow(3,(MOD-1)/(1<<i));
        for(int i=0;i<=18;++i) inv[1<<i]=fpow(1<<i,MOD-2);
    }
    inline void NTT(int a[],int n,int t){
        if(t==-1) reverse(a+1,a+n);
        for(int i=0,pos=0;i<n;++i){
            if(i<pos) swap(a[i],a[pos]);
            for(int p=n>>1;(pos^=p)<p;p>>=1);
        }
        static int p[MAX_N]; p[0]=1;
        for(int step=1;step<n;step<<=1){
            int stepx=step<<1,wn=w[stepx];
            for(int i=1;i<step;++i) p[i]=(ll)p[i-1]*wn%MOD;
            for(int i=0;i<n;i+=stepx)
                for(int j=i;j<i+step;++j){
                    int x=a[j],y=(ll)a[j+step]*p[j-i]%MOD;
                    a[j]=x+y>=MOD?x+y-MOD:x+y;
                    a[j+step]=x-y<0?x-y+MOD:x-y;
                }
        }
        if(t==-1) for(int i=0;i<n;++i) a[i]=(ll)a[i]*inv[n]%MOD;
    }
    void get_inv(int a[],int b[],int n){
        static int g[MAX_N];
        b[0]=fpow(a[0],MOD-2);
        for(int sz=2;sz<=n;sz<<=1){
            for(int i=sz>>1;i<sz+sz;++i) b[i]=g[i]=0;
            for(int i=0;i<sz;++i) g[i]=a[i];
            NTT(b,sz+sz,1); NTT(g,sz+sz,1);
            for(int i=0;i<sz+sz;++i)
                b[i]=(ll)b[i]*(2ll+MOD-(ll)b[i]*g[i]%MOD)%MOD;
            NTT(b,sz+sz,-1);
        }
    }
    void times(int a[],int b[],int n){
        NTT(a,n,1); NTT(b,n,1);
        for(int i=0;i<n;++i) a[i]=(ll)a[i]*b[i]%MOD;
        NTT(a,n,-1);
    }
}
void fac_init(int n){
    fac[0]=inv[0]=inv[1]=fac_inv[0]=1;
    for(int i=1;i<=n;++i){
        fac[i]=(ll)fac[i-1]*i%MOD;
        if(i>1) inv[i]=(ll)(MOD-MOD/i)*inv[MOD%i]%MOD;
        fac_inv[i]=(ll)fac_inv[i-1]*inv[i]%MOD;
    }
}
int a[MAX_N],b[MAX_N],c[MAX_N],p[MAX_N];
int main(){
    POLY::init();
    int n,top; scanf("%d",&n);
    fac_init(n);
    for(top=1;top<=n+1;top<<=1);
    for(int i=1;i<=n;++i) p[i]=fpow(2,(ll)i*(i-1)/2);
    for(int i=1;i<=n;++i){
        b[i]=(ll)p[i]*fac_inv[i-1]%MOD;
        c[i]=(ll)p[i]*fac_inv[i]%MOD;
    }
    c[0]=(c[0]+1)%MOD;  
    POLY::get_inv(c,a,top);
    POLY::times(a,b,top+top);
    printf("%lld",(ll)a[n]*fac[n-1]%MOD);
    return 0;
}