本博客部分参考自dtz的博客，Troywar的博客
前言：
为什么要来学这个鬼畜的东西呢……因为模拟赛出了一道板题（好吧这实际上就是BZOJ4162）：
看到的时候我非常兴奋，这不是一道sb矩阵快速幂题吗？连邻接矩阵都直接给出来了；
结果一看数据范围，算了算$n^3logk=1.25e10$，emmm我咋不知道矩阵快速幂有什么优化方法？
于是我就自闭了，写了50滚粗（然后一个点1008ms一个点1114ms少了20，气傻了）
然后才知道这是个常系数齐次线性递推模板题……
常系数齐次线性递推
0.预备知识：
有一个数列递推式为：
$$S_i=sumlimits_{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|=egin{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&lambdaend{bmatrix}-egin{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_1end{bmatrix}$$
$$=egin{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_1end{bmatrix}$$
把$phi(lambda)$按最后一列$(j=k)$拉普拉斯展开得
$$phi(lambda)=sumlimits_{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-sumlimits_{i=1}^{k}a_ilambda^{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)=sumlimits_{i=0}^{k-1}b_ix^i$，则
$$r(A)=sumlimits_{i=0}^{k-1}b_iA^i$$
$$A^n=sumlimits_{i=0}^{k-1}b_iA^i$$
设$T_i$为表示数列中连续$k$项的向量$egin{bmatrix}S_i&S_{i+1}&cdots&S_{i+k-1}end{bmatrix}$（就是矩阵乘法每次乘出来的那个列向量），把$n-k+1$代回去得：
$$A^{n-k+1}T_0=T_0sumlimits_{i=0}^{k-1}b_iA^i=sumlimits_{i=0}^{k-1}b_iA^iT_i=sumlimits_{i=0}^{k-1}b_iT^i$$
答案就是$A^{n-k+1}T_0$的最后一位，即：
$$S_n=sumlimits_{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
",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 }