常係數齊次線性遞推

本博客部分參考自dtz的博客，Troywar的博客

前言：

爲何要來學這個鬼畜的東西呢……由於模擬賽出了一道板題（好吧這實際上就是BZOJ4162）：

看到的時候我很是興奮，這不是一道sb矩陣快速冪題嗎？連鄰接矩陣都直接給出來了；

結果一看數據範圍，算了算$n^3logk=1.25e10$，emmm我咋不知道矩陣快速冪有什麼優化方法？

因而我就自閉了，寫了50滾粗（而後一個點1008ms一個點1114ms少了20，氣傻了）

而後才知道這是個常係數齊次線性遞推模板題……

常係數齊次線性遞推

0.預備知識：

有一個數列遞推式爲：

$$S_i=\sum\limits_{j=1}^{k}a_jS_{i-j}$$

這是一個常係數齊次線性遞推的形式，經典問題是給出前$k$項，要求數列第$n$項。樸素的矩陣快速冪相信你們都會……可是矩陣快速冪時間複雜度是$O(k^3logn)$的（據說有些神祕方法能夠優化到$O(k^{log_27}n)$可是沒啥用），當$k$比較大的時候就GG了，所以須要一些優化。

1.特徵多項式和Cayley-Hamilton定理

設$A$是一個$n$階矩陣，且數$\lambda$和非零列向量$x$知足

$$Ax=\lambda x \;(1)$$

則稱$\lambda$爲矩陣$A$的特徵值，$x$爲矩陣$A$的特徵向量；

稍微變式能夠獲得

$$(A-\lambda I)x=0$$

其中$I$表示單位矩陣；

此式有非零解當且僅當其行列式

$$|A-\lambda I|=0 \;(2)$$

(1)式能夠當作一個關於$\lambda$的一元$n$次方程，稱爲矩陣$A$的特徵方程，(2)式是關於$\lambda$的$n$次多項式，稱爲矩陣$A$的特徵多項式，記爲$\phi(\lambda)$；

注意這裏的「特徵多項式」指的是廣義的多項式，未知數不必定是一個數或字母，也能夠是向量或者矩陣；

咱們稱未知量爲一個矩陣的多項式爲矩陣的多項式（不一樣於矩陣多項式）；

那麼對於任意一個矩陣$A$的特徵多項式$\phi(\lambda)$，都有

$$\phi(A)=0（Caylay-Hamilton定理）$$

證實：$\phi(A)=|A-AI|=|0|=0$

2.求解常係數齊次線性遞推矩陣的特徵多項式

沿用上面的定義，設數列$S$的遞推矩陣$A$的特徵多項式爲$\phi(\lambda)$，則：

$$\phi(\lambda)=|A-\lambda I|=\begin{bmatrix}\lambda&0&\cdots&0&0\\ 0&\lambda&\cdots&0&0\\ \cdots&\cdots&\cdots&\cdots&\cdots\\ 0&0&\cdots&\lambda&0\\ 0&0&\cdots&0&\lambda\end{bmatrix}-\begin{bmatrix}0&0&\cdots&0&a_k\\ 1&0&\cdots&0&a_{k-1}\\ 0&1&\cdots&0&a_{k-2}\\ \cdots&\cdots&\cdots&\cdots&\cdots\\ 0&0&\cdots&1&a_1\end{bmatrix}$$

$$=\begin{bmatrix}\lambda&0&\cdots&0&-a_k\\ -1&\lambda&\cdots&0&-a_{k-1}\\ 0&-1&\cdots&0&-a_{k-2}\\ \cdots&\cdots&\cdots&\cdots&\cdots\\ 0&0&\cdots&-1&\lambda-a_1\end{bmatrix}$$

把$\phi(\lambda)$按最後一列$(j=k)$拉普拉斯展開得

$$\phi(\lambda)=\sum\limits_{i=1}^{k}a_{k-i+1}(-1)^{i+k}\phi(\lambda)_{i,k}$$

其中$\phi(\lambda)_{i,k}$表示$(i,k)$的代數餘子式；

化簡得：

$$\phi(\lambda)=(-1)^k(\lambda^k-\sum\limits_{i=1}^{k}a_i\lambda^{k-i})$$

3.常係數齊次線性遞推

回顧最初的問題，咱們想快速求出$A^{n-k+1}$，但直接矩陣乘法時間複雜度會沒法承受，因此能夠考慮從答案出發，用若干個簡化的式子來分解$A^n$（這裏用$n$代替了$n-k+1$，規模上是等價的）；

根據多項式取模的定義，咱們有：

$$x^n(\mod \phi(x))=\phi(x)g(x)+r(x) \;(3)$$

其中$g(x),r(x)$爲兩個多項式；

把(3)式推廣到矩陣的多項式，有：

$$A^n(\mod \phi(A))=\phi(A)g(A)+r(A)$$

因爲$\phi(A)=0$，$A^n \mod \phi(A)=A^n$，因此：

$$A^n=r(A)$$

