定義

對于數列$f$，如果有遞推式

$$f_n=\sum _{i=1}^k{a}_i\times f_{n?i}(n\ge k)$$

并且已知前$k$項，可稱其為$k$ 階齊次線性遞推數列。

這個算法可以快速求出$f$的第$n$項，復雜度為$O(k^2\log n)$,可優化至$O(k\log k\log n)$

前置知識

矩陣快速冪

不難看出矩陣快速冪可以解決這個問題

我們構造$k$轉移矩陣

$$A=?\begin{array}{cccccc}{a}_1& a_2& a_3& ?& a_{k?1}& a_k\\ 1& 0& 0& ?& 0& 0\\ 0& 1& 0& ?& 0& 0\\ 0& 0& 1& ?& 0& 0\\ ?& ?& ?& ?& ?& ?\\ 0& 0& 0& ?& 1& 0\end{array}?$$

對于一個狀態

$$S=?\begin{array}{c}{f}_{i+k?1}\\ f_{i+k?2}\\ ?\\ ?\\ f_{i+1}\\ f_i\end{array}?$$

左乘上轉移矩陣即可得到下一個狀態

$$AS=?\begin{array}{c}{f}_{i+k}\\ f_{i+k?1}\\ ?\\ ?\\ f_{i+2}\\ f_{i+1}\end{array}?$$

我們只需要求出$A^{n}S$就可以算出答案

直接求的復雜度是$O(k^3\log n)$，不可接受

接下來將對矩陣快速冪進行優化

特征向量、值、多項式、方程

對于一個$k$階方陣$A$，如果有列向量$x$和實數$\lambda$，滿足

$$Ax=\lambda x$$

我們稱$x$為$A$的特征向量，$\lambda$為$A$對于$x$的特征值。

對這個式子變形

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

其中$I$為$k$階單位矩陣

鬼知道我為什么要寫這個

我們記關于$\lambda$的多項式

$$f(\lambda )=|\lambda I?A|$$

為$A$的特征多項式

$f(\lambda )=0$稱為特征方程，不過我們用不上

Caylay-Hamilton 定理

定理內容：對于方陣$A$,設$f(\lambda )$為其特征多項式，有$f(A)=0$

你沒看錯，$f(A)$是矩陣的多項式，就是把矩陣$A$代進去，計算$a_{0}I+a_{1}A+a_2{A}^2+?+a_k{A}^k$得到一個零矩陣

證明有太多前置知識，我還不會，希望有生之年能補上。

代數余子式

余子式：一個矩陣$A$去掉$A_{i,j}$所在的第$i$行第$j$列后的行列式稱為余子式（對就是矩陣樹那樣），記為$a_{i,j}$,是個具體的數。

代數余子式：$(?1{)}^{i+j}{a}_{i,j}$

引理 矩陣的行列式等于它第一行的代數余子式之和。

對于第$1$行第$i$個元素，因為是第一行，所以貢獻的逆序對個數為$i?1$,剛好和代數余子式的$(?1{)}^{i+j}$消掉，去掉第$1$行第$i$列之后順序沒有任何影響，所以乘一個$a_{i,j}$即可

算法流程

考慮之前的矩陣快速冪，最后推到了$A^{n}S$

我們考慮怎么快速求出$A^n$

我們發現$A^n$是一個矩陣的多項式，考慮前面的CH定理

我們設$A$的特征多項式為$f(\lambda )$

那么我們可以把$A^n$寫成

$$A^n=f(A)g(A)+h(A)$$

因為$f(A)=0$,所以

$$A^n=h(A)$$

也就是說，$A^n$和$h(A)$在矩陣意義上等價

那么如果我們知道$A$的特征多項式，我們就可以用一般意義下的多項式取模（具體做法后面講）算出$h(A)$的表達式,并且最高項次數不超過$k?1$

我們設求出的

$$h(x)=\sum _{i=0}^{k?1}{b}_i{x}^i$$

考慮我們要求的

$$A^{n}S$$

因為$A^n=h(A)$,所以等于

$$\sum _{i=0}^{k?1}{b}_i{A}^{i}S$$

$A^{i}S$相當于是第$i+1$個狀態，我們寫成矩陣形式

$$\sum _{i=0}^{k?1}{b}_i?\begin{array}{c}{f}_{i+k?1}\\ f_{i+k?2}\\ ?\\ ?\\ f_{i+1}\\ f_i\end{array}?=A^{n}S=?\begin{array}{c}{f}_{n+k?1}\\ f_{n+k?2}\\ ?\\ ?\\ f_{n+1}\\ f_n\end{array}?$$

