特征多项式和常系数线性齐次递推学习笔记

特征多项式

\(A\) 为给定的 \(n\times n\) 矩阵,\(I_n\)\(n\times n\) 单位矩阵,\(A\) 的特征多项式定义为

\[ p(\lambda )=\det(\lambda I_{n}-A) \]

其中 \(\det\) 表示行列式

Cayley–Hamilton theorem

根据 凯莱–哈密顿定理\(A\) 满足方程

\[ p(A)=O \]

其中 \(O\) 是零矩阵

因此我们可以利用这个 \(n\) 次的多项式 \(p(A)\) 来降低 \(A\) 的高次幂

求特征多项式

代入 \(k+1\) 个值,求行列式,插值得到多项式系数

高斯消元求行列式是 \(\mathcal O(n^3)\),总复杂度 \(\mathcal O(n^4)\)

\(\mathcal O(n^3)\) 的做法我就不学了

这里有

求矩阵的高次幂

BZOJ 4162: shlw loves matrix II

\(k\times k\) 的矩阵 \(M\)\(n\) 次幂

\(k\le 50, n\le 2^{10000}\)

暴力

直接暴力不卡常本机 2s,复杂度\(\mathcal O(k^3\log n)\)

正常做法

\(\mathcal O(k^4)\) 求出矩阵的特征多项式 \(p(x)\)

然后就要求出 \(x^n \bmod p(x)\),直接快速幂,其中乘法和取模都是 \(\mathcal O(k^2)\)

最后暴力乘出 \(M^0, M^1, \dotsc, M^{k-1}\),按系数计算答案,复杂度也是 \(\mathcal O(k^4)\)

总复杂度 \(\mathcal O(k^4+k^2\log n)\)

代码在例题里贴

常系数线性齐次递推

给出系数 \(c_1, c_2, \dotsc, c_k\) 和数列 \(f\) 的前 \(k\)\(f_0, f_1, \dotsc, f_{k-1}\)

对于 \(n\ge k\)\(f_n=\sum\limits_{i=1}^k f_{n-i}c_i\)

\(f_n\)

暴力

构造出转移矩阵

\[ \begin{bmatrix} c_1 & c_2 & \cdots & c_{k-1} & c_k \\ 1 & 0 & \cdots & 0 & 0 \\ 0 & 1 & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \cdots & 1 & 0 \end{bmatrix} \begin{pmatrix} f_{k-1} \\ f_{k-2} \\ f_{k-3} \\ \vdots \\ f_0 \end{pmatrix} = \begin{pmatrix} f_k \\ f_{k-1} \\ f_{k-2} \\ \vdots \\ f_1 \end{pmatrix} \]

可以做到 \(\mathcal O(k^3\log n)\)

特征多项式优化

上述转移矩阵设为 \(A\),特征多项式是

\[ \begin{align} p(\lambda) & =\det(\lambda I-A) \\ & = \det\left( \begin{bmatrix} \lambda -c_1 & -c_2 & \cdots & -c_{k-1} & -c_k \\ -1 & \lambda & \cdots & 0 & 0 \\ 0 & -1 & \cdots & \lambda & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \cdots & -1 & \lambda \end{bmatrix} \right) \end{align} \]

对第一行展开

可以发现 \(k\) 个代数余子式 \(C_{1,1},C_{1,2},\dotsc,C_{1,k}\) 分别是 \(\lambda^{k-1}, \lambda^{k-2}, \dotsc, \lambda^0\)

于是 \(p(\lambda) = \lambda^k - \sum\limits_{i=1}^k c_i\lambda^{k-i}\)

我们令 \(g(x)=x^n \bmod p(x)\),于是 \(A^n=\sum\limits_{i=0}^{k-1} g_i A^i\)

\(\vec{f}= \begin{pmatrix} f_{k-1} \\ f_{k-2} \\ f_{k-3} \\ \vdots \\ f_0 \end{pmatrix}\)\(\left[\vec{a}\right]_ k\) 表示列向量 \(\vec{a}\) 的第 \(k\)

\[ \begin{align} f(n) & =\left[A^n \vec{f}\right]_ k \\ & = \left[\sum_{i=0}^{k-1} g_i A^i \vec{f}\right]_ k \\ & = \sum_{i=0}^{k-1} g_i \left[A^i \vec{f}\right]_ k \\ & = \sum_{i=0}^{k-1} g_i f_i \end{align} \]

