常系数齐次线性递推

发表于 更新于 本文字数： 4.7k 阅读时长 ≈ 4 分钟
常系数齐次线性递推

建议学习：> here <

一个线性变换会使一个向量在方向上发生偏移，但是如果能找到变换后方向不会发生偏移的向量，将其当作基向量，就可以快速计算递推式的任意项。
将这种基向量称作特征向量 $ \vec{v} $ ，每次变换后伸长或缩短的倍数称作特征值 $ \lambda $ ，转移矩阵为$ A $，其中 $ A $ 有 $ n $ 列（转移长度为 $ n $）。
于是：
$$
A \vec{v} = \lambda \vec{v}\\
(\lambda I - A) \vec{v} = 0
$$

其中 $ \vec{v} $ 取零向量是无意义的。
若要使 $ \vec{v} $ 为非零解，则 $ \det(\lambda I - A) = 0 $ ，即将空间降维。
其中 $ \det(\lambda I - A) = 0 $ 是次数为 $ A $ 的列数的特征多项式 $ f(\lambda) $ 。
根据Cayley-Hamilton定理， $ f(A) = 0 $ ，证明莫得感兴趣的话百度一下。
转移矩阵大概是长这样的：
$$
\begin{bmatrix}
a_1&a_2&a_3&\dots&a_{n-1}&a_n\\
1 &0 &0 &\dots&0&0 \\
0 &1 &0 &\dots&0&0 \\
\dots&\dots&\dots&\dots&\dots&\dots\\
0 &0 &0 &\dots&1&0
\end{bmatrix}
$$

特征多项式是长这样的：
$$
\det(\begin{bmatrix}
\lambda-a_1&-a_2&-a_3&\dots&-a_{n-1}&-a_n\\
-1 &\lambda &0 &\dots&0&0 \\
0 &-1 &\lambda &\dots&0&0 \\
\dots&\dots&\dots&\dots&\dots&\dots\\
0 &0 &0 &\dots&-1&\lambda
\end{bmatrix})
$$

这个东西的行列式可以手算（判断第一行选哪一个）：
$$
f(\lambda) = \lambda^n - \sum_{i = 1}^{n}a_i\lambda^{n-i}
$$

设初始项 $ H $，求第 $ m $ 项。
求 $ (A^{m} \times H)_0 $ 。

$$
\because f(A) = 0\\
\therefore A^m \bmod f(A) = A^m
$$

直接多项式快速幂 + 多项式取模计算出 $ A^m \bmod f(A) $ 。
设之后得到的多项式
$$
g(A) = \sum_{i = 0}^{n - 1} c_i A^i
$$

最后求
$$
\sum_{i = 0}^{n - 1} (c_{i} A^{i} H){0} = \sum{i = 0}^{n - 1} c_{i} h_{i}
$$

授之以渔不如授之以鱼 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
#include <bits/stdc++.h>
#define rep(i, a, b) for(int i = (a); i <= (b); i++)
#define per(i, a, b) for(int i = (a); i >= (b); i--)
#define pb push_back
using namespace std;
typedef long long ll;
const int mod = 998244353;

