常系数齐次线性递推

$\boldsymbol{k}$阶常系数齐次线性递推指的是以下问题:

给出$k$阶线性递推式$\boldsymbol{f}$以及初始项$\boldsymbol{a_0,a_1…a_{k-1}}$,并满足:
$$
a_{n}=\sum_{i=1}^{k} f_{i} \times a_{n-i}
$$
对于某个$o$,求解$f_o$在某个固定模数下的值。

然后以下是shadowice1984的人话版本:

1.它是常系数的,换句话说递推系数和下标n无关

2.它是线性的,换句话说递推式中每一项的次数都是1,没有乱七八糟的2次3次项

3.它是齐次的,换句话说递推式的常数项等于0

这个问题的弱智解法是$O(k^3\log n)$的矩阵快速幂。但是通过多项式加速可以将复杂度优化至$O(k^2\log n)$甚至$O(k\log k \log n)$.

于是这个世界上就再也没有矩阵快速幂了

$1$ 特征向量与特征多项式

首先,给定$n$阶矩阵$\boldsymbol{A}$,若存在向量$\boldsymbol{x}$使得$\exists \lambda \in \mathbb {R}$ ,满足:
$$
\tag{1} \boldsymbol{A}\cdot \boldsymbol{x} = \lambda\boldsymbol{x}
$$
则称$\boldsymbol{x}$是$\boldsymbol{A}$的特征向量,$\lambda$称为$\boldsymbol{A}$的特征值,$(1)$式称之为特征方程

同时也可以写作
$$
(\boldsymbol{A}-\lambda\boldsymbol{I})\cdot \boldsymbol{x}= \boldsymbol{0}
$$
有解的条件是
$$
\tag{2} \det(\boldsymbol{A}-\lambda\boldsymbol{I})=0
$$
同时定义特征多项式类形式幂级数,即不考虑未知数$x$的数学形态,$x$可以为矩阵、向量$~etc.$

我们观察,$(2)$中右边的行列式计算的性质保证了行列式展开后每一项一定至多是$\lambda$的$n$次(主对角线相乘)、$n-1$次$\cdots$ $0$次,那么特征多项式即为关于$\lambda$的$\det(\boldsymbol{A}-\lambda\boldsymbol{I})$展开,记为$\boldsymbol{\phi(\lambda)}$.

接着观察$\phi(\lambda)$:
$$
\begin{aligned}
\phi(\lambda)= |\boldsymbol{A}-\lambda \boldsymbol{I}| & = \left[\begin{array}{ccccc}{~\lambda} & {0} & {\cdots} & {0} & {0} \\ {0} & {\lambda} & {\cdots} & {0} & {0} \\ {\cdots} & {\cdots} & {\cdots} & {\cdots} & {\cdots} \\ {0} & {0} & {\cdots} & {\lambda} & {0} \\ {0} & {0} & {\cdots} & {0} & {\lambda}\end{array}\right]-\left[\begin{array}{ccccc}{~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} & {0} & {\lambda}\end{array}\right]\\ & =\left[\begin{array}{ccccc}{\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_{1}}\end{array}\right]
\end{aligned}
$$

然后对于这个我们可以进行拉普拉斯展开(就是用代数余子式展开一个行列式),我们选择最后一列:

$$
\phi(\lambda)=\sum\limits_{i=1}^{k} a_{k-i+1}(-1)^{i+k}\phi(\lambda)_{i, k}
$$

其中$\phi(\lambda)_{i,k}$表示丢掉$(i,k)$同行同列之后的代数余子式。

这东西可以暴力迭代式化简/直接良序归纳得出他的背诵化简形式:
$$
\phi(\lambda)=(-1)^{n}\left(\lambda^{n}-\sum_{i=1}^{n} a_{i} \lambda^{n-i}\right)
$$

好像胡扯了半天也没扯到点子上,递推在哪??

$2$ 一步转化/构造

現在は大力构造時間です!!