由多項式取模的定義知$r(A)$的次數確定小於$\phi(A)$的次數，即小於$k$，因此能夠直接計算；

考慮多項式取模，若$A(x)\mod P(x)=B(x)$，則$A^2(x)\mod P(x)=B^2(x)$；

所以能夠用快速冪的方法來計算$A^n \mod \phi(A)$的結果，即$r(x)$的每項係數；

設已求出$r(x)=\sum\limits_{i=0}^{k-1}b_ix^i$，則

$$r(A)=\sum\limits_{i=0}^{k-1}b_iA^i$$

$$A^n=\sum\limits_{i=0}^{k-1}b_iA^i$$

設$T_i$爲表示數列中連續$k$項的向量$\begin{bmatrix}S_i&S_{i+1}&\cdots&S_{i+k-1}\end{bmatrix}$（就是矩陣乘法每次乘出來的那個列向量），把$n-k+1$代回去得：

$$A^{n-k+1}T_0=T_0\sum\limits_{i=0}^{k-1}b_iA^i=\sum\limits_{i=0}^{k-1}b_iA^iT_i=\sum\limits_{i=0}^{k-1}b_iT^i$$

答案就是$A^{n-k+1}T_0$的最後一位，即：

$$S_n=\sum\limits_{i=0}^{k-1}b_iS_{i+k}$$

顯然$S_k$到$S_{2k-1}$能夠暴力計算，利用上式能夠直接求出答案。

至此問題終於解決了！（完結撒花）

4.時間複雜度

因爲$k$通常不超過2000，因此實現時多項式乘法和取模能夠直接暴力，這樣子時間複雜度是$O(k^2logn)$的；

若是遇到毒瘤題固然也能夠用多項式全家桶，時間複雜度就是$O(klogklogn)$的。

例題

【BZOJ4161】shlw loves matrix I

模板題，就是上面提到的例子

 1 #include<algorithm>
 2 #include<iostream>
 3 #include<cstring>
 4 #include<cstdio>
 5 #include<cmath>
 6 #include<queue>
 7 #define inf 2147483647
 8 #define eps 1e-9
 9 #define mod 1000000007
10 using namespace std; 11 typedef long long ll; 12 typedef double db; 13 int n,k,ans=0,gm[5001],a[5001],h[5001],x[5001],ret[5001]; 14 void add(int &x,int y){ 15     if(x+y>=mod)x=x+y-mod; 16     else x+=y; 17 } 18 void mul(int *s,int *x,int *y){ 19     static int t[5001]; 20     for(int i=0;i<=k*2;i++)t[i]=0; 21     //mul
22     for(int i=0;i<k;i++){ 23         for(int j=0;j<k;j++){ 24             add(t[i+j],(ll)x[i]*y[j]%mod); 25  } 26  } 27     //mod
28     for(int i=k*2-2;i>=k;i--){ 29         for(int j=k-1;j>=0;j--){ 30             add(t[i+j-k],mod-(ll)t[i]*gm[j]%mod); 31  } 32  } 33     for(int i=0;i<k;i++){ 34         s[i]=t[i]; 35  } 36 } 37 void getpw(int y){ 38     for(;y;y>>=1,mul(x,x,x)){ 39         if(y&1)mul(ret,ret,x); 40  } 41 } 42 int main(){ 43     scanf("%d%d",&n,&k); 44     n++; 45     for(int i=1;i<=k;i++){ 46         scanf("%d",&a[i]); 47         if(a[i]<0)a[i]+=mod; 48         gm[k-i]=mod-a[i]; 49  } 50     gm[k]=1; 51     for(int i=1;i<=k;i++){ 52         scanf("%d",&h[i]); 53         if(h[i]<0)h[i]+=mod; 54  } 55     if(n<=k){ 56         printf("%d\n",h[n]); 57         return 0; 58  } 59     for(int i=k+1;i<=k*2;i++){ 60         for(int j=1;j<=k;j++){ 61             add(h[i],(ll)a[j]*h[i-j]%mod); 62  } 63  } 64     x[1]=ret[0]=1; 65     getpw(n-k); 66     for(int i=0;i<k;i++){ 67         add(ans,(ll)ret[i]*h[i+k]%mod); 68  } 69     printf("%d",ans); 70     return 0; 71 } 

【BZOJ4612】shlw loves matrix II