inline int Max(const int &x, const int &y) { return x > y ? x : y; }
inline int Min(const int &x, const int &y) { return x < y ? x : y; }
inline int add(const int &x, const int &y) { return x + y < mod ? x + y : x + y - mod; }
inline int sub(const int &x, const int &y) { return x - y < 0 ? x - y + mod : x - y; }
inline int mul(const int &x, const int &y) { return (int)((ll)x * y % mod); }
inline int ksm(int x, int y = mod - 2) {
	int ss = 1; for(; y; y >>= 1, x = mul(x, x)) if(y & 1) ss = mul(ss, x); return ss;
}
namespace Poly {
	inline int Get(int x) { int ss = 1; for(; ss <= x; ss <<= 1); return ss; }
	void ntt(vector<int> &A, int lmt, int opt) {
		A.resize(lmt + 5);
		for(int i = 0, j = 0; i < lmt; i++) {
			if(i > j) swap(A[i], A[j]);
			for(int k = lmt >> 1; (j ^= k) < k; k >>= 1);
		}
		vector<int> w(lmt >> 1);
		for(int mid = 1; mid < lmt; mid <<= 1) {
			w[0] = 1;
			int w0 = ksm(opt == 1 ? 3 : (mod + 1) / 3, (mod - 1) / mid / 2);
			for(int j = 1; j < mid; j++) w[j] = mul(w[j - 1], w0);
			for(int R = mid << 1, j = 0; j < lmt; j += R)
				for(int k = 0; k < mid; k++) {
					int x = A[j + k], y = mul(w[k], A[j + mid + k]);
					A[j + k] = add(x, y), A[j + mid + k] = sub(x, y);
				}
		}
		if(opt == -1)
			for(int i = 0, inv = ksm(lmt); i < lmt; i++) A[i] = mul(A[i], inv);
	}
	vector<int> Mul(const vector<int> &a, const vector<int> &b) {
		vector<int> A = a, B = b; int lmt = Get(a.size() + b.size() - 2);
		ntt(A, lmt, 1), ntt(B, lmt, 1);
		for(int i = 0; i < lmt; i++) A[i] = mul(A[i], B[i]);
		ntt(A, lmt, -1); return A.resize(a.size() + b.size() - 1), A;
	}
	vector<int> Inv(const vector<int> &A, int sz = -1) {
		if(sz == -1) sz = A.size();
		vector<int> res; if(sz == 1) return res.pb(ksm(A[0])), res;
		res = Inv(A, (sz + 1) / 2);
		vector<int> tmp(A.begin(), A.begin() + sz);
		int lmt = Get(sz * 2 - 2);
		ntt(tmp, lmt, 1), ntt(res, lmt, 1);
		for(int i = 0; i < lmt; i++) res[i] = mul(sub(2, mul(res[i], tmp[i])), res[i]);
		ntt(res, lmt, -1); return res.resize(sz), res;
	}
	void Div(const vector<int> &A, const vector<int> &B, vector<int> &D, vector<int> &R) {
		if(B.size() > A.size()) return (void)(D.clear(), D.pb(0), R = A);
		vector<int> a = A, b = B, iB; int n = A.size(), m = B.size();
		reverse(a.begin(), a.end()), reverse(b.begin(), b.end());
		b.resize(n - m + 1), iB = Inv(b, n - m + 1);
		D = Mul(a, iB), D.resize(n - m + 1), reverse(D.begin(), D.end());
		R = Mul(B, D);
		for(int i = 0; i < m; i++) R[i] = (mod + A[i] - R[i]) % mod;
		R.resize(m);
	}
}
int n, m, ans;
vector<int> f, h, res, b, tmp;

vector<int> Mul(const vector<int> &a, const vector<int> &b) {
	vector<int> t = Poly::Mul(a, b), r;
	Poly::Div(t, f, tmp, r); return r;
}
int main() {
	scanf("%d%d", &n, &m);
	if(n <= m) {
		rep(i, 1, m) {
			int x; scanf("%d", &x);
			if(i == n) { printf("%d\n", (x % mod + mod) % mod); return 0; }
		}
		return 0;
	}
	f.resize(m + 1), h.resize(m), f[m] = 1;
	rep(i, 1, m) {
		scanf("%d", &f[m - i]);
		f[m - i] = ((mod - f[m - i]) % mod + mod) % mod;
	}
	rep(i, 0, m - 1) scanf("%d", &h[i]), h[i] = (h[i] % mod + mod) % mod;
	res.pb(1), b.pb(0), b.pb(1);
	for(; n; n >>= 1, b = Mul(b, b)) if(n & 1) res = Mul(res, b);
	rep(i, 0, m - 1) ans = (ans + 1ll * res[i] * h[i]) % mod;
	printf("%d\n", ans);
	return 0;
}