常系数齐次线性递推

题意：

$h(i)=sum_{j=1}^{k} {h(i-j)*f(j)}$

给出$h(0..k-1)$ 求$h(n)$

题解：

学了无限久。。

首先矩阵快速幂求复杂度是$logn*k^3$的

因为转移矩阵具有一些特性，所以可以优化复杂度

所以这个算法也只能对于线性递推

$$f(alpha)=det(alpha*I-M)$$ 这个东西是它的特征多项式

然后去计算这个东西的行列式（可以用数学归纳法）

得到$$f(alpha)=alpha^k-sum_{i=1}^{k} {M[i]* alpha^{k-i}}$$

而这个东西叫做它的特征值

然后根据Cayley-Hamilton定理矩阵的特征多项式是它的零化多项式

而矩阵的零化多项式满足$M'=0$ 即$f(M)=0$

于是我们可以把当前多项式对$f(M)$取模

然后我们等价于求

$$x^n \% (x^k-sum_{i=0}^{k-1}{f[k-i]*x^i} )$$

取模的话就是

dep(i,2*k-2,k)
      rep(j,1,k)
        z.a[i-j]=(z.a[i-j]+1ll*z.a[i]*f[j])%mo;

因为第一项为1，后面的为负号，所以就直接这么加了

这样一次取模是$k^2$的可以利用fft优化多项式取模优化到$klogk$

这样我们可以得到$M^n=sum_{i=0}^{k-1} {M^i*c[i]}$

关键在于如何快速求这个式子

如果求完整的那么复杂度反而更高了

我们注意到我们只需要$a[n]$这项

所以我们两边同乘$vo$ $$vo*M^n=sum_{i=0}^{k-1} {vo*M^i*c[i]}$$

而$vo*M^k$的第一项为$a[k]$ (网上有代码用了最后一项这样写起来会麻烦一点而且得特判n<k）

所以就是$$a[n]=sum_{i=0}^{k-1} {a[i]*c[i]}$$

这样就可以在$logn*k^2$内解决了

代码:

#include <bits/stdc++.h>
using namespace std;
#define rint register int
#define IL inline
#define rep(i,h,t) for(int i=h;i<=t;i++)
#define dep(i,t,h) for(int i=t;i>=h;i--)
#define ll long long
#define me(x) memset(x,0,sizeof(x))
#define mep(x,y) memcpy(x,y,sizeof(y))
#define mid ((h+t+1)>>1)
namespace IO{
    char ss[1<<24],*A=ss,*B=ss;
    IL char gc()
    {
        return A==B&&(B=(A=ss)+fread(ss,1,1<<24,stdin),A==B)?EOF:*A++;
    }
    template<class T> void read(T &x)
    {
        rint f=1,c; while (c=gc(),c<48||c>57) if (c=='-') f=-1; x=(c^48);
        while (c=gc(),c>47&&c<58) x=(x<<3)+(x<<1)+(c^48); x*=f; 
    }
    char sr[1<<24],z[20]; int Z,C1=-1;
    template<class T>void wer(T x)
    {
        if (x<0) sr[++C1]='-',x=-x;
        while (z[++Z]=x%10+48,x/=10);
        while (sr[++C1]=z[Z],--Z);
    }
    IL void wer1()
    {
        sr[++C1]=' ';
    }
    IL void wer2()
    {
        sr[++C1]='
';
    }
    template<class T>IL void maxa(T &x,T y) {if (x<y) x=y;}
    template<class T>IL void mina(T &x,T y) {if (x>y) x=y;} 
    template<class T>IL T MAX(T x,T y){return x>y?x:y;}
    template<class T>IL T MIN(T x,T y){return x<y?x:y;}
};
using namespace IO;
const int mo=1e9+7;
const int N=4010;
int f[N],h[N],n,k;
struct re{
    int a[N];
    re() {me(a);}
}a;
re js(re x,re y)
{
    re z;
    rep(i,0,k-1)
      rep(j,0,k-1)
        z.a[i+j]=(z.a[i+j]+1ll*x.a[i]*y.a[j])%mo;
    dep(i,2*k-2,k)
      rep(j,1,k)
        z.a[i-j]=(z.a[i-j]+1ll*z.a[i]*f[j])%mo;
    return z;
}
re ksm(int x)
{
    if (x==1) return(a);
    re now=ksm(x/2);
    now=js(now,now);
    if (x%2==1) now=js(now,a);
    return now;
}
int main()
{
    //转移矩阵为f 初始为h
    freopen("1.in","r",stdin);
//    freopen("1.out","w",stdout);
    read(n); read(k);
    rep(i,1,k) 
    {
      read(f[i]);
      if (f[i]<0) f[i]+=mo;
    }
    rep(i,0,k-1)
    {
        read(h[i]);
        if (h[i]<0) h[i]+=mo;
    }
    a.a[1]=1;
    re now=ksm(n);
    ll ans=0;
    rep(i,0,k-1)
      ans=(ans+1ll*now.a[i]*h[i])%mo; 
   /* re now=ksm(n-k+1);
    rep(i,k,2*k-1)
      rep(j,1,k)
        h[i]=(h[i]+1ll*h[i-j]*f[j])%mo;
    ll ans=0;
    rep(i,0,k-1)
      ans=(ans+1ll*now.a[i]*h[k+i-1])%mo; */ 
    cout<<ans<<endl;
    return 0;
}