题解 P4723 【【模板】线性递推】
shadowice1984
2018-07-07 21:59:47
十分好奇同一个复杂度的代码为啥我跑了780ms有人跑了7800ms
这里介绍一个线性代数黑科技——常系数齐次线性递推式的优化
_________________
### 前置知识：多项式求逆
蛤你不会多项式求逆？出门左转模板区多项式求逆
### 前置知识：多项式除法，多项式取模
蛤你不会多项式取模?出门左转模板区多项式除法
## 本题题解
# 常系数齐次线性递推式第n项的快速计算
### 翻译
首先这个标题讲的不是人话，我们来尝试翻译一下
我们要快速的求一个递推式的第n项
而这个递推式满足以下几个条件
1.它是常系数的，换句话说递推系数和下标n无关
2.它是线性的，换句话说递推式中每一项的次数都是1，没有乱七八糟的2次3次项
3.它是齐次的，换句话说递推式的常数项等于0
(事实上带常数项的东西我们也能做，不过这里暂时不提及了)
____________________
好了如果看不懂文字说明也没关系，我们要算数，算的是$f_{n}$
其中满足当$n \geq k$时
## $f_{n}=\sum_{i=1}^{k}a_{i}f_{n-i}$
如果$n \leq k$的话我们会给你输入这个递推序列的初值
### 暴力
我会矩阵快速幂！(什么你不会矩阵快速幂？，同是黄牌题还是要写一写的)
那么我们发现直接计算的复杂度是$O(k^3logn)$的，这顶多算是将复杂度向n倾斜，根本算不上什么优化
但是我们发现矩阵乘法是有启发性意义的，换句话说我们只需要以一种奇技淫巧计算出转移矩阵的n次幂即可
但是让我们仔细想想，我们真的需要转移矩阵的n次幂吗?
_____________
#### 魔法的开始
其实你并不关心转移矩阵的n次幂，你只关心初值向量$St$（就是那k个初始值构成的数组）乘上这个转移矩阵的n次幂所得到的向量
进一步的讲我们甚至不关心这个向量长什么样，我们只关心这个向量的第一项是什么
所以说矩阵乘法之所以慢是因为你知道的太多了……
(以下为了方便我们假设转移矩阵为A)
那么假设我们构造了一个奥妙重重的序列c使得
## $A^{n}=\sum_{i=0}^{k-1}c_{i}A^{i}$
你可能会说，这有个皮皮用，$O(k^4)$还不如我矩阵快速幂快
但是我们冷静一下，又没让我们求$A^{n}$,我们真正感兴趣的是$St×A^{n}$
所以我们让等式两边同时乘上一个向量$St$
## $St×A^{n}=St×\sum_{i=0}^{k-1}c_{i}A^{i}$
此时矩阵乘法对于矩阵加法是具有分配律的，所以我们可以将向量$St$分配进去
## $St×A^{n}=\sum_{i=0}^{k-1}c_{i}St×A^{i}$
此时还是没个什么用，但是让我们继续冷静一下，我们真正关心的是这个向量的第0项，所以我们发现刚才的等式在取每个向量的第0项的意义下也是成立的
## $Ans=(St×A^{n})_{0}=\sum_{i=0}^{k-1}c_{i}(St×A^{i})_{0}$
等会，$(St×A^{i})_{0}$是什么啊，转移了i次的向量的第0项？每转移一次这个向量的第0项会变成上一个向量的第1项，第1项变成第2项，……以此类推
那$(St×A^{i})_{0}$不就是$St_{i}$吗？
然后请擦亮你的双眼
## $Ans=\sum_{i=0}^{k-1}c_{i}St_{i}$
换句话说只要我们可以求出$c$我们就可以在$O(k)$时间内计算出$Ans$
__________________
#### 多项式科技-构造序列c
前面的推倒过程中我们发现了c序列十分的神奇，有了c序列我们就可以快速计算答案
但是问题来了我们怎么计算c序列呢？
回想一遍c的定义，是满足这个条件的数组
## $A^{n}=\sum_{i=0}^{k-1}c_{i}A^{i}$
那么我们发现两边的次数不等，换句话说这样的c很有可能不存在……
所以我们先将$A^{n}$写成这样的形式
## $A^{n}=Q(A)G(A)+R(A)$
其中$Q,G,R$是一些以矩阵为参数的多项式，然后满足$Q(A)G(A)$的次数为n
然后此时我们发现$R(A)$的次数小于$G(A)$的次数，那么如果我们钦定$G(A)$的次数为k的话，那么我们就可以将$A^{n}$写成这样的形式
## $A^{n}=Q(A)G(A)+\sum_{i=0}^{k-1}c_{i}A^{i}$
此时我们最好祈祷这个多项式$G$十分神奇，它满足$G(A)=0$，假设我们足够幸运的话，$G(A)$的系数会是这样的一个数列，满足
## $\sum_{i=0}^{k}g_{i}A^{i}=0$
那么我们就可以奇迹般的消掉$Q(A)G(A)$这一项，此时
## $A^{n}=\sum_{i=0}^{k-1}c_{i}A^{i}=R(A)$
换句话说c数组就是多项式$R(A)$的系数,求多项式$R(A)$就行了
刚才的关于$R(A)$的方程，我好像在哪里见过……，对了，是多项式取模模板题的等式
换句话说
## $R(A)=A^{n}modG(A)$
好了问题来了怎么求$A^{n}modG(A)$呢？
直接快速幂就行了，只是把平常中快速幂的取模换成了多项式取模就行了
如果不会多项式取模的话可以去隔壁模板区学习
_________________
### 特征多项式
根据刚才的推倒过程我们知道了只要构造一个数列g，使得
## $\sum_{i=0}^{k}g_{i}A^{i}=0$
我们就可以在$O(klogklogn)$时间内构造一个数组c，使得
## $A^{n}=\sum_{i=0}^{k-1}c_{i}A^{i}$
然后就可以在$O(k)$的时间内计算出答案
## $Ans=\sum_{i=0}^{k-1}c_{i}St_{i}$
但是问题是如何构造一个序列g呢，到刚才为止我们都没有用到关于常系数齐次线性递推式的任何知识，换句话说刚才的推倒对于任意矩阵都成立
如果对于任意一个矩阵我们都可以平凡的构造出一个序列g的话那矩乘基本就废了
所以g的构造肯定没那么简单
但是事实上根据一些线性代数的黑科技我们还是可以构造出g的，下面是推倒过程~
## 如果你不愿意进行数学证明，那么可以直接记住这个结论
假设递推系数是$a_{1}...a_{k}$
那么$g_{k-i}=a_{i}$,$g_{k}=1$，记住这个结论之后就可以直接去写代码了
____________________
#### 前置知识:特征值与特征向量
如果等式
## $ (\lambda I-A)v=0$
成立，那么$\lambda$称为A的特征值，v称为A的特征向量
以下是几个我也不会证的结论
1.如果大小为n×n的矩阵A满秩(行列式不为0)，那么A有n组**线性无关的**特征向量
2.如果$Det( \lambda I-A)=0$那么存在向量v使得等式成立，否则不存在这样的向量v
__________________
#### 前置知识:Cayley-Hamilton定理
这里给出一个C-H定理十分不标准的表述，下面的等式在方阵A满秩时恒成立
## $\prod_{k}(\lambda_{k} I-A)=0$
其中$\lambda_{k}$是A的第k个特征值
证明的话我们采取这样的思路
我们不直接证明这个矩阵为0矩阵，而是证明任意向量乘上这个矩阵为0
然后我们要证明任意向量乘上这个东西是0矩阵还是有点难度
但是由于A有n组线性无关的特征向量，因此我们可以将这个任意向量表示为这n个向量的若干倍的和的形式
如果我们可以证明这个矩阵乘上任意特征向量为0矩阵的话，自然任意向量乘这个矩阵就是0矩阵了
问题来了怎么证明它乘任意特征向量为0矩阵呢？
先说明一个小结论
### $(aI-A)(bI-A)=abI^2-bAI-aAI-A^2=(bI-A)(aI-A)$
然后我们就可以证明C-H定理
## $v_{i}\prod_{k}(\lambda_{k} I-A)=v_{i}(\lambda_{i}I-A)\prod_{k \neq i}(\lambda_{k} I-A)$
## $=0×\prod_{k \neq i}(\lambda_{k} I-A)=0$
具体来说就是因为像上面那种括号相乘的矩阵是可交换的，所以我们可以把这个特征向量对应的矩阵从$\Pi$里交换出来，那么这两个东西相乘自然等于0了
_______________________
换句话说多项式
## $\prod_{k}(\lambda_{k} I-A)$
就是我们想要找的$G(A)$
问题来了怎么算这个多项式啊……
这里有个奇怪的结论是上面这个多项式和下面这个多项式的系数是一样的
## $f(\lambda)=Det(\lambda I-A)$
喂这两个函数主元不一样甚至传的参数都不一样你让我怎么相信啊……
但是仔细思考一下上下两个多项式，(这里不强调函数是因为主元不一致，你可以理解成一个二元函数)是否有相同的0点呢?
答案是肯定的，$\lambda_{1},\lambda_{2}...\lambda_{k}$都是这两个多项式的0点，而且这两个多项式中$\lambda$的次数都是一样的(为什么一样会构造证明)，所以这两个多项式的系数自然是一致的
_________________________
那么问题转化为了求函数
## $f(\lambda)=Det(\lambda I-A)$
的系数
我们发现对于一般的矩阵来讲这个多项式的系数会非常的不好求，
但是我们可以将这个常系数齐次线性递推式的矩阵画出来，手玩一下行列式会发现
## $f(\lambda)=(-1)^n(\lambda^n-\sum_{i-1}^{n}a_{i}\lambda^{n-i})$
由于将函数取反之后函数的0点依旧不变，因此实际写代码的时候可以省略那个$(-1)^n$项
____________________________
## 实现
没啥好说的，不要想着卡常数，取模和求逆的时候不要吝啬你的ntt
一个优化是因为每次都是模的同一个多项式，因此可以预处理出$G(A)$和$G_{reverse}^{-1}(A)$的ntt形式，然后可以只求逆一次
然后就是数组无脑清空吧……
上代码(最慢的点796ms)
```C
// luogu-judger-enable-o2
#include<cstdio>
#include<algorithm>
using namespace std;const int N=65536+10;typedef long long ll;const ll mod=998244353;
int n;int k;int rv[20][N];ll rt[20][20];int Len;ll tr1[N];ll tr2[N];long long st[N];long long xs[N];
ll sg[N];ll a[N];ll res[N];ll irg[N];ll q[N];ll rf[N];int DL=-1;ll ans=0;ll ret[N];
inline ll po(ll a,ll p){ll r=1;for(;p;p>>=1,a=a*a%mod)if(p&1)r=r*a%mod;return r;}
inline void ntt(ll* a,int o,int len,int d)//ntt
{
for(int i=0;i<len;i++)if(i<rv[d][i])swap(a[i],a[rv[d][i]]);
for(int k=1,j=1;k<len;k<<=1,j++)
for(int s=0;s<len;s+=(k<<1))
for(int i=s,w=1;i<s+k;i++,w=w*rt[o][j]%mod)
{ll a0=a[i];ll a1=a[i+k]*w%mod;a[i]=(a0+a1)%mod,a[i+k]=(a0+mod-a1)%mod;}
if(o==1){ll inv=po(len,mod-2);for(int i=0;i<len;i++)(a[i]*=inv)%=mod;}
}
inline void poly_inv(ll* a,ll* b,int len)//求逆
{
b[0]=po(a[0],mod-2);
for(int k=1,j=0;k<=len;k<<=1,j++)
{
for(int i=0;i<k;i++)tr1[i]=a[i];for(int i=0;i<k;i++)tr2[i]=b[i];
ntt(tr1,0,k<<1,j);ntt(tr2,0,k<<1,j);
for(int i=0;i<(k<<1);i++)b[i]=tr2[i]*(2+mod-tr1[i]*tr2[i]%mod)%mod;
ntt(b,1,k<<1,j);for(int i=k;i<(k<<1);i++)b[i]=0;
}
}
inline void poly_mod(ll* a)//取模
{
int mi=(k<<1);while(a[--mi]==0);if(mi<k)return;
for(int i=0;i<(Len<<1);i++)rf[i]=0;for(int i=0;i<=mi;i++)rf[i]=a[i];
reverse(rf,rf+mi+1);for(int i=mi-k+1;i<=mi;i++)rf[i]=0;ntt(rf,0,Len<<1,DL+1);
for(int i=0;i<(Len<<1);i++)q[i]=(rf[i]*irg[i])%mod;ntt(q,1,(Len<<1),DL+1);
for(int i=mi-k+1;i<=(Len<<1);i++)q[i]=0;reverse(q,q+mi-k+1);ntt(q,0,(Len<<1),DL+1);
for(int i=0;i<(Len<<1);i++)(q[i]*=sg[i])%=mod;ntt(q,1,(Len<<1),DL+1);
for(int i=0;i<k;i++)(a[i]+=mod-q[i])%=mod;for(int i=k;i<=mi;i++)a[i]=0;
}
int main()
{
for(int i=0;i<=15;i++)
for(int j=0;j<(1<<(i+1));j++)rv[i][j]=(rv[i][j>>1]>>1)|((j&1)<<i);
for(int t=2,j=1;j<=18;t<<=1,j++)rt[0][j]=po(3,(mod-1)/t);
for(int t=2,j=1;j<=18;t<<=1,j++)rt[1][j]=po(332748118,(mod-1)/t);
scanf("%d%d",&n,&k);
for(Len=1;Len<=k;Len<<=1,DL++); //预处理
for(int i=1;i<=k;i++){scanf("%lld",&xs[i]);xs[i]=xs[i]<0?xs[i]+mod:xs[i];}
for(int i=0;i<k;i++){scanf("%lld",&st[i]);st[i]=st[i]<0?st[i]+mod:st[i];}
for(int i=1;i<=k;i++)sg[k-i]=mod-xs[i];sg[k]=1;for(int i=0;i<=k;i++)ret[i]=sg[i];
for(int i=0;i<=k;i++)rf[i]=sg[i];reverse(rf,rf+k+1);poly_inv(rf,irg,Len);
for(int i=0;i<=k;i++)rf[i]=0;ntt(sg,0,Len<<1,DL+1);ntt(irg,0,Len<<1,DL+1);a[1]=1;res[0]=1;
while(n)//快速幂
{
if(n&1)
{
ntt(res,0,Len<<1,DL+1);ntt(a,0,Len<<1,DL+1);
for(int i=0;i<(Len<<1);i++)(res[i]*=a[i])%=mod;
ntt(res,1,Len<<1,DL+1);ntt(a,1,Len<<1,DL+1);poly_mod(res);
}ntt(a,0,Len<<1,DL+1);for(int i=0;i<(Len<<1);i++)(a[i]*=a[i])%=mod;
ntt(a,1,Len<<1,DL+1);poly_mod(a);n>>=1;
}for(int i=0;i<k;i++)(ans+=res[i]*st[i])%=mod;printf("%lld",ans);return 0;
}
```