手動博客搬家: 本文發表于20181212 09:37:21, 原地址https://blog.csdn.net/suncongbo/article/details/84962727

嗚啊怎么又是數學了啊。。。數學比例\(\frac{16}{33}=0.4848\)

orz yhx-12243神仙

題目鏈接: https://codeforces.com/contest/947/problem/E

題意:
有一個\([0,n]\)的隨機數\(x\)初始為\(i\)的概率為\(p_i\). \(m\)次操作每次從\([0,x]\)中等概率隨機選擇一個數\(y\)把\(x\)變成\(y\). 對每個\(i=0,1,...,n\)求\(m\)次之后\(x=i\)的概率。

題解:
嗯，是個線代神題。本題解很長，希望大家都能耐心看懂。(當時我看題解也看了三小時多)一、na?ve的dp和矩乘優化
首先，一個\(O(nm)\)的\(dp\)很好想: 設\(f[i][j]\)表示\(i\)輪后\(x=j\)的概率，則輕易列出方程: \(f[i][j]=\sum^{n}_{k=i} \frac{dp[i-1][k]}{k+1}\).
按此\(dp\)時間復雜度\(O(nm)\), 太慢了。考慮優化。
寫出轉移矩陣: 令\((n+1)\)維列向量為\(dp[i]\)數組, 則\(\bm f= \begin{bmatrix} 1 & \frac{1}{2} & \frac{1}{3} & \frac{1}{4} & ... & \frac{1}{n+1}\\ 0 & \frac{1}{2} & \frac{1}{3} & \frac{1}{4} & ... & \frac{1}{n+1}\\ 0 & 0 & \frac{1}{3} & \frac{1}{4} & ... & \frac{1}{n+1} \\ & & & & ... \\ 0 & 0 & 0 & 0 & ... & \frac{1}{n+1} \end{bmatrix}\times \bm f_{i-1}\)
令上面的轉移矩陣為\(\bm A\). 如果直接利用矩陣乘法優化，時間復雜度\(O(n^3\log m)\), 還是太慢了。
即使是使用Hamilton-Cayley定理，復雜度也難以降到\(O(n^2)\)以下。
那么我們該怎么辦呢？于是一個神奇的做法出現了——矩陣對角化。二、對轉移矩陣特征系統的深入了解
在對矩陣進行對角化之前，我們先看一個概念——特征系統(包括特征根和特征向量)。
在這里特征系統相關知識不再闡述，隨便一本線代書上就會講。
那我們觀察轉移矩陣\(\bm A\), 發現一個事實——\((n+1)\)階矩陣\(\rm A\)共有\((n+1)\)個特征根，它們分別是\(1,\frac{1}{2},\frac{1}{3},...,\frac{1}{n+1}\). 這是因為\(\rm A\)是上對角矩陣，其特征多項式為\(f(\lambda)=|\lambda \bm I-\bm A|=\prod^{n+1}_{i=1} (\lambda-\frac{1}{i})\).
求出這個還不夠，要對矩陣進行對角化，我們還需要求出它的特征向量。(什么你說為什么要求？一會就知道了T_T)
很好，我們對較小的矩陣進行手算，得出結論: 第\(i\)個特征根\(\lambda_i=\frac{1}{i+1}\)的特征向量為$$\bm v_i=\begin{bmatrix} (-1)^i & (-1)^{i+1}{i\choose 1} & (-1)^{i+2}{i\choose 2} & ... & (-1)^{i+n}{i\choose n}\end{bmatrix}^T$$即\(v_{i,j}=(-1)^{i+j}{i\choose j}\). (注意這里\(v_{i,j}\)表示的是第\(i\)個特征(列)向量\(\bm v_i\)的第\(j\)個元素，如果放到特征向量構成的矩陣中應該是第\(j\)行第\(i\)列. 為了避免混亂，下文將用\(v_{i,j}\)表示第\(i\)個特征列向量的第\(j\)個元素，\(V_{i,j}\)表示特征向量構成的矩陣中第\(i\)行第\(j\)列的元素，即\(V_{i,j}=v_{j,i}\).)
即: $$\bm V=\begin{bmatrix} 1 & -1 & 1 & -1 & ... & (-1)^n \ 0 & 1 & -2 & 3 & ... & (-1)^{n+1}{n\choose 1} \ 0 & 0 & 1 & -3 & ... & (-1)^{n+2}{n\choose 2} \ & & & & ... \ 0 & 0 & 0 & 0 & ... & (-1)^{n+n}{n\choose n} \end{bmatrix}$$
其中\(\bm V\)為每一個\(\bm v\)列向量一起構成的矩陣。一會我們將看到，這個矩陣將起到非常重要的作用。
嗯，我們剛剛通過對較小的矩陣進行手算的方式得到了特征向量\(\bm V\), 現在我們嘗試通過理論推導去證明這個公式。
我們先形式化特征向量的表達式: \(v_{i,j}=(-1)^{i+j}{i\choose j}, V_{i,j}=(-1)^{i+j} {j\choose i}\)
根據特征向量的定義，對于特征(列)向量\(\bm v_i\)我們需要證明\(\bm A\bm v_i=\lambda_i\bm v_i\).
我們考慮寫出式子: \(\bm A\bm v_i\)是一個\(n\times n\)方陣和\(n\)階列向量的乘積，考慮乘積的第\(r\)個元素就是矩陣\(\bm A\)的第\(r\)行和\(\bm v_i\)的內積，則$$\sum^{n}{j=0} A{r,j}v_{i,j} =\sum^{n}{j=r} \frac{1}{j+1}(-1)^{i+j}{i\choose j} \ =\sum^{i}{j=r} \frac{1}{j+1} (-1)^{i+j}\frac{i!}{j!(i-j)!}\ =\frac{1}{i+1}\sum^{i}{j=r}(-1)^{i+j}\frac{(i+1)!}{(j+1)!(i-j)!}\ =\frac{1}{i+1}\sum^{i}{j=r}(-1)^{i+j}{i+1\choose j+1}\ =\frac{1}{i+1}\sum^{i}{j=r}(-1)^{i+1}{j-i-1\choose j+1}\ =\frac{1}{i+1}(-1)^{i+1}\sum^{i}{j=r}{j-i-1\choose j+1}\ = \frac{1}{i+1}(-1)^{i+1}({i-i-1+1\choose i+1}-{r-i-1\choose r})\ =\frac{1}{i+1}(-1)^i{r-i-1\choose r}\ =(-1)^{i+r}{i\choose r}=\lambda_i v_{i,r}$$
嗯，這一步推導值得仔細體會。一共包含九行，其中第5和9行的原理是\({-a\choose b}={a+b-1\choose b}(-1)^b\); 第7行的原理是對楊輝三角任何一斜列求和，等于最右下角元素的下一個減去最左上角元素的左一個：\(\sum^{n}_{j=0} {a+j\choose b+j}={a+n+1\choose b}-{a+n\choose b-1}\).
下面仔細審視如此奇怪的推導的意圖: 我們發現第三行把\(\frac{1}{j+1}\)巧妙地轉化成了\(\frac{1}{i+1}\), 第五行把\((-1)^{i+j}\)變成了\((-1)^{i+1}\), 然后這樣一個組合數求和的如此難以處理的問題被去掉了雜質，變成了一個純粹的楊輝三角一斜列的求和。這就是要百費周折地各種化成負的再化回來的原因。這一段推導真的有很多精妙的地方，希望自己能夠借鑒。好了，恢復正題。
剛才我們成功證明了特征向量\(\bm v_i\)的表達式，然后這有什么用嗎？