我们考虑构造一个这样的数列$\boldsymbol{g}$,满足:
$$
\boldsymbol{A}^n=\sum_{i=0}^{k-1}g_i\boldsymbol{A}^i
$$

然后我们两边同时左乘一个初始状态$\boldsymbol{F_0}$:
$$
\boldsymbol{F_0\cdot A}^n=\sum_{i=0}^{k-1}g_i\boldsymbol{F_0\cdot A}^i
$$

我们考虑对于矩阵乘法的$\boldsymbol{res}$矩阵,我们关心的是它的$(0,0)$处的值,不妨记为$\boldsymbol{res_0}$,那么就有:
$$
\boldsymbol{res}_0=(\boldsymbol{F_0\cdot A} ^ n)_0 =\sum_{i=0}^{k-1}g_i(\boldsymbol{F_0\cdot A}^i)_0
$$
其中$\boldsymbol{X}_p$表示一维向量$\boldsymbol{X}$的第$p$项。

然后使用观察法,我们发现$\boldsymbol{F_0\cdot A}^i$这东西就是转移了$i$次之后的第$0$项,换言之也就是$\boldsymbol{F_0}$的第$i$项,那么就会有:
$$
\boldsymbol{res_0}=\sum_{i=0}^{k-1}g_i\boldsymbol{(F_0)}_i
$$
所以只要我们可以构造出一组$g$,我们就可以在$O(k)$的时间内求出答案。

$3~$ 快速线性递推

现在我们的任务就是构造$\boldsymbol{g}$序列了。

首先我们随便设一个关于$\boldsymbol{A}$的多项式:
$$
\boldsymbol{A}^n=P(\boldsymbol{A})Q(\boldsymbol{A})+R(\boldsymbol{A})
$$
于是我们发现,似乎根据带余除法的性质$\deg (R(\boldsymbol{A}))\leq \deg(\boldsymbol{A}^n)=n$,正好满足$\deg(\sum_{i=0}^{k-1}g_i\boldsymbol{(F_0)}_i)\leq n$的性质,所以我们考虑当$P(\boldsymbol{A})Q(\boldsymbol{A})$这一项为$\boldsymbol{0}$时,恰好构造出一组可行的答案。

我们结合特征多项式考虑,$\phi(\lambda)=0$,恰好满足我们的需求。于是我们直接令$\boldsymbol{A}^n \bmod \phi(\lambda)$即可求出$R(\boldsymbol{A})=g$。

然后普通的就可以直接$k^2$暴力NTT,如果很闲并且模数是NTT模数的话还可以用$\log^2$的复杂度搞出来。

以下是LuoguP4723

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
#define Gp 3
#define rr register
#define MAXN 251003
#define LL long long
#define Mod 998244353

using namespace std ;

LL base[MAXN], p[MAXN] ;
int N, M, L, K, R[MAXN] ;
LL F[MAXN], G[MAXN], P[MAXN], Ans ;

