5322 Hope [CDQ分治FFT加速計算dp]

---

題意

每一個每一個排列,排列中每個數(shù)向它后面第一個比它大的數(shù)連一條邊.

每個排列對于答案的貢獻是這個排列所生成的圖中的每一個聯(lián)通量中點的個數(shù)的平方之積.

例如:排列

$1,2,3,6,4,5$

其中

$1,2,3,6$形成一個大小為$4$的聯(lián)通分量.

$4,5$形成一個大小為$2$的聯(lián)通分量.

那么這個排列的貢獻就是$4^2?2^2=64$.

我們需要求所有排列的貢獻.

題解

這種題做法一般就是$dp$.

我們枚舉排列中最大的數(shù)$n$的所在的位置,當(dāng)$n$的位置固定以后,所有在$n$前面的數(shù)都會與$n$形成一個連通分量.而后面的所有的數(shù)就組成了一個子問題.

定義$dp[n]$表示$n$的所有排列所形成的貢獻之和.

我們可以得到遞推公式.

$dp[n]={\sum }_{i=1}^n{C}_{n?1}^{i?1}?(i?1)!?i^2?dp[n?i]$

化簡得到

$dp[n]=(n?1)!{\sum }_{i=1}^n{i}^2?\frac{dp[n?i]}{(n?i)!}$

記$F[i]=\frac{dp[i]}{i!},G[i]=i^2$

那么

$dp[n]=(n?1)!{\sum }_{i=1}^{n}G[i]F[n?i]$

觀察到$dp[n]$的遞推公式是一個卷積的形式.

求卷積我們顯然想到了$FFT/NTT$來加速$dp$的計算,對于要求出所有的$dp[n]$,因此我們需要用$CDQ$分治來輔助計算$FFT/NTT$.

采用分治策略,即將區(qū)間$F[l,r]$分成區(qū)間$F[l,mid]$和$F[mid+1,r]$兩個區(qū)間.

用$F[l,mid]$區(qū)間的內(nèi)容和$G$多項式做卷積,即可計算出$F[l,mid]$部分對于$F[mid+1,r]$的影響.注意!這里要求$F[l,mid]$必須已經(jīng)完整的被計算出來了.隨后再地歸計算$F[mid+1,r]$部分,解決它們內(nèi)部的貢獻以來.

因此代碼的大致框架就是:

void solve(int l,int r) {if(l == r) return ;int mid = (l + r) >> 1;solve(l,mid);//保證區(qū)間F[l,mid]已經(jīng)被完整的計算出來了Conv();//求卷積,計算F[l,mid]對F[mid+1,r]的貢獻.注意此時F[1,l-1]對F[mid+1,r]的貢獻早已被求過了.//注意在此時,F[mid+1,r]內(nèi)部依賴于F[1,mid]的貢獻已經(jīng)全部被計算出來了,因此下一步求內(nèi)部F[mid+1,r]對F[mid+1,r]的貢獻.solve(mid+1,r); }

舉個栗子

當(dāng)$n=4$時刻
初始$[1,1]$就算完整的被計算出來了.
先求$[1,1]$對$[2,2]$的貢獻,從而得到了完整的$[1,2]$.
然后計算$[1,2]$對$[3,4]$的貢獻,此時$[3,3]$已經(jīng)被完整的計算出來了.
然后算$[3,3]$對于$[4,4]$的貢獻,導(dǎo)致$[4,4]$被完整的計算出來.
至此$[1,4]$都被完整的計算出來了.

注意

其中計算$[l,mid]$對$[mid+1,r]$的貢獻的時候,需要注意:

$x\in [mid+1,r]$

$F[x]\leftarrow F[l+0]?G[x?l]+...+F[l+(mid?l)]?G[x?mid]$

將$F[l],F[l+1],...,F[mid]$填在數(shù)組$a[0],a[1],...,a[mid?l]$的位置.

將$G[0],G[1],G[2],...,G[n]$填在數(shù)組$b[0],b[1],b[2],..,b[n]$的位置.
然后求個卷積$c=a?b$

最后卷積數(shù)組的$c[x?l]$位置的值就是$F[x]$.

描述上稍微有點問題,$F$和$dp$的關(guān)系有點混淆,諸位看我代碼就ok了.

代碼

#include <iostream> #include <algorithm> #include <cstring> #define pr(x) std::cout << #x << ':' << x << std::endl #define rep(i,a,b) for(int i = a;i <= b;++i) #define clr(x) memset(x,0,sizeof(x)) #define setinf(x) memset(x,0x3f,sizeof(x)) #define Max(x,y) x = std::max(x,y) #define Min(x,y) x = std::min(x,y) #define Add(x,y) x = (((x)+(y))%P) #define Sub(x,y) x = (((x)-(y)+P%P) #define Mul(x,y) x = ((x)*(y)%P) typedef long long LL; const int N = 1 << 20; const int P = 998244353; const int G = 3; const int NUM = 20;LL wn[NUM]; LL a[N], b[N];LL quick_mod(LL a, LL b, LL m) {LL ans = 1;a %= m;while(b){if(b & 1){ans = ans * a % m;b--;}b >>= 1;a = a * a % m;}return ans; }void GetWn() {for(int i = 0; i < NUM; i++){int t = 1 << i;wn[i] = quick_mod(G, (P - 1) / t, P);} } void Rader(LL a[], int len) {int j = len >> 1;for(int i = 1; i < len - 1; i++){if(i < j) std::swap(a[i], a[j]);int k = len >> 1;while(j >= k){j -= k;k >>= 1;}if(j < k) j += k;} }void NTT(LL a[], int len, int on) {Rader(a, len);int id = 0;for(int h = 2; h <= len; h <<= 1){id++;for(int j = 0; j < len; j += h){LL w = 1;for(int k = j; k < j + h / 2; k++){LL u = a[k] % P;LL t = w * a[k + h / 2] % P;a[k] = (u + t) % P;a[k + h / 2] = (u - t + P) % P;w = w * wn[id] % P;}}}if(on == -1){for(int i = 1; i < len / 2; i++)std::swap(a[i], a[len - i]);LL inv = quick_mod(len, P - 2, P);for(int i = 0; i < len; i++)a[i] = a[i] * inv % P;} }void Conv(LL a[], LL b[], int n) {NTT(a, n, 1);NTT(b, n, 1);for(int i = 0; i < n; i++)a[i] = a[i] * b[i] % P;NTT(a, n, -1); }LL dp[N],Fac[N],iFac[N];void solve(int l,int r) {if(l == r) return ;int mid = (l + r) >> 1;solve(l,mid);int len = 1;while(len <= (r-l+1)) len <<= 1;rep(i,0,len) {b[i] = 1LL*i*i%P;}rep(i,l,mid) {a[i-l] = dp[i]*iFac[i]%P;}rep(i,mid-l+1,len) a[i] = 0;Conv(a,b,len);rep(i,mid+1,r) {dp[i] = (dp[i] + (a[i-l]*Fac[i-1]%P))%P;}solve(mid+1,r); }int main() {Fac[0] = iFac[0] = 1;rep(i,1,N-1) {Fac[i] = Fac[i-1]*i % P;iFac[i] = quick_mod(Fac[i],P-2,P);}GetWn();dp[0] = 1;solve(0,100000);int n;while(std::cin >> n) {std::cout << dp[n] << std::endl;}return 0; }

總結(jié)

以上是生活随笔為你收集整理的HDU5322 - cdq分治FFT加速dp的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。