常系数齐次线性递推-矩阵快速幂的优化-学习笔记
i207M
2019-01-26 20:48:57
常系数：系数都是常数
齐次：没有常数项
线性：没有高次项
通常的解法：矩阵快速幂，$O(k^3\log k)$
我们有优化的方法。
考虑通常的解法：
根据矩阵乘法的结合律，易得到：
$$\left( \begin{bmatrix} f_1 &f_2 &f_3 &f_4 & \cdots &f_{k-2} &f_{k-1} \\ 1 &0 &0 &0 & \cdots &0 &0 \\ 0 &1 &0 &0 & \cdots &0 &0\\ \cdots & \cdots& \cdots & \cdots & \cdots & \cdots & \cdots\\ 0 &0 &0 &0 & \cdots &1 &0 \end{bmatrix} \right) ^n \times \begin{bmatrix} a_{k-1} \\ a_{k-2} \\ \cdots \\ a_{1} \\ a_{0}\end{bmatrix} =\begin{bmatrix} a_{n+k-1} \\ a_{n+k-2} \\ \cdots \\ a_{n+1} \\ a_{n}\end{bmatrix}$$
所以我们只需要算出 $M^n \times A$，然后取最后一个数就可以了。
我们除了快速幂还有什么方法吗？
我们能不能有一个取巧的办法，使得我们只需要算低次的矩阵乘法？
我们可以将 $M^n$分解为 $M^n=f(M) \times g(M) +R(M)$，且$f(M)=O$ ，所以 $M^n=R(M)$ ， $R(M)$ 是一个次数不超过 k-1 的多项式。
也就是说，我们只要求出来 $M^n \mod f(M)$就万事大吉了。
我们可以多项式快速幂+多项式取模来做，但是，怎么求出$f(M)$呢？
------
我们先准备一些知识：
则有
### Cayley-Hamilton定理
$f(A)=O$（$f$指特征多项式）
证明的话很简单啊，对一个0矩阵求行列式当然是0。
好的，我们会构造了，但是如何优雅地求出它每一项的系数？
注意我们的矩阵是一个上海森堡矩阵。（但这并没有什么用）
我们把矩阵的第一行拉普拉斯展开一下
$f(\lambda)=(\lambda-a_1)A_{1,1}+(-a_2)A_{1,2}+(-a_3)A_{1.3}...+(-a_n)A_{1,n}$
其中$A_{i,j}$表示A的代数**余子式**
余子式的简单理解就是把$(i,j)$所在的行列去掉，剩下的$n-1$阶矩阵的行列式。
显然所有余子式都是**下三角矩阵**，行列式很容易求得
那么我们得到了$f(\lambda)=\lambda^n-\sum_{i=1}^{n}a_i\lambda^{n-i}$
------
于是就当我们算出了$M^n \mod f(M)$，然后怎么计算呢？暴力算出矩阵？复杂度又爆炸了。考虑我们最后要求什么？
我们把多项式的每一项单独考虑，对于$k$次项。我们有一个初值的列向量，还有一个$A^k$乘在一起，还有一个$c_k$作为系数；我们只要求出计算后的第一个位置的值即可。注意前两者的乘积，就是$a_k$！所以我们计算一下$\sum a_kc_k$即可！
~~打算WC回去的火车上写~~
```cpp
void _Inv(int n,int a[],int b[])
{
if(n==1)
{
b[0]=qpow(a[0],md-2);
return;
}
_Inv(n>>1,a,b);
static int c[N];
for(ri i=0; i<n; ++i) c[i]=a[i];
int m=pre(n<<1);
dft(b,m,1),dft(c,m,1);
for(ri i=0; i<m; ++i) b[i]=(2-(LL)b[i]*c[i]%md+md)*b[i]%md;
dft(b,m,-1);
for(ri i=n; i<m; ++i) b[i]=0;
for(ri i=0; i<m; ++i) c[i]=0;
}
void Inv(int n,int a[],int b[])
{
int m=1;
for(; m<n; m<<=1);
_Inv(m,a,b);
for(ri i=n; i<m; ++i) b[i]=0;
}
void Div(int n,int a[],int m,int b[],int d[])
{
static int tmpa[N],tmpb[N],tmpc[N];
int len=n-m+1;
for(ri i=0; i<len; ++i) tmpa[i]=a[n-i-1];
for(ri i=0; i<m; ++i) tmpb[i]=b[m-i-1];
Inv(len,tmpb,tmpc);
int t=pre(len<<1);
dft(tmpa,t,1),dft(tmpc,t,1);
for(ri i=0; i<t; ++i) tmpc[i]=(LL)tmpa[i]*tmpc[i]%md;
dft(tmpc,t,-1);
reverse(tmpc,tmpc+len);
for(ri i=0; i<t; ++i) tmpa[i]=0;
for(ri i=0; i<m; ++i) tmpb[i]=b[i];
for(ri i=0; i<len; ++i) d[i]=tmpc[i];
for(ri i=len; i<t; ++i) d[i]=0;
for(ri i=0; i<t; ++i) tmpc[i]=0;
t=pre(n<<1);
dft(tmpb,t,1),dft(d,t,1);
for(ri i=0; i<t; ++i) d[i]=(LL)tmpb[i]*d[i]%md;
dft(d,t,-1);
for(ri i=0; i<m-1; ++i) d[i]=(a[i]-d[i]+md)%md;
for(ri i=m-1; i<t; ++i) d[i]=0;
for(ri i=0; i<t; ++i) tmpb[i]=0;
}
void mul(int n,int a[],int b[],int Md[])
{
static int c[N];
for(ri i=0; i<n; ++i) c[i]=b[i];
int m=pre(n<<1);
dft(a,m,1),dft(c,m,1);
for(ri i=0; i<m; ++i) a[i]=(LL)a[i]*c[i]%md;
dft(a,m,-1);
for(ri i=0; i<m; ++i) c[i]=0;
Div(m,a,n+1,Md,c);
for(ri i=0; i<n; ++i) a[i]=c[i];
for(ri i=n; i<m; ++i) a[i]=0;
}
void Qpow(int a[],int b,int res[],int m,int Md[])
{
res[0]=1;
for(; b; b>>=1,mul(m,a,a,Md)) if(b&1) mul(m,res,a,Md);
}
int n,k;
int a[N],f[N],c[N],d[N];
signed main()
{
#ifdef M207
freopen("in.in","r",stdin);
// freopen("out.out","w",stdout);
#endif
in(n,k);
a[k]=1;
for(ri i=1; i<=k; ++i) in(a[k-i]);
for(ri i=0; i<k; ++i) a[i]=(md-a[i]%md)%md;
for(ri i=0; i<k; ++i) in(f[i]),f[i]=(f[i]%md+md)%md;
c[1]=1;
Qpow(c,n,d,k,a);
int ans=0;
for(ri i=0; i<k; ++i) ans=(ans+(LL)d[i]*f[i])%md;
out(ans);
return 0;
}
```