常系数齐次线性递推
具体记在笔记本上了，以后可能补照片，这里稍微写一下，主要贴代码。
概述
形式：
[h_n = a_1 h_{n-1}+a_2h_{n-2}+...+a_kh_{n-k}
]
矩阵乘法是(O(k^3 log n))
利用特征多项式可以做到(O(k^2log n))
特征多项式
[f(lambda) = mid M - lambda Imid
]
是关于(lambda)的(n)次多项式
根据(Cayley-hamilton)定理得到 (f(M)=0)(zero matrix)，所以就能得到(M^k)与(M^0...M^{k-1})的线性递推关系，(M^n)也可以用(M^0...M^{k-1})线性表示，多项式乘法快速幂求这个递推关系的系数。最后代入就行了。
其实就是：
[M^n = M^n mod f(M)
]
所以还需要多项式求余。
这玩意卡我好长时间，然后发现就是手算的原理。当然用毒瘤算法可以做到nlogn
两种形式
对于常系数齐次线性递推关系，我们可以一眼看出它的特征多项式
[f(lambda) = lambda^k - a_1 lambda^{k-1} -...-a_k
]
并且最高次系数为1，求余很简单
对于一般的矩阵，需要求行列式得到n+1个点的值，然后拉格朗日插值，这一步复杂度(O(n^4))
拉格朗日插值公式：
[sum_{j=0}^n y_j prod_{i=0 , i
eq j}^n frac{x-x_i}{x_j - x_i}
]
代码
例题为bzoj 4161 & bzoj 4162
常系数齐次线性递推
#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <queue>
using namespace std;
typedef long long ll;
const int N = 4005, mo = 1e9+7;
inline int read(){
char c=getchar(); int x=0,f=1;
while(c<'0'||c>'9') {if(c=='-')f=-1;c=getchar();}
while(c>='0'&&c<='9') {x=x*10+c-'0';c=getchar();}
return x*f;
}
int n, k, a[N], f[N], b[N], c[N], ans;
void mul_mod(int *x, int *y) {
static int c[N];
for(int i=0; i<k<<1; i++) c[i] = 0;
for(int i=0; i<k; i++)
for(int j=0; j<k; j++) c[i+j] = (c[i+j] + (ll) x[i] * y[j]) %mo;
// mod f(M)
for(int i=2*k-2; i>=k; i--)
for(int j=1; j<=k; j++) c[i-j] = (c[i-j] + (ll) c[i] * a[j]) %mo;
for(int i=0; i<k; i++) x[i] = c[i];
}
void Pow(int *a, int b, int *ans) {
for(; b; b>>=1, mul_mod(a, a)) if(b&1) mul_mod(ans, a);
}
int main() {
freopen("in", "r", stdin);
n=read(); k=read();
for(int i=1; i<=k; i++) a[i] = read(), a[i] += a[i] < 0 ? mo : 0;
for(int i=0; i<k; i++) f[i] = read(), f[i] += f[i] < 0 ? mo : 0;
c[1] = 1; b[0] = 1;
Pow(c, n, b);
for(int i=0; i<k; i++) ans = (ans + (ll) b[i] * f[i]) %mo;
printf("%d
", ans);
}
一般矩阵快速幂
#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
typedef long long ll;
#define pll pair<ll, ll>
#define fir first
#define sec second
const int N = 52, mo = 1e9+7;
inline int read(){
char c=getchar(); int x=0,f=1;
while(c<'0'||c>'9') {if(c=='-')f=-1;c=getchar();}
while(c>='0'&&c<='9') {x=x*10+c-'0';c=getchar();}
return x*f;
}
ll Pow(ll a, int b) {
ll ans = 1;
for(; b; b>>=1, a=a*a%mo)
if(b&1) ans=ans*a%mo;
return ans;
}
int n, k, a[N][N], t[N][N], t2[N][N], ans[N][N];
char s[10005];
void mul(int a[N][N], int b[N][N], int c[N][N]) {
for(int i=1; i<=n; i++)
for(int k=1; k<=n; k++) if(a[i][k])
for(int j=1; j<=n; j++) if(b[k][j])
c[i][j] = (c[i][j] + (ll) a[i][k] * b[k][j]) %mo;
}
int det(int a[N][N]) {
int s = 0;
for(int i=1; i<=n; i++) {
int r;
for(r=i; r<=n; r++) if(a[r][i]) break;
if(r == n+1) return 0;
if(r != i) {
s ^= 1;
for(int j=1; j<=n; j++) swap(a[r][j], a[i][j]);
}
ll inv = Pow(a[i][i], mo-2);
for(int k=i+1; k<=n; k++) {
ll t = (ll) (mo - a[k][i]) * inv %mo;
for(int j=i; j<=n; j++) a[k][j] = (a[k][j] + t * a[i][j]) %mo;
}
}
ll ans = 1;
for(int i=1; i<=n; i++) ans = ans * a[i][i] %mo;
return (s&1) ? mo - ans : ans;
}
struct poly {
int a[N<<1];
poly(int x=0) {memset(a, 0, sizeof(a)); a[0] = x;}
int& operator [](int x) {return a[x];}
poly operator + (poly b) {
poly c;
for(int i=0; i<=n; i++) c[i] = (a[i] + b[i]) %mo;
return c;
}
poly operator * (ll b) {
poly c;
for(int i=0; i<=n; i++) c[i] = (ll) a[i] * b %mo;
return c;
}
poly operator * (poly b) {
poly c;
for(int i=0; i<=n; i++)
for(int j=0; j<=n; j++) c[i+j] = (c[i+j] + (ll) a[i] * b[j]) %mo;
return c;
}
poly operator * (pll t) {
ll k = t.fir, b = t.sec;
poly c(a[0] * b %mo);
for(int i=1; i<=n; i++) c[i] = (a[i-1] * k + a[i] * b) %mo;
return c;
}
friend poly operator % (poly a, poly b) {
for(int i=n; i>=0; i--) {
ll t = (ll) (mo - a[i+n]) * Pow(b[n], mo-2) %mo;
for(int j=0; j<=n; j++) a[i+j] = (a[i+j] + b[j] * t) %mo;
}
return a;
}
} ;
int y[N];
int main() {
freopen("in", "r", stdin);
scanf("%s",s);
n=read();
for(int i=1; i<=n; i++)
for(int j=1; j<=n; j++) a[i][j] = read();
for(int x=0; x<=n; x++) {
memcpy(t, a, sizeof(a));
for(int i=1; i<=n; i++) t[i][i] = (t[i][i] + mo - x) %mo;
y[x] = det(t);
}
poly f;
for(int j=0; j<=n; j++) {
poly t(1);
for(int i=0; i<=n; i++) if(i != j)
t = (t * make_pair(1, mo-i) ) * Pow(j-i+mo, mo-2);
t = t * y[j];
f = f + t;
}
poly p, b(1); p[1] = 1;
for(int i=strlen(s)-1; i>=0; i--, p = p * p % f)
if(s[i] == '1') b = b * p % f;
memset(t, 0, sizeof(t));
for(int i=1; i<=n; i++) t[i][i] = 1;
for(int p=0; p<n; p++) {
for(int i=1; i<=n; i++)
for(int j=1; j<=n; j++)
ans[i][j] = (ans[i][j] + (ll) t[i][j] * b[p]) %mo;
memset(t2, 0, sizeof(t2));
mul(t, a, t2);
memcpy(t, t2, sizeof(t2));
}
for(int i=1; i<=n; i++)
for(int j=1; j<=n; j++) printf("%d%c", ans[i][j], "
"[j==n]);
return 0;
}