三、矩陣對角化加速運算——從\(O(n^3\log m)\)到\(O(n^2)\)
這里是本題的重點。
首先，我們考慮現在要做的事情：快速計算該矩陣的\(m\)次冪。假設這個矩陣是一個對角矩陣\(\rm diag(x_0,x_1,...,x_n)\)，則它的\(m\)次方非常容易計算: \(\rm diag(x_0^m,x_1^m,...,x_n^m)\).
我們現在希望快速計算矩陣\(\bm f\)的\(m\)次冪，于是我們可以考慮把它轉化為對角矩陣的冪運算。
我們找到一個可逆矩陣\(\bm \Phi\)滿足\(\bm\Phi^{-1}\bm A\bm\Phi=\bm B\)其中\(\bm B\)為對角矩陣。有一個定理告訴我們，\(n\times n\)矩陣\(\bm A\)可對角化當且僅當 \(\bm A\)有\(n\)個線性無關的特征向量。
然后我們考慮如何計算\(\bm A^m\): \(\bm \Phi^{-1}\bm A\bm \Phi=\bm B\)移項可得\(\bm A=\bm\Phi\bm B\bm\Phi^{-1}\), 則\(\bm A^m=(\bm \Phi\bm B\bm \Phi^{-1})^m=\bm \Phi\bm B(\bm \Phi^{-1}\bm \Phi\bm A)^{m-1}\bm\Phi^{-1}=\bm\Phi\bm B^m\bm\Phi^{-1}\), 因此只要求出\(\bm\Phi\bm B^m\bm\Phi^{-1}\)即可。
下面我們考慮如何求\(\bm\Phi\): 有一個定理告訴我們，\(\bm\Phi=\begin{bmatrix}\bm v_1 & \bm v_2 & \bm v_3 & ... & \bm v_n \end{bmatrix}\), 即為\(n\)個列特征向量構成的矩陣。我們可以驗證一下這一點: \(\bm A\bm \Phi=\bm A \begin{bmatrix} \bm v_1 & \bm v_2 & \bm v_3 & ... & \bm v_n \end{bmatrix}=\begin{bmatrix} \bm A\bm v_1 & \bm A\bm v_2 & \bm A\bm v_3 & ... & \bm A\bm v_n \end{bmatrix}=\begin{bmatrix}\lambda_1 \bm v_1 & \lambda_2 \bm v_2 & \lambda_3 \bm v_3 & ... & \lambda_n \bm v_n\end{bmatrix}=\bm \Phi \rm diag(\lambda_1 ,\lambda_2 , \lambda_3 , ... , \lambda_n)\), 從而\(\bm\Phi^{-1}\bm A\bm\Phi=\rm diag(\lambda_1,\lambda_2,\lambda_3,...,\lambda_n)\).
于是，在這道題中，我們構造出橋接矩陣\(\bm \Phi\)和它的逆矩陣，然后進行計算。
根據上面的推導，\(\bm \Phi_{i,j}=(-1)^{i+j} {j\choose i}\), 那如何求出\(\bm \Phi^{-1}\)呢？我們發現其實\(\bm \Phi^{-1}\)就是\(\bm\Phi\)的每一項取絕對值的結果, \(\bm\Phi^{-1}_{i,j}={j\choose i}.\)$$\bm\Phi^{-1}=\begin{bmatrix} 1 & 1 & 1 & 1 & ... & {n\choose 0} \ 0 & 1 & 2 & 3 & ... & {n\choose1} & \ 0 & 0 & 1 & 3 & ... & {n\choose 2} \ & && & ... \ 0 & 0 & 0 & 0 & ... &{n\choose n}\end{bmatrix}$$ 我們可以用二項式反演推出這一點。$$\sum^{n}{k=0}\bm \Phi^{-1}{i,k}\bm\Phi_{k,j}=\sum^{n}{k=0}\bm\Phi^{-1}{i,k}(-1)^{k+j}{j\choose k}\ =\sum^n_{k=0}\bm\Phi^{-1}{i,k}(-1)^{j-k}{j\choose k}=[i=j]\ \bm\Phi^{-1}{i,j}=\sum^{n}_{k=0}[i=k]{j\choose k}={j\choose k}$$
最后一步用到了二項式反演。
就這樣，我們求出了該矩陣的橋接矩陣\(\bm\Phi\)及其逆矩陣\(\bm\Phi^{-1}\), 這樣我們的問題轉化為了將一個矩陣\(\bm f_0\)依次左乘\(\bm\Phi^{-1}, \bm B^m, \bm\Phi\)得到\(\bm f_m\).暴力做是\(O(n^2)\)的。