因为 \(A\vec{f}\) 相当于转移了一次,转移 \(i\) 次后向量的第 \(k\) 项即 \(f_i\)

实现

其实只需要求 \(g(x)=x^n \bmod p(x)\) 就好了

可以快速幂,暴力取模复杂度 \(\mathcal O(k^2\log n)\),用 NTT 可以优化到 \(\mathcal O(k\log k\log n)\)

例题

shlw loves matrix II

BZOJ 4162: shlw loves matrix II

上面讲过

代码

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
#include<cstdio>
#include<algorithm>
#include<cctype>
#include<string.h>
#include<cmath>

using namespace std;
#define ll long long

const int N = 55, M = 10005, P = 1000000007;
int k, l, f[N], g[N], d[N], s[N<<1], a[N][N], b[N][N], ans[N][N];
char n[M];
inline int Pow(ll x, int y=P-2){
int ans=1;
for(; y; y>>=1, x=x*x%P) if(y&1) ans=ans*x%P;
return ans;
}
inline int det(int x){
memcpy(b, a, sizeof a);
for(int i=0; i<k; ++i) b[i][i]=(b[i][i]+P-x)%P;
int ans=1;
for(int i=0; i<k; ++i){
if(!b[i][i]){
ans=P-ans;
for(int j=i+1; j<k; ++j) if(b[j][i]){ swap(b[i], b[j]); break;}
if(!b[i][i]) return 0;
}
ans=(ll)ans*b[i][i]%P;
int inv=Pow(b[i][i]);
for(int j=i; j<k; ++j) b[i][j]=(ll)b[i][j]*inv%P;
for(int j=i+1; j<k; ++j) for(int t=k-1; t>=i; --t)
b[j][t]=(b[j][t]+(ll)(P-b[j][i])*b[i][t])%P;
}
return ans;
}
inline void mul(int (*a)[N], int (*b)[N]){
static int tmp[N][N];
memset(tmp, 0, sizeof tmp);
for(int i=0; i<k; ++i) for(int t=0; t<k; ++t) for(int j=0; j<k; ++j)
tmp[i][j]=(tmp[i][j]+(ll)a[i][t]*b[t][j])%P;
memcpy(a, tmp, sizeof tmp);
}
int main() {
scanf("%s%d", n, &k);
for(int i=0; i<k; ++i) for(int j=0; j<k; ++j) scanf("%d", a[i]+j);
f[0]=1;
for(int i=0; i<=k; ++i) for(int j=i+1; ~j; --j) f[j]=((ll)(P-f[j])*i+(j?f[j-1]:0))%P;
for(int i=0; i<=k; ++i){
int x=1;
memcpy(g, f, (k+2)<<2);
for(int j=k+1; j; --j) g[j-1]=(g[j-1]+(ll)g[j]*i)%P;
for(int j=0; j<=k; ++j) if(i!=j) x=(ll)x*(i-j+P)%P;
x=(ll)Pow(x)*det(i)%P;
for(int j=0; j<=k; ++j) d[j]=(d[j]+(ll)g[j+1]*x)%P;
}
for(int i=0, inv=Pow(d[k]); i<=k; ++i) d[i]=(ll)d[i]*inv%P;
s[0]=1;
for(int x=0; n[x]; ++x){
static int tmp[N];
memcpy(tmp, s, k<<2), memset(s, 0, sizeof s);
for(int i=0; i<k; ++i) for(int j=0; j<k; ++j) s[i+j]=(s[i+j]+(ll)tmp[i]*tmp[j])%P;
if(n[x]&1){
for(int i=k*2-1; i; --i) s[i]=s[i-1];
s[0]=0;
}
for(int i=k*2-1; i>=k; --i)
for(int j=k; ~j; --j) s[i-j]=(s[i-j]+(ll)(P-d[k-j])*s[i])%P;
}
memset(b, 0, sizeof b);
for(int i=0; i<k; ++i) b[i][i]=1;
for(int x=0; x<k; ++x){
for(int i=0; i<k; ++i) for(int j=0; j<k; ++j)
ans[i][j]=(ans[i][j]+(ll)s[x]*b[i][j])%P;
mul(b, a);
}
for(int i=0; i<k; ++i) for(int j=0; j<k; ++j) printf("%d%c", ans[i][j], " \n"[j==k-1]);
return 0;
}

【模板】线性递推

Luogu P4723 【模板】线性递推

模板,需要 NTT 优化

