问题引入
给出数列(g)，满足当(n>m)时
当(n<=m)时,(g_n=c_n)
m比较小，n特别大，快速计算(g_n)
Newbie的解法
暴力递推计算
时间复杂度(O(nm))
Pupil的解法
可以将转移和数列都写成(m imes m)的矩阵的形式，矩阵快速幂即可
时间复杂度(O(m^3log n))
Master的解法
我们需要一些数学知识进行铺垫：
Part 1 矩阵的特征值与特征多项式
我们知道一个矩阵乘一个列向量仍然是一个列向量。
若对于m阶矩阵A，有常数(lambda)，非零列向量(vec v)，满足$$lambdavec v=Avec v$$则称(lambda)为矩阵A的特征值，(vec v)为矩阵的特征向量
上式也可以写作$$(lambda I-A)vec v=0$$其中(I)为单位矩阵
此式有解的充要条件是(|lambda I-A|=0)，即矩阵(lambda I-A)的行列式为0
(|lambda I-A|)可以看做是关于(lambda)的一个m次多项式，记作(f(lambda))，(f(lambda))称作矩阵A的特征多项式，对于矩阵A的任意一个特征值(lambda_0)，都有(f(lambda_0)=0)。
Part 2 Hamilton-Cayley theorem
对于矩阵，也一样的定义多项式运算（把多项式中的x换乘矩阵A），加法就是直接对应相加，常数乘法就按位相乘，乘法是矩阵乘法，0次方是单位矩阵，它的结果仍然是一个矩阵。
显然，矩阵多项式满足交换律，即(f(A)g(A)=g(A)f(A))成立。
简单证明：考虑某两项相乘的结果(A^x imes A^y)，由于前后都是A，矩阵乘法满足结合律，因此指数可以任意分配，换成(A^y imes A^x)也是可以的
哈密顿—凯莱定理：对于矩阵A的特征多项式(f(x))，满足(f(A)=0)
证明网上到处都有，此处就不赘述了。
Part 3 求解转移矩阵的特征多项式
回到原题，我们对于Pupil解法的转移矩阵A，求解它的特征多项式
考虑矩阵(lambda I-A)
它长这样：
根据行列式的定义，将第一行展开
(|lambda I-A|=(lambda-a_1)A_{1,1}+a_2 imes A_{1,2}+cdots+a_m imes A_{1,m})
其中(A{i,j})表示矩阵A的代数余子式，即挖掉第i行和第j列以后剩下的矩阵的行列式。
我们发现所有的余子矩阵都是下三角矩阵，行列式就是对角线乘积。
化简整理，可得(f(lambda)=|lambda I-A|=lambda^m-sumlimits_{i=0}^{m-1}a_{m-i}lambda ^i)
Part 4 计算答案
我们设要求的数列(g)的初始矩阵为(G)，它是一个m行1列的矩阵（列向量），从第m行到第1行分别为(g_{1dots m})（注意顺序是反的）
实际上我们想知道的(g_n)就是矩阵(A^{n-1}G)的第m行第一列的值。
此时的关键就是(A^{n-1})，因为(n-1)非常大，无法直接计算
然而根据前面的铺垫，我们有(f(A)=0)，(A^{n-1})我们可以看做只有一项的一个关于A的多项式
那么根据多项式除法相关知识，可以得到(A^{n-1}=P(A)f(A)+Q(A))，其中(Q(A))的次数是小于(f(A))的次数也就是小于m的，(Q(A))相当于多项式(A^{n-1})对多项式(f(A))取模
可能会有这样的疑问，(f(A)=0)怎么能作除数呢？
其实并不要紧，我们并不需要知道(f(A))的实际值，我们相当于将(A^{n-1})减去了若干个(f(A))，将次数降低了，而结果不变。
实现上来说，由于(f)的系数已知，我们可以先将式子里的矩阵A换成变量(x)，代入，利用多项式取模算出Q的系数，然后再将x换回A，这样得出来的Q的系数是相同的。并且计算(Q(A) imes G)与(A^{n-1} imes G)的结果是一样的。
为了求出(Q(x))的系数，我们可以采用快速幂的做法，初始(Q_0(x)=x^1)，然后不断的自己与自己相乘，乘完对多项式(f(x))取模
这一部分如果暴力取模，时间复杂度为(O(m^2log n))
如果采用NTT优化多项式取模，时间复杂度为(O(mlog mlog n))
这样求出了(Q(A))的系数，不妨设(Q(A)=sumlimits_{i=0}^{m-1}d_iA^i)
要求矩阵(Q(A) imes G)的第m行第一列的值
也就是(sumlimits_{i=0}^{m-1}d_iA^iG)的第m行第一列
然而(A^iG)的第m行第一列的值就是(g_{i+1})
所以$$g_n=sumlimits_{i=0}^{m-1}d_ig_{i+1}=sumlimits_{i=0}^{m-1}d_ic_{i+1}$$
还有一种情况，前m项并没有直接给出，也是通过递推得出的，暴力递推求前m项的复杂度是(O(m^2))的
考虑优化
考虑数列(g)的一般生成函数(G(x))（与矩阵G不同）
转移序列(a)的一般生成函数(A(x))
由于(G(x))是无限长的一个序列，我们可以得到(G(x)=G(x)A(x)+r)
其中(r)是一个常数，相当于第0项
移项，可以得到$$G(x)={rover 1-A(x)}$$
在模(x^{m+1})意义下多项式求逆即可
时间复杂度是(O(mlog m))的
模板题（[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
",ans);
}