四、最后的優化——從\(O(n^2)\)到\(O(n\log n)\)
先考慮如何快速左乘\(\bm B^m\). 顯然因為是對角矩陣，直接\(O(n)\)做即可。
然后考慮如何快速左乘\(\bm \Phi^{-1}\): 推一發設\(b_k=\sum^{n}_{i=k} {i\choose k}a_i=\sum^{n}_{i=k}\frac{i!}{k!(i-k)!}a_i\), 則\(k!b_k=\sum^{n}_{i=k}\frac{i!a_i}{(i-k)!}\). 令\(\alpha_k=k!a_k, \beta_k=k!b_k, \gamma_k=k!\), \(\beta_k=\sum^{n}_{i=k}\alpha_i\gamma_{i-k}\)把數組\(\beta\)和\(\alpha\)倒過來可得\(\beta_{n-k}=\sum^{n}_{i=k}\alpha_{n-i}\gamma_{i-k}=\sum_{i+j=n-k}\alpha_i\gamma_j\). 看出來了吧？是個卷積。愉快地使用FFT解決，\(O(n\log n)\).
快速左乘\(\bm \Phi\)同理\(k!b_k=\sum^{n}_{i=0}\frac{(-1)^{i-k}}{(i-k)!}(i!a_i)\)。總時間復雜度\(O(n\log n)\).
問題至此解決！

代碼實現

(說了這么多其實就是fft一下)

