【模板】线性递推

求 $f_n,f_n={\sum }_{i=1}^k{a}_i{f}_{n-i}$
$k\le 3e4,n\le 1e9$
首先暴力的矩阵乘是 $k^{3}logn$的
引入概念：特征多项式
$f(\lambda )=|\lambda I-A|$
$I$ 为单位矩阵
然后对于这道题来讲，转移矩阵是这样的
$\lambda I-A=$
$\begin{array}{ccccc}\lambda -a_1& -a_2& ...& -a_{k-1}& -a_k\\ -1& \lambda & 0& ...& 0\\ 0& -1& \lambda & ...& 0\\ & & ...& & \\ 0& 0& ...& -1& \lambda \end{array}$
对第一行展开然后分别求行列式
$f(\lambda )=(\lambda -a_1)\ast A_{1,1}+(-a_2)\ast A_{1,2}+...+(-a_k)\ast A_{1,k}$
然后发现去除第一行和某一列后主对角线都是满的，消成下三角矩阵，行列式是主对角线的乘积
$f(\lambda )={\lambda }^k-{\sum }_{i=0}^{k-1}{a}_{k-i}{\lambda }^i$
然后由 $Cayley-Hamilton$ 定理
$f(A)=0$，一个巧妙的伪证是 $f(A)=|AI-A|=0$
设 $R(A)=A^n\%f(A)$，即 $g(A)\ast f(A)+R(A)=A^n$
又因为 $f(A)=0$ ，所以 $A^n=R(A)$
看看求的是什么，是向量 $v$，乘上 $A^n$ 的第 k 列，令$f(A)$ 的各项系数为 $c$
即 $(v\ast A^n{)}_k=(v\ast {\sum }_{i=0}^k{c}_i{A}^i{)}_k=({\sum }_{i=0}^k{c}_i\ast (v\ast A^i){)}_k$
初始向量 $v$ 乘上转移矩阵 $i$ 次，即为给出的 $f_i$
于是 $ans={\sum }_{i=0}^k{c}_i\ast f_i$
$c$即特征多项式的系数可以用 $A^n$对$f(\lambda )={\lambda }^k-{\sum }_{i=0}^{k-1}{a}_{k-i}{\lambda }^i$ 多项式取模得出
复杂度 $O(k\ast log(k)\ast log(n))$
主要的思想就是不必求出 $A^n$，将$A^n$ 转换成${\sum }_{i=0}^k{c}_i{A}^i$
并且也不求 $A^i$，而是用初始向量去乘它${\sum }_{i=0}^k{c}_i(A_i\ast v)$，得到了我们想要的答案

#include<bits/stdc++.h>
#define N 400050
using namespace std;
typedef long long ll;
const int Mod = 998244353, G = 3;
ll add(ll a, ll b){ return a + b >= Mod ? a + b - Mod : a + b;}
ll mul(ll a, ll b){ return a * b % Mod;}
ll power(ll a, ll b){ ll ans = 1;
	for(;b;b>>=1){if(b&1) ans = mul(ans, a); a = mul(a, a);}
	return ans;
}
ll inv[N]; int n, k;
#define poly vector<ll>
#define C 20
poly w[C + 1];
int up, bit, rev[N]; ll a[N];
void Init(int len){ up = 1, bit = 0;	
	while(up < len) up <<= 1, bit++;
	for(int i = 0; i < up; i++) rev[i] = (rev[i>>1]>>1)|((i&1)<<(bit-1));
}
void prework(){
	inv[0] = inv[1] = 1;
	for(int i = 2; i <= N-50; i++) inv[i] = mul(Mod-Mod/i, inv[Mod%i]);
	for(int i = 1; i <= C; i++) w[i].resize(1<<(i-1));
	ll wn = power(G, (Mod-1)/(1<<C)); w[C][0] = 1; 
	for(int i = 1; i < (1<<(C-1)); i++) w[C][i] = mul(w[C][i-1], wn);
	for(int i = C-1; i; i--)
		for(int j = 0; j < (1<<(i-1)); j++)
			w[i][j] = w[i+1][j<<1];
}
void NTT(poly &a, int flag){
	for(int i = 0; i < up; i++) if(i < rev[i]) swap(a[i], a[rev[i]]);
	for(int i = 1, l = 1; i < up; i <<= 1, l++)
		for(int j = 0; j < up; j += (i<<1))
			for(int k = 0; k < i; k++){
				ll x = a[k + j], y = mul(w[l][k], a[k + j + i]);
				a[k + j] = add(x, y); a[k + j + i] = add(x, Mod-y);
			}
	if(flag == -1){
		reverse(a.begin() + 1, a.begin() + up);
		for(int i = 0; i < up; i++) a[i] = mul(a[i], inv[up]);
	}
}
poly operator * (poly a, poly b){
	int len = a.size() + b.size() - 1;
	Init(len); a.resize(up); b.resize(up);
	NTT(a, 1); NTT(b, 1);
	for(int i = 0; i < up; i++) a[i] = mul(a[i], b[i]); 
	NTT(a, -1); a.resize(len); return a;
}
poly Inv(poly a, int len){
	poly b(1, power(a[0], Mod - 2)), c;
	for(int lim = 4; lim < (len << 2); lim <<= 1){
		Init(lim); c = a; c.resize(lim >> 1);
		c.resize(up); b.resize(up); NTT(c, 1); NTT(b, 1);
		for(int i = 0; i < up; i++) b[i] = mul(b[i], add(2, Mod - mul(b[i], c[i])));
		NTT(b, -1); b.resize(lim >> 1);
	} b.resize(len); return b;
}
poly operator - (poly a, poly b){
	poly c; int len = max(a.size(), b.size()); c.resize(len);
	a.resize(len); b.resize(len);
	for(int i = 0; i < len; i++) c[i] = add(a[i], Mod - b[i]);
	return c;
}
poly operator / (poly a, poly b){
	int lim = 1, len = a.size() - b.size() + 1;
	reverse(a.begin(), a.end()); reverse(b.begin(), b.end());
	while(lim <= len) lim <<= 1;
	b = Inv(b, lim); b.resize(len);
	a = a * b; a.resize(len);
	reverse(a.begin(), a.end()); 
	return a;
}
poly operator % (poly a, poly b){ 
	int len = a.size() - b.size() + 1;
	if(len < 0) return a;
	poly c = a - (a / b) * b;
	c.resize(b.size() - 1); return c;
}
int main(){
	scanf("%d%d", &n, &k); prework();
	poly f(k + 1);
	for(int i = 1; i <= k; i++) scanf("%lld", &f[k - i]), f[k - i] = Mod - (f[k - i] % Mod + Mod) % Mod;
	f[k] = 1;
	for(int i = 0; i < k; i++) scanf("%lld", &a[i]), a[i] = (a[i] % Mod + Mod) % Mod;
	poly res, g;
	res.resize(k + 1); res[0] = 1;
	g.resize(k + 1); g[1] = 1;
	for(;n;n>>=1, g = g * g % f) if(n & 1) res = res * g % f;
	ll ans = 0;
	for(int i = 0; i < res.size(); i++) ans = add(ans, mul(res[i], a[i])); 
	cout << ans; 
	return 0;
}