代码

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
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
#include<cstdio>
#include<algorithm>
#include<cctype>
#include<string.h>
#include<cmath>
#include<vector>

using namespace std;
#define ull unsigned long long

inline char read() {
static const int IN_LEN = 1000000;
static char buf[IN_LEN], *s, *t;
return (s==t?t=(s=buf)+fread(buf,1,IN_LEN,stdin),(s==t?-1:*s++):*s++);
}
template<class T>
inline void read(T &x) {
static bool iosig;
static char c;
for (iosig=false, c=read(); !isdigit(c); c=read()) {
if (c == '-') iosig=true;
if (c == -1) return;
}
for (x=0; isdigit(c); c=read()) x=((x+(x<<2))<<1)+(c^'0');
if (iosig) x=-x;
}
const int OUT_LEN = 1000000;
char obuf[OUT_LEN], *ooh=obuf;
inline void print(char c) {
if (ooh==obuf+OUT_LEN) fwrite(obuf, 1, OUT_LEN, stdout), ooh=obuf;
*ooh++=c;
}
template<class T>
inline void print(T x) {
static int buf[30], cnt;
if (x==0) print('0');
else {
if (x<0) print('-'), x=-x;
for (cnt=0; x; x/=10) buf[++cnt]=x%10+48;
while(cnt) print((char)buf[cnt--]);
}
}
inline void flush() { fwrite(obuf, 1, ooh - obuf, stdout); }

const int N = 1<<17, P = 998244353;
struct Z{
unsigned x;
Z(const unsigned _x=0):x(_x){}
inline Z operator +(const Z &rhs)const{ return x+rhs.x<P?x+rhs.x:x+rhs.x-P;}
inline Z operator -(const Z &rhs)const{ return x<rhs.x?x-rhs.x+P:x-rhs.x;}
inline Z operator -()const{ return x?P-x:0;}
inline Z operator *(const Z &rhs)const{ return static_cast<ull>(x)*rhs.x%P;}
inline Z operator +=(const Z &rhs){ return x=x+rhs.x<P?x+rhs.x:x+rhs.x-P, *this;}
inline Z operator -=(const Z &rhs){ return x=x<rhs.x?x-rhs.x+P:x-rhs.x, *this;}
inline Z operator *=(const Z &rhs){ return x=static_cast<ull>(x)*rhs.x%P, *this;}
};

int n, k;
Z ans;
vector<Z> f;