inline LL expow(LL a, int b){
rr LL res = 1 ;
while (b){
if (b & 1) (res *= a) %= Mod ;
a = 1ll * a * a % Mod, b >>= 1 ;
}
return res ;
}
void NTT(LL *J, const int &l, const int &flag){
rr LL Gn, Gi = 1 ;
for (rr int i = 0 ; i < l ; ++ i)
if (i < R[i]) J[i] ^= J[R[i]] ^= J[i] ^= J[R[i]] ;
for (rr int i = 1 ; i < l ; i <<= 1){
Gn = expow(Gp, (Mod - 1) / (i << 1)) ;
for (int j = 0 ; j < l ; Gi = 1, j += (i << 1)){
for (int k = 0 ; k < i ; ++ k, Gi = 1ll * Gi * Gn % Mod){
int real = J[j + k], iroot = 1ll * J[j + k + i] * Gi % Mod ;
J[j + k] = (real + iroot) % Mod, J[j + k + i] = (real - iroot + Mod) % Mod ;
}
}
}
if (flag > 0) return ;
rr int Inv = expow(1ll * l, Mod - 2) ; reverse(J + 1, J + l) ;
for (rr int i = 0 ; i < l ; ++ i) J[i] = 1ll * J[i] * Inv % Mod ;
}
LL t[MAXN], T[MAXN] ;
void _Inv(LL *f, LL *g, int len){
if (len <= 1) {
g[0] = expow(f[0], Mod - 2) ;
return ;
}
rr int i, l = 0, Len = 1 ;
_Inv(f, g, (len + 1) >> 1) ;
while (Len < (len << 1)) Len <<= 1, ++ l ;
for (i = 0 ; i < Len ; ++ i) R[i] = (R[i >> 1] >> 1) | ((i & 1) << (l - 1)) ;
for (i = 0 ; i < len ; ++ i) t[i] = f[i] ; NTT(g, Len, 1), NTT(t, Len, 1) ;
for (i = 0 ; i < Len ; ++ i) g[i] = (2ll - t[i] * g[i] % Mod + Mod) % Mod * g[i] % Mod ;
NTT(g, Len, -1) ; for (i = len ; i < Len ; ++ i) g[i] = 0 ;
}
LL IG[MAXN], Ft[MAXN], Gt[MAXN] ;
void _Div(LL *f, const int &Len){
rr int L1 = (K << 1), D, i ;
while (!f[-- L1]) ; if (L1 < K) return ;
for (i = 0 ; i < Len ; ++ i) Ft[i] = 0 ;
for (i = 0 ; i <= L1 ; ++ i) Ft[i] = f[i] ; D = L1 - K + 1 ;
reverse(Ft, Ft + L1 + 1) ; for (i = D ; i <= L1 ; ++ i) Ft[i] = 0 ;
NTT(Ft, Len, 1) ; for (i = 0 ; i < Len ; ++ i) P[i] = Ft[i] * IG[i] % Mod ;
NTT(P, Len, - 1) ; for (i = D ; i <= L1 ; ++ i) P[i] = 0 ; reverse(P, P + D) ;
NTT(P, Len, 1) ; for (i = 0 ; i < Len ; ++ i) (P[i] *= G[i]) %= Mod ; NTT(P, Len, -1) ;
for (i = 0 ; i < K ; ++ i) f[i] = (f[i] - P[i] + Mod) % Mod ; for (i = K ; i <= L1 ; ++ i) f[i] = 0 ;
}
int main(){
// freopen("2.in", "r", stdin) ;
rr int Len = 1, l = 0, i ; cin >> N >> K ;
while (Len < (K << 1)) Len <<= 1, ++ l ; F[0] = 1 ; Len <<= 1, ++ l ;
for (i = 0 ; i < Len ; ++ i) R[i] = (R[i >> 1] >> 1) | ((i & 1) << (l - 1)) ;
for (i = 1 ; i <= K ; ++ i) scanf("%lld", &p[i]), p[i] = p[i] < 0 ? p[i] + Mod : p[i] ;
for (i = 0 ; i < K ; ++ i) scanf("%lld", &base[i]), base[i] = base[i] < 0 ? base[i] + Mod : base[i] ;
for (i = 1 ; i <= K ; ++ i) Gt[K - i] = G[K - i] = Mod - p[i] ; T[1] = G[K] = Gt[K] = 1 ;reverse(Gt, Gt + K + 1) ;
NTT(G, Len, 1) ; _Inv(Gt, IG, Len >> 1), NTT(IG, Len, 1) ;
while(N){
NTT(T, Len, 1) ;
if (N & 1) {
NTT(F, Len, 1) ;
for (i = 0 ; i < Len ; ++ i) F[i] = F[i] * T[i] % Mod ;
NTT(F, Len, -1) ; _Div(F, Len) ;
}
for (i = 0 ; i < Len ; ++ i) (T[i] *= T[i]) %= Mod ; NTT(T, Len, -1) ; _Div(T, Len), N >>= 1 ;
}
for (i = 0 ; i < K ; ++ i) (Ans += (F[i] * base[i]) % Mod) %= Mod ; printf("%lld ", Ans) ; return 0 ;
}