仔細觀察，發現我們只需要最下面的元素

所以

$$\sum _{i=0}^{k?1}{b}_i{f}_i=f_n$$

也就是把$b$和題目給的前$k$項對應乘起來就好了

---

現在的問題是如何求特征多項式

$$\lambda I?A=?\begin{array}{ccccccc}\lambda ?a_1& ?a_2& ?a_3& ?a_4& ?& ?a_{k?1}& ?a_k\\ ?1& \lambda & 0& 0& ?& 0& 0\\ 0& ?1& \lambda & 0& ?& 0& 0\\ 0& 0& ?1& \lambda & ?& 0& 0\\ ?& ?& ?& ?& ?& ?& ?\\ 0& 0& 0& 0& ?& \lambda & 0\\ 0& 0& 0& 0& ?& ?1& \lambda \end{array}?$$

根據上面的引理，

$$|\lambda I?A|=(\lambda ?a_1)a_{1,1}+\sum _{i=2}^k(?1{)}^{i+1}(?a_i)a_{1,i}$$

其中$a_{i,j}$為余子式

因為去掉第$1$行和任意一列后剩下的是一個下三角矩陣，并且對角線上刪除的列的左邊是$?1$,右邊是$\lambda$,所以

$$a_{1,i}=(?1{)}^{i?1}{\lambda }^{k?i}$$

所以

$$\begin{gathered} |\lambda I?A|=(\lambda ?a_1){\lambda }^{k?1}+\sum _{i=2}^k{a}_i{\lambda }^{k?i} \\ ={\lambda }^k?\sum _{i=1}^k{a}_i{\lambda }^{k?i} \end{gathered}$$

即

$$f(\lambda )={\lambda }^k?\sum _{i=1}^k{a}_i{\lambda }^{k?i}$$

就是特征多項式

實現

前面的多項式取模部分是挖了坑的

$$A^n=f(A)g(A)+h(A)$$

我們已知$f(A)$為特征多項式，要求出余數$h(A)$

因為$n$很大，所以我們倍增計算并順便取模

常用的暴力取模好寫的多，復雜度$O(k^2\log n)$，具體實現見代碼注釋

BZOJ4161 Shlw loves matrixI

#include <iostream> #include <cstdio> #include <cstring> #include <cctype> #define MAXN 4005//數組開兩倍 using namespace std; int n,k,F[MAXN],R[MAXN],H[MAXN]; const int MOD=1e9+7; typedef long long ll; inline int add(const int& x,const int& y){return x+y>=MOD? x+y-MOD:x+y;} inline void mul(int* f,int* g,int* to) {static int t[MAXN];memset(t,0,sizeof(t));for (int i=0;i<k;i++)for (int j=0;j<k;j++)t[i+j]=add(t[i+j],(ll)f[i]*g[j]%MOD);//求卷積for (int i=(k<<1);i>=k;t[i--]=0)for (int j=1;j<=k;j++) t[i-j]=add(t[i-j],(ll)t[i]*F[j]%MOD);//因為最高項次數為1，后面都是讀入的數的相反數，所以讀入不取負，這里的減改成加即可for (int i=0;i<k;i++) to[i]=t[i]; } inline void qpow(int* f,int p) {static int ans[MAXN];memset(ans,0,sizeof(ans));ans[0]=1;while (p){if (p&1) mul(f,ans,ans);mul(f,f,f);p>>=1;}for (int i=0;i<k;i++) f[i]=ans[i]; } int main() {scanf("%d%d",&n,&k);for (int i=1;i<=k;i++) scanf("%d",&F[i]),F[i]<0&&(F[i]+=MOD);for (int i=0;i<k;i++) scanf("%d",&H[i]),H[i]<0&&(H[i]+=MOD);R[1]=1;qpow(R,n);int ans=0;for (int i=0;i<k;i++) ans=add(ans,(ll)R[i]*H[i]%MOD);printf("%d\n",ans);return 0; }

---

如果不幸遇到了毒瘤出題人，則需要用多項式全家桶中的多項式取模

復雜度為$O(k\log k\log n)$,常數較大

Luogu P4723