namespace Poly{
Z w[N];// for DFT

inline Z Pow(Z x, int y=P-2){
Z ans=1;
for(; y; y>>=1, x=x*x) if(y&1) ans=ans*x;
return ans;
}
inline void Init(){
for(int i=1; i<N; i<<=1){
w[i]=1;
Z t=Pow(3, (P-1)/i/2);
for(int j=1; j<i; ++j) w[i+j]=w[i+j-1]*t;
}
}
inline int Get(int x){ int n=1; while(n<=x) n<<=1; return n;}
inline void DFT(vector<Z> &f, int n){
static ull F[N];
if((int)f.size()!=n) f.resize(n);
for(int i=0, j=0; i<n; ++i){
F[i]=f[j].x;
for(int k=n>>1; (j^=k)<k; k>>=1);
}
for(int i=1; i<n; i<<=1) for(int j=0; j<n; j+=i<<1){
Z *W=w+i;
ull *F0=F+j, *F1=F+j+i;
for(int k=j; k<j+i; ++k, ++W, ++F0, ++F1){
ull t=(*F1)*(W->x)%P;
(*F1)=*F0+P-t, (*F0)+=t;
}
}
for(int i=0; i<n; ++i) f[i]=F[i]%P;
}
inline void IDFT(vector<Z> &f, int n){
f.resize(n), reverse(f.begin()+1, f.end());
DFT(f, n);
Z I=Pow(n);
for(int i=0; i<n; ++i) f[i]=f[i]*I;
}
inline vector<Z> operator *(const vector<Z> &f, const vector<Z> &g){
if((ull)f.size()*g.size()<=1000){
vector<Z> ans;
ans.resize(f.size()+g.size()-1);
for(unsigned i=0; i<f.size(); ++i) for(unsigned j=0; j<g.size(); ++j)
ans[i+j]+=f[i]*g[j];
return ans;
}
static vector<Z> F, G;
F=f, G=g;
int p=Get(f.size()+g.size()-2);
DFT(F, p), DFT(G, p);
for(int i=0; i<p; ++i) F[i]*=G[i];
IDFT(F, p);
return F.resize(f.size()+g.size()-1), F;
}
vector<Z> PolyInv(const vector<Z> &f, int n=-1){
if(n==-1) n=f.size();
if(n==1) return {Pow(f[0])};
vector<Z> ans=PolyInv(f, (n+1)/2), tmp(&f[0], &f[0]+n);
int p=Get(n*2-2);
DFT(tmp, p), DFT(ans, p);
for(int i=0; i<p; ++i) ans[i]=((Z)2-ans[i]*tmp[i])*ans[i];
IDFT(ans, p);
return ans.resize(n), ans;
}
// a=d*b+r
inline void PolyDiv(const vector<Z> &a, const vector<Z> &b, vector<Z> &d, vector<Z> &r){
if(b.size()>a.size()) return d.clear(), (void)(r=a);

vector<Z> 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=PolyInv(B, n-m+1);
d=A*iB;
d.resize(n-m+1), reverse(d.begin(), d.end());

r=b*d, r.resize(m-1);
for(int i=0; i<m-1; ++i) r[i]=a[i]-r[i];
}
inline void Sqr(vector<Z> &f){
int n=f.size();
if((ull)n*n<=1000){
vector<Z> ans(n*2-1);
for(int i=0; i<n; ++i) for(int j=0; j<n; ++j) ans[i+j]+=f[i]*f[j];
f=ans;
return;
}
int p=Get(n*2-2);
DFT(f, p);
for(int i=0; i<p; ++i) f[i]*=f[i];
IDFT(f, p);
f.resize(n*2-1);
}
inline vector<Z> solve(int n, const vector<Z> &a){
if(n==1) return {0, 1};
vector<Z> ans=solve(n>>1, a), t;
Sqr(ans);
if(n&1) ans.insert(ans.begin(), 0);
PolyDiv(ans, a, t, t);
return t;
}
}
int main() {
Poly::Init();
read(n), read(k), f.resize(k+1), f[k]=1;
for(int i=k, x; i--;) read(x), f[i].x=(P-x%P)%P;
f=Poly::solve(n, f);
for(int i=0, x; i<k; ++i) read(x), ans+=f[i]*(x%P+P);
return printf("%d", ans.x), 0;
}

Shlw loves matrixI

BZOJ 4161: Shlw loves matrixI

可以暴力取模,由于模数原因也可以用 MTT

代码

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
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
#include<cstdio>
#include<algorithm>
#include<cctype>
#include<string.h>
#include<cmath>
#include<vector>

using namespace std;
#define ll long long

const double Pi = acos(-1);
const int N = 1<<13, P = 1000000007;
int n, k, ans;
vector<int> a;