//Wrong Coding: //FFT to (dgr<<1) #include<cstdio> #include<cstdlib> #include<cstring> #include<algorithm> #define llong long long #define modinc(x) {if(x>=P) x-=P;} using namespace std; const int N = 1<<19; const int LGN = 19; const int P = 998244353; const int G = 3; llong fact[N+3]; llong finv[N+3]; llong a[N+3]; llong b[N+3]; llong c[N+3]; llong d[N+3]; llong e[N+3]; llong f[N+3]; llong g[N+3]; llong fftid[N+3]; llong tmp1[N+3],tmp2[N+3]; int n; llong m; llong quickpow(llong x,llong y) {llong cur = x,ret = 1ll;for(int i=0; y; i++){if(y&(1ll<<i)) {ret = ret*cur%P; y-=(1ll<<i);}cur = cur*cur%P;}return ret; } llong mulinv(llong x) {return x<=N ? finv[x]*fact[x-1]%P : quickpow(x,P-2);} void init_fftid(int dgr) {int len = 0; for(int i=0; i<=LGN; i++) if(dgr==(1<<i)) {len = i; break;}fftid[0] = 0ll;for(int i=1; i<dgr; i++) fftid[i] = ((fftid[i>>1])>>1)|((i&1)<<(len-1)); } int getdgr(int x) {int ret = 1; while(ret<x) ret<<=1;return ret; } void ntt(int dgr,int coe,llong poly[],llong ret[]) {init_fftid(dgr);for(int i=0; i<dgr; i++) ret[i] = poly[i];for(int i=0; i<dgr; i++) if(i<fftid[i]) swap(ret[i],ret[fftid[i]]);for(int i=1; i<=(dgr>>1); i<<=1){llong tmp = quickpow(G,(P-1)/(i<<1));if(coe==-1) tmp = mulinv(tmp);for(int j=0; j<dgr; j+=(i<<1)){llong expn = 1ll;for(int k=0; k<i; k++){llong x = ret[k+j],y = ret[k+i+j]*expn%P;ret[k+j] = x+y; modinc(ret[k+j]);ret[k+i+j] = x-y+P; modinc(ret[k+i+j]);expn = expn*tmp%P;}}}if(coe==-1){llong tmp = mulinv(dgr);for(int i=0; i<dgr; i++) ret[i] = ret[i]*tmp%P;} } int main() {fact[0] = 1ll;for(int i=1; i<=N; i++) fact[i] = fact[i-1]*i%P;finv[N] = quickpow(fact[N],P-2);for(int i=N-1; i>=0; i--) finv[i] = finv[i+1]*(i+1)%P;scanf("%d%I64d",&n,&m);for(int i=0; i<=n; i++) scanf("%I64d",&a[i]);for(int i=0; i<=n; i++){b[i] = mulinv(i+1);b[i] = quickpow(b[i],m);}int _dgr = n+1,dgr = getdgr(_dgr);for(int i=0; i<_dgr; i++) d[i] = fact[i]*a[i]%P;for(int i=0; i<_dgr-1-i; i++) swap(d[i],d[_dgr-1-i]);for(int i=0; i<_dgr; i++) e[i] = finv[i];for(int i=0; i<_dgr; i++) f[i] = (i&1)==0 ? finv[i] : (P-finv[i])%P;//calculate C = INVPHI*Antt(dgr<<1,1,d,tmp1); ntt(dgr<<1,1,e,tmp2);for(int i=0; i<(dgr<<1); i++) tmp1[i] = tmp1[i]*tmp2[i]%P;ntt(dgr<<1,-1,tmp1,c);for(int i=_dgr; i<(dgr<<1); i++) c[i] = 0ll;for(int i=0; i<_dgr-1-i; i++) swap(c[i],c[_dgr-1-i]);for(int i=0; i<dgr; i++) c[i] = c[i]*finv[i]%P;//calculate C = B*Cfor(int i=0; i<_dgr; i++) c[i] = c[i]*b[i]%P;//calculate G = PHI*Cfor(int i=0; i<_dgr; i++) d[i] = fact[i]*c[i]%P;for(int i=0; i<_dgr-1-i; i++) swap(d[i],d[_dgr-1-i]);ntt(dgr<<1,1,d,tmp1); ntt(dgr<<1,1,f,tmp2);for(int i=0; i<(dgr<<1); i++) tmp1[i] = tmp1[i]*tmp2[i]%P;ntt(dgr<<1,-1,tmp1,g);for(int i=0; i<_dgr-1-i; i++) swap(g[i],g[_dgr-1-i]);for(int i=0; i<dgr; i++) g[i] = g[i]*finv[i]%P;for(int i=0; i<=n; i++) printf("%I64d ",g[i]);return 0; }

后事
——臥槽這E怎么這么難……
——其實這場還有F。

總結

以上是生活随笔為你收集整理的Codeforces 947E/923E Perpetual Subtraction (线性代数、矩阵对角化、DP)的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。