问题引入
给出数列
g
g
g,满足当
n
>
m
n>m
n>m时
g
n
=
∑
i
=
1
m
g
n
−
i
×
a
i
g_n=\sum\limits_{i=1}^{m}g_{n-i}\times a_i
gn=i=1∑mgn−i×ai
当
n
<
=
m
n<=m
n<=m时,
g
n
=
c
n
g_n=c_n
gn=cn
m比较小,n特别大,快速计算 g n g_n gn
Newbie的解法
暴力递推计算
时间复杂度 O ( n m ) O(nm) O(nm)
Pupil的解法
可以将转移和数列都写成 m × m m\times m m×m的矩阵的形式,矩阵快速幂即可
时间复杂度 O ( m 3 log n ) O(m^3\log n) O(m3logn)
Master的解法
我们需要一些数学知识进行铺垫:
Part 1 矩阵的特征值与特征多项式
我们知道一个矩阵乘一个列向量仍然是一个列向量。
若对于m阶矩阵A,有常数 λ \lambda λ,非零列向量 v ⃗ \vec v v,满足 λ v ⃗ = A v ⃗ \lambda\vec v=A\vec v λv=Av则称 λ \lambda λ为矩阵A的特征值, v ⃗ \vec v v为矩阵的特征向量
上式也可以写作
(
λ
I
−
A
)
v
⃗
=
0
(\lambda I-A)\vec v=0
(λI−A)v=0其中
I
I
I为单位矩阵
此式有解的充要条件是
∣
λ
I
−
A
∣
=
0
|\lambda I-A|=0
∣λI−A∣=0,即矩阵
λ
I
−
A
\lambda I-A
λI−A的行列式为0
∣ λ I − A ∣ |\lambda I-A| ∣λI−A∣可以看做是关于 λ \lambda λ的一个m次多项式,记作 f ( λ ) f(\lambda) f(λ), f ( λ ) f(\lambda) f(λ)称作矩阵A的特征多项式,对于矩阵A的任意一个特征值 λ 0 \lambda_0 λ0,都有 f ( λ 0 ) = 0 f(\lambda_0)=0 f(λ0)=0。
Part 2 Hamilton-Cayley theorem
对于矩阵,也一样的定义多项式运算(把多项式中的x换乘矩阵A),加法就是直接对应相加,常数乘法就按位相乘,乘法是矩阵乘法,0次方是单位矩阵,它的结果仍然是一个矩阵。
显然,矩阵多项式满足交换律,即
f
(
A
)
g
(
A
)
=
g
(
A
)
f
(
A
)
f(A)g(A)=g(A)f(A)
f(A)g(A)=g(A)f(A)成立。
简单证明:考虑某两项相乘的结果
A
x
×
A
y
A^x\times A^y
Ax×Ay,由于前后都是A,矩阵乘法满足结合律,因此指数可以任意分配,换成
A
y
×
A
x
A^y \times A^x
Ay×Ax也是可以的
哈密顿—凯莱定理:对于矩阵A的特征多项式 f ( x ) f(x) f(x),满足 f ( A ) = 0 f(A)=0 f(A)=0
证明网上到处都有,此处就不赘述了。
Part 3 求解转移矩阵的特征多项式
回到原题,我们对于Pupil解法的转移矩阵A,求解它的特征多项式
考虑矩阵
λ
I
−
A
\lambda I-A
λI−A
它长这样:
(1)
λ
I
−
A
=
(
λ
−
a
1
−
a
2
⋯
−
a
m
−
1
−
a
m
−
1
λ
⋯
0
0
0
−
1
⋯
0
0
⋮
⋮
⋱
⋮
⋮
0
0
⋯
−
1
λ
)
\lambda I-A= \left( { \begin{matrix} \lambda-a_1 & -a_2 & \cdots &-a_{m-1} & -a_m \\ -1 & \lambda & \cdots & 0 &0 \\ 0 & -1 &\cdots & 0 & 0\\ \vdots & \vdots & \ddots & \vdots &\vdots \\ 0 & 0 & \cdots & -1 & \lambda \end{matrix} \tag{1} } \right)
λI−A=⎝⎜⎜⎜⎜⎜⎛λ−a1−10⋮0−a2λ−1⋮0⋯⋯⋯⋱⋯−am−100⋮−1−am00⋮λ⎠⎟⎟⎟⎟⎟⎞(1)
根据行列式的定义,将第一行展开
∣
λ
I
−
A
∣
=
(
λ
−
a
1
)
A
1
,
1
+
a
2
×
A
1
,
2
+
⋯
+
a
m
×
A
1
,
m
|\lambda I-A|=(\lambda-a_1)A_{1,1}+a_2\times A_{1,2}+\cdots+a_m\times A_{1,m}
∣λI−A∣=(λ−a1)A1,1+a2×A1,2+⋯+am×A1,m
其中
A
i
,
j
A{i,j}
Ai,j表示矩阵A的代数余子式,即挖掉第i行和第j列以后剩下的矩阵的行列式。
我们发现所有的余子矩阵都是下三角矩阵,行列式就是对角线乘积。
化简整理,可得
f
(
λ
)
=
∣
λ
I
−
A
∣
=
λ
m
−
∑
i
=
0
m
−
1
a
m
−
i
λ
i
f(\lambda)=|\lambda I-A|=\lambda^m-\sum\limits_{i=0}^{m-1}a_{m-i}\lambda ^i
f(λ)=∣λI−A∣=λm−i=0∑m−1am−iλi
负号都被行列式里面逆序对个数的负号消掉了。
Part 4 计算答案
我们设要求的数列
g
g
g的初始矩阵为
G
G
G,它是一个m行1列的矩阵(列向量),从第m行到第1行分别为
g
1
…
m
g_{1\dots m}
g1…m(注意顺序是反的)
实际上我们想知道的
g
n
g_n
gn就是矩阵
A
n
−
1
G
A^{n-1}G
An−1G的第m行第一列的值。
此时的关键就是 A n − 1 A^{n-1} An−1,因为 n − 1 n-1 n−1非常大,无法直接计算
然而根据前面的铺垫,我们有 f ( A ) = 0 f(A)=0 f(A)=0, A n − 1 A^{n-1} An−1我们可以看做只有一项的一个关于A的多项式
那么根据多项式除法相关知识,可以得到 A n − 1 = P ( A ) f ( A ) + Q ( A ) A^{n-1}=P(A)f(A)+Q(A) An−1=P(A)f(A)+Q(A),其中 Q ( A ) Q(A) Q(A)的次数是小于 f ( A ) f(A) f(A)的次数也就是小于m的, Q ( A ) Q(A) Q(A)相当于多项式 A n − 1 A^{n-1} An−1对多项式 f ( A ) f(A) f(A)取模
可能会有这样的疑问,
f
(
A
)
=
0
f(A)=0
f(A)=0怎么能作除数呢?
其实并不要紧,我们并不需要知道
f
(
A
)
f(A)
f(A)的实际值,我们相当于将
A
n
−
1
A^{n-1}
An−1减去了若干个
f
(
A
)
f(A)
f(A),将次数降低了,而结果不变。
实现上来说,由于 f f f的系数已知,我们可以先将式子里的矩阵A换成变量 x x x,代入,利用多项式取模算出Q的系数,然后再将x换回A,这样得出来的Q的系数是相同的。并且计算 Q ( A ) × G Q(A)\times G Q(A)×G与 A n − 1 × G A^{n-1}\times G An−1×G的结果是一样的。
为了求出
Q
(
x
)
Q(x)
Q(x)的系数,我们可以采用快速幂的做法,初始
Q
0
(
x
)
=
x
1
Q_0(x)=x^1
Q0(x)=x1,然后不断的自己与自己相乘,乘完对多项式
f
(
x
)
f(x)
f(x)取模
这一部分如果暴力取模,时间复杂度为
O
(
m
2
log
n
)
O(m^2\log n)
O(m2logn)
如果采用NTT优化多项式取模,时间复杂度为
O
(
m
log
m
log
n
)
O(m\log m\log n)
O(mlogmlogn)
这样求出了
Q
(
A
)
Q(A)
Q(A)的系数,不妨设
Q
(
A
)
=
∑
i
=
0
m
−
1
d
i
A
i
Q(A)=\sum\limits_{i=0}^{m-1}d_iA^i
Q(A)=i=0∑m−1diAi
要求矩阵
Q
(
A
)
×
G
Q(A)\times G
Q(A)×G的第m行第一列的值
也就是
∑
i
=
0
m
−
1
d
i
A
i
G
\sum\limits_{i=0}^{m-1}d_iA^iG
i=0∑m−1diAiG的第m行第一列
然而
A
i
G
A^iG
AiG的第m行第一列的值就是
g
i
+
1
g_{i+1}
gi+1
所以 g n = ∑ i = 0 m − 1 d i g i + 1 = ∑ i = 0 m − 1 d i c i + 1 g_n=\sum\limits_{i=0}^{m-1}d_ig_{i+1}=\sum\limits_{i=0}^{m-1}d_ic_{i+1} gn=i=0∑m−1digi+1=i=0∑m−1dici+1
还有一种情况,前m项并没有直接给出,也是通过递推得出的,暴力递推求前m项的复杂度是
O
(
m
2
)
O(m^2)
O(m2)的
考虑优化
考虑数列
g
g
g的一般生成函数
G
(
x
)
G(x)
G(x)(与矩阵G不同)
转移序列
a
a
a的一般生成函数
A
(
x
)
A(x)
A(x)
由于
G
(
x
)
G(x)
G(x)是无限长的一个序列,我们可以得到
G
(
x
)
=
G
(
x
)
A
(
x
)
+
r
G(x)=G(x)A(x)+r
G(x)=G(x)A(x)+r
其中
r
r
r是一个常数,相当于第0项
移项,可以得到
G
(
x
)
=
r
1
−
A
(
x
)
G(x)={r\over 1-A(x)}
G(x)=1−A(x)r
在模
x
m
+
1
x^{m+1}
xm+1意义下多项式求逆即可
时间复杂度是
O
(
m
log
m
)
O(m\log m)
O(mlogm)的
模板题([BZOJ4161] Shlw loves matrixI)
Code
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <iostream>
#include <algorithm>
#define fo(i,a,b) for(int i=a;i<=b;++i)
#define fod(i,a,b) for(int i=a;i>=b;--i)
#define N 4005
#define mo 1000000007
#define LL long long
using namespace std;
LL f[N],g[N],h[N],s1[N],a[N],u1[N];
int n,m;
void mul(LL *x,LL *y,LL *z)
{
fo(i,0,2*m-2) u1[i]=0;
fo(i,0,m-1) fo(j,0,m-1) u1[i+j]=(u1[i+j]+x[i]*y[j])%mo;
fod(i,2*m-2,m)
{
fo(j,0,m) u1[i-m+j]=(u1[i-m+j]-f[j]*u1[i])%mo;
}
fo(i,0,m-1) z[i]=u1[i];
}
int main()
{
cin>>n>>m;
fo(i,1,m) scanf("%lld",&a[i]),f[m-i]=-a[i];
f[m]=1;
g[1]=1;
s1[0]=1;
for(int t=n;t;t>>=1)
{
if(t&1) mul(s1,g,s1);
mul(g,g,g);
}
fo(i,0,m-1) scanf("%lld",&h[i]);
LL ans=0;
fo(i,0,m-1) ans=(ans+s1[i]*h[i]%mo+mo)%mo;
printf("%lld\n",ans);
}