struct cp{
double a, b;
inline void operator +=(const cp &rhs){ a+=rhs.a, b+=rhs.b;}
inline cp operator +(const cp &rhs)const{ return (cp){a+rhs.a, b+rhs.b};}
inline cp operator -(const cp &rhs)const{ return (cp){a-rhs.a, b-rhs.b};}
inline cp operator *(const cp &rhs)const{ return (cp){a*rhs.a-b*rhs.b, a*rhs.b+b*rhs.a};}
inline cp operator *(const double rhs)const{ return (cp){a*rhs, b*rhs};}
inline cp operator ~()const{ return (cp){a, -b};}
} w[N];
inline void Init(){
for(int i=1; i<N; i<<=1){
w[i]=(cp){1, 0};
for(int j=1; j<i; ++j)
w[i+j]=((j&31)==1?(cp){cos(Pi*j/i), sin(Pi*j/i)}:w[i+j-1]*w[i+1]);
}
}
inline int Pow(ll x, int y=P-2){
int ans=1;
for(; y; y>>=1, x=x*x%P) if(y&1) ans=ans*x%P;
return ans;
}
inline int Get(int x){ int n=1; while(n<=x) n<<=1; return n;}
inline void DFT_(cp *f, int n){
for(register int i=0, j=0; i<n; ++i){
if(i>j) swap(f[i], f[j]);
for(register int k=n>>1; (j^=k)<k; k>>=1);
}
for(register int i=1; i<n; i<<=1) for(register int j=0; j<n; j+=i<<1)
for(register int k=j; k<j+i; ++k){
cp t=w[i+k-j]*f[k+i];
f[k+i]=f[k]-t, f[k]+=t;
}
}
inline void DFT(cp *f, int n){
if(n==1) return;
n>>=1;
static cp a[N/2];
for(register int i=0; i<n; ++i) a[i]=(cp){f[i<<1].a, f[i<<1|1].a};
DFT_(a, n);
for(register int i=0; i<n; ++i){
cp q=~a[(n-i)&(n-1)], x=(a[i]+q)*.5, y=(a[i]-q)*(cp){0, -.5}, t=y*w[n+i];
f[i]=x+t, f[n+i]=x-t;
}
}
inline void IDFT(cp *f, int n){
if(n==1) return;
reverse(f+1, f+n), n>>=1;
static cp a[N/2];
for(register int i=0; i<n; ++i)
a[i]=(f[i]+f[i+n])*.5 + (f[i]-f[i+n])*(cp){0, .5}*w[n+i];
DFT_(a, n);
double k=1./n;
for(register int i=0; i<n; ++i) f[i<<1]=(cp){a[i].a*k, 0}, f[i<<1|1]=(cp){a[i].b*k, 0};
}
inline vector<int> operator *(const vector<int> &f, const vector<int> &g){
vector<int> ans(f.size()+g.size()-1);
if((ll)f.size()*g.size()<=1000){
for(unsigned i=0; i<f.size(); ++i) for(unsigned j=0; j<g.size(); ++j)
ans[i+j]=(ans[i+j]+(ll)f[i]*g[j])%P;
return ans;
}
int l=Get(f.size()+g.size()-2);
static cp f0[N], f1[N], g0[N], g1[N], A[N], B[N], C[N];
memset(f0, 0, sizeof(cp)*l), memset(f1, 0, sizeof(cp)*l);
memset(g0, 0, sizeof(cp)*l), memset(g1, 0, sizeof(cp)*l);
for(unsigned i=0; i<f.size(); ++i) f0[i].a=f[i]&32767, f1[i].a=f[i]>>15;
for(unsigned i=0; i<g.size(); ++i) g0[i].a=g[i]&32767, g1[i].a=g[i]>>15;
DFT(f0, l), DFT(f1, l), DFT(g0, l), DFT(g1, l);
for(int i=0; i<l; ++i)
A[i]=f1[i]*g1[i], B[i]=f1[i]*g0[i]+f0[i]*g1[i], C[i]=f0[i]*g0[i];
IDFT(A, l), IDFT(B, l), IDFT(C, l);
for(unsigned i=0; i<ans.size(); ++i)
ans[i]=(((ll)(A[i].a+.5)%P<<30)+((ll)(B[i].a+.5)<<15)+(ll)(C[i].a+.5))%P;
return ans;
}
vector<int> PolyInv(const vector<int> &f, int n=-1){
if(n==-1) n=f.size();
if(n==1){
vector<int> ans;
return ans.push_back(Pow(f[0])), ans;
}
vector<int> ans=PolyInv(f, (n+1)/2), tmp(&f[0], &f[0]+n);
tmp=tmp*ans*ans, tmp.resize(n);
for(unsigned i=0; i<ans.size(); ++i)
tmp[i]=(2*ans[i]-tmp[i])%P, tmp[i]=(tmp[i]<0?tmp[i]+P:tmp[i]);
for(int i=ans.size(); i<n; ++i) tmp[i]=(tmp[i]?P-tmp[i]:0);
return tmp;
}
inline void PolyDiv(const vector<int> &a, const vector<int> &b, vector<int> &r){
if(b.size()>a.size()) return (void)(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=PolyInv(B, n-m+1);
r=A*iB;
r.resize(n-m+1), reverse(r.begin(), r.end());

r=b*r, r.resize(m-1);
for(int i=0; i<m-1; ++i) r[i]=(a[i]-r[i]<0?a[i]-r[i]+P:a[i]-r[i]);
}
inline vector<int> Solve(int n){
if(n==1) return {0, 1};
vector<int> ans=Solve(n>>1), t;
ans=ans*ans;
if(n&1) ans.insert(ans.begin(), 0);
return PolyDiv(ans, a, t), t;
}

int main() {
Init();
scanf("%d%d", &n, &k), a.resize(k+1), a[k]=1;
for(int i=k; i--;) scanf("%d", &a[i]), a[i]=(P-a[i])%P;
a=Solve(n);
for(int i=0, x; i<k; ++i) scanf("%d", &x), ans=(ans+(ll)a[i]*(x+P))%P;
return printf("%d", ans), 0;
}

【NOI2017】泳池

咕咕咕

0%