#include <iostream> #include <cstdio> #include <cstring> #include <cctype> #include <algorithm> #define MAXN 131072 using namespace std; inline int read() {int ans=0,f=1;char c=getchar();while (!isdigit(c)) (c=='-')&&(f=-1),c=getchar();while (isdigit(c)) ans=(ans<<3)+(ans<<1)+(c^48),c=getchar();return f*ans; } const int MOD=998244353; typedef long long ll; inline int add(const int& x,const int& y){return x+y>=MOD? x+y-MOD:x+y;} inline int dec(const int& x,const int& y){return x<y? x-y+MOD:x-y;} inline int qpow(int a,int p) {int ans=1;while (p){if (p&1) ans=(ll)ans*a%MOD;a=(ll)a*a%MOD;p>>=1;}return ans; } #define inv(x) qpow(x,MOD-2) int rt[2][24]; int l,lim,r[MAXN]; inline void init(){lim=1<<l;for (int i=0;i<lim;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));} inline void NTT(int* a,int type) {for (int i=0;i<lim;i++) if (i<r[i]) swap(a[i],a[r[i]]);for (int L=0;L<l;L++){int mid=1<<L,len=mid<<1,Wn=rt[type][L+1];for (int s=0;s<lim;s+=len)for (int k=0,w=1;k<mid;k++,w=(ll)w*Wn%MOD){int x=a[s+k],y=(ll)w*a[s+mid+k]%MOD;a[s+k]=add(x,y);a[s+mid+k]=dec(x,y);}}if (type){int t=inv(lim);for (int i=0;i<lim;i++) a[i]=(ll)a[i]*t%MOD;} } void getinv(int* A,int* B,int n) {static int f[MAXN],t[MAXN];if (n==1) return (void)(*B=inv(*A));getinv(A,t,(n+1)>>1);for (l=0;(1<<l)<(n<<1);l++);init();for (int i=0;i<n;i++) f[i]=A[i];for (int i=n;i<lim;i++) f[i]=t[i]=0;NTT(f,0);NTT(t,0);for (int i=0;i<lim;i++) B[i]=(ll)t[i]*dec(2,(ll)f[i]*t[i]%MOD)%MOD;NTT(B,1);for (int i=n;i<lim;i++) B[i]=0; } void getmod(int* A,int* B,int n,int m) {static int f[MAXN],g[MAXN],h[MAXN],t[MAXN];for (int i=0;i<=n;i++) f[i]=A[n-i];for (int i=0;i<=m;i++) g[i]=B[m-i];getinv(g,t,n-m+1);for (int i=n-m+1;i<lim;i++) f[i]=0;NTT(f,0);NTT(t,0);for (int i=0;i<lim;i++) h[i]=(ll)f[i]*t[i]%MOD;NTT(h,1);for (int i=n-m+1;i<lim;i++) h[i]=0;reverse(h,h+n-m+1),reverse(g,g+m+1);for (l=0;(1<<l)<=n;l++);init();for (int i=m+1;i<lim;i++) g[i]=0;for (int i=n-m+1;i<lim;i++) h[i]=0;NTT(g,0);NTT(h,0);for (int i=0;i<lim;i++) g[i]=(ll)g[i]*h[i]%MOD;NTT(g,1);for (int i=0;i<m;i++) A[i]=dec(A[i],g[i]);for (int i=m;i<lim;i++) A[i]=0; } int n,k; inline void mul(int* A,int* B,int* C,int* M) {static int f[MAXN],g[MAXN];for (int i=0;i<k;i++) f[i]=A[i],g[i]=B[i];for (l=0;(1<<l)<(k<<1);l++);init();for (int i=k;i<lim;i++) f[i]=g[i]=0;NTT(f,0);NTT(g,0);for (int i=0;i<lim;i++) C[i]=(ll)f[i]*g[i]%MOD;NTT(C,1);getmod(C,M,(k-1)<<1,k); } inline void qpow(int* A,int* R,int p) {static int ans[MAXN];memset(ans,0,sizeof(ans));ans[0]=1;while (p){if (p&1) mul(ans,A,ans,R);mul(A,A,A,R);p>>=1;}for (int i=0;i<k;i++) A[i]=ans[i]; } int f[MAXN],a[MAXN],h[MAXN]; int main() {rt[0][23]=qpow(3,119);rt[1][23]=inv(rt[0][23]);for (int i=23;i>=1;i--) rt[0][i-1]=(ll)rt[0][i]*rt[0][i]%MOD,rt[1][i-1]=(ll)rt[1][i]*rt[1][i]%MOD;n=read(),k=read();for (int i=1;i<=k;i++) (f[k-i]=-read())<0&&(f[k-i]+=MOD);for (int i=0;i<k;i++) (a[i]=read())<0&&(a[i]+=MOD);f[k]=1;h[1]=1;qpow(h,f,n);int ans=0;for (int i=0;i<k;i++) ans=add(ans,(ll)h[i]*a[i]%MOD);printf("%d\n",ans);return 0; }

總結

以上是生活随笔為你收集整理的常系数齐次线性递推学习笔记的全部內容，希望文章能夠幫你解決所遇到的問題。

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