就是開頭的題目，這題的矩陣是任意給出的，所以不能直接求出特徵多項式，要帶入$k+1$個$\lambda$求值以後用拉格朗日插值插出特徵多項式，求出$r(x)$以後再把$A$帶進去暴力算就好了，時間複雜度$O(n^4+n^2logk)$

 1 #include<algorithm>
 2 #include<iostream>
 3 #include<cstring>
 4 #include<cstdio>
 5 #include<cmath>
 6 #include<queue>
 7 #define inf 2147483647
 8 #define eps 1e-9
 9 #define mod 1000000007
 10 using namespace std;  11 typedef long long ll;  12 typedef double db;  13 struct lambda{  14     int id,x;  15  lambda(){}  16     lambda(int _id,int _x){  17         id=_id,x=_x;  18  }  19 }p[55];  20 int n,a[55][55],sq[55][55],tmp[55][55],dt[55][55],x[110],pw[110],t[110],f[55],g[55],ans[55][55];  21 char s[10001];  22 int fastpow(int x,int y){  23     int ret=1;  24     for(;y;y>>=1,x=(ll)x*x%mod){  25         if(y&1)ret=(ll)ret*x%mod;  26  }  27     return ret;  28 }  29 int det(){  30     int ret=1,nw,tmp;  31     for(int i=1;i<=n;i++){  32         nw=i;  33         for(int j=i;j<=n;j++){  34             if(dt[j][i]){  35                 nw=j;  36                 break;  37  }  38  }  39         if(nw!=i){  40             for(int j=i;j<=n;j++){  41  swap(dt[i][j],dt[nw][j]);  42  }  43             ret=-ret;  44  }  45         for(int j=i+1;j<=n;j++){  46             if(dt[j][i]){  47                 tmp=(ll)dt[j][i]*fastpow(dt[i][i],mod-2)%mod;  48                 for(int k=i;k<=n;k++){  49                     dt[j][k]=(dt[j][k]-(ll)dt[i][k]*tmp%mod)%mod;  50  }  51  }  52  }  53         ret=(ll)ret*dt[i][i]%mod;  54  }  55     return (ret+mod)%mod;  56 }  57 void mul(int *s,int *x,int *y){  58     for(int i=0;i<=n*2;i++)t[i]=0;  59     for(int i=0;i<n;i++){  60         for(int j=0;j<n;j++){  61             //add(t[i+j],(ll)x[i]*y[j]%mod);
 62             t[i+j]=(t[i+j]+(ll)x[i]*y[j]%mod)%mod;  63  }  64  }  65     int mo=fastpow(f[n],mod-2);  66     for(int i=n*2-2;i>=n;i--){  67         for(int j=n-1;j>=0;j--){  68             //add(t[i+j-n],mod-(ll)f[j]*t[i]%mod*mo%mod);
 69             t[i+j-n]=(t[i+j-n]-(ll)f[j]*t[i]%mod*mo%mod)%mod;  70  }  71  }  72     for(int i=0;i<n;i++){  73         s[i]=t[i];  74  }  75 }  76 int main(){  77     scanf("%s%d",s,&n);  78     for(int i=1;i<=n;i++){  79         for(int j=1;j<=n;j++){  80             scanf("%d",&a[i][j]);  81  }  82  }  83     for(int i=0;i<=n;i++){  84         memcpy(dt,a,sizeof(a));  85         for(int j=1;j<=n;j++)dt[j][j]-=i;  86         p[i]=lambda(i,det());  87  }  88     for(int i=0;i<=n;i++){  89         for(int j=0;j<=n;j++)g[j]=0;  90         g[0]=p[i].x;  91         for(int j=0;j<=n;j++){  92             if(i!=j){  93                 g[0]=(ll)g[0]*fastpow(p[j].id-p[i].id,mod-2)%mod;  94  }  95  }  96         for(int j=0;j<=n;j++){  97             if(i!=j){  98                 for(int k=n;k;k--){  99                     g[k]=((ll)g[k]*p[j].id%mod-g[k-1])%mod; 100  } 101                 g[0]=(ll)g[0]*p[j].id%mod; 102  } 103  } 104         for(int j=0;j<=n;j++){ 105             f[j]=(f[j]+g[j])%mod; 106  } 107  } 108     x[1]=pw[0]=1; 109     for(int i=strlen(s)-1;i>=0;i--){ 110         if(s[i]=='1')mul(pw,pw,x); 111  mul(x,x,x); 112  } 113     for(int i=1;i<=n;i++){ 114         sq[i][i]=1; 115  } 116     for(int i=0;i<n;i++){ 117         for(int j=1;j<=n;j++){ 118             for(int k=1;k<=n;k++){ 119                 //add(ans[j][k],(ll)sq[j][k]*pw[i]%mod);
120                 ans[j][k]=(ans[j][k]+(ll)sq[j][k]*pw[i]%mod)%mod; 121  } 122  } 123         memset(tmp,0,sizeof(tmp)); 124         for(int j=1;j<=n;j++){ 125             for(int k=1;k<=n;k++){ 126                 for(int l=1;l<=n;l++){ 127                     //add(tmp[j][k],(ll)sq[j][l]*a[l][k]%mod);
128                     tmp[j][k]=(tmp[j][k]+(ll)sq[j][l]*a[l][k]%mod)%mod; 129  } 130  } 131  } 132         memcpy(sq,tmp,sizeof(tmp)); 133  } 134     for(int i=1;i<=n;i++){ 135         for(int j=1;j<=n;j++){ 136             printf("%d ",(ans[i][j]%mod+mod)%mod); 137  } 138         puts(""); 139  } 140     return 0; 141 }