bzoj4161 Shlw loves matrixI

Description 给定数列 {hn}前k项，其后每一项满足 hn = a1*h(n-1) + a2*h(n-2) + … +
ak*h(n-k) 其中 a1,a2…ak 为给定数列。请计算 h(n)，并将结果对 1000000007 取模输出。

Input 第 1 行包含两个整数 n,k 第 2 行包含 k 个整数 a1,a2…ak 第 3 行包含 k 个整数
h1,h2…hk

Output

一行一个整数 hn mod 1000000007

普通的矩阵乘法是$O(k^3\log n)$的，显然会T。
观察转移矩阵$M$，枚举第一行选什么可以看出它的特征多项式是$f(x)=x^k-\sum _ {i=1} ^ k a_i x^{k-i}$，也就是$M^i=\sum_{j=1}^ka_jM^{i-j}$。这样的话，任何一个$M$的幂都可以表示成$M^0,M^1,\dots,M^{k-1}$的线性组合，只要不断往前滚动就行了。于是矩阵乘法变成了卷积运算，复杂度降到了$O(k^2)$。
得到了$M^n=\sum_{i=0}^{k-1}c_iM^i$以后，答案就是$M^nH=\sum_{i=0}^{k-1}c_iM^iH$的第一行元素。于是，只需要把后面的元素加起来就行了，最后的答案是$\sum_{i=0}^{k-1}c_ih_{k+i-1}$，$h_{k-1}\cdots k_{2k-2}$可以暴力$O(k^2)$算出。
总的复杂度$O(k^2\log n)$，需要一些常数优化。

#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
#define LL long long
const int p=1000000007;
int n,k,f[2010],h[4010];
int rd()
{
    int x=0,f=1;
    char c=getchar();
    while (c<'0'||c>'9')
    {
        if (c=='-') f=-1;
        c=getchar();
    }
    while (c>='0'&&c<='9')
    {
        x=x*10+c-'0';
        c=getchar();
    }
    return x*f;
}
struct mat
{
    int a[4010];
    void init()
    {
        a[0]=1;
        for (int i=1;i<k;i++) a[i]=0;
    }
    mat operator * (const mat &mm) const
    {
        mat ret;
        for (int i=0;i<=2*k-2;i++) ret.a[i]=0;
        for (int i=0;i<k;i++)
            for (int j=0;j<k;j++)
                ret.a[i+j]=(ret.a[i+j]+(LL)a[i]*mm.a[j])%p;
        for (int i=2*k-2;i>=k;i--)
            for (int j=1;j<=k;j++)
                ret.a[i-j]=(ret.a[i-j]+(LL)ret.a[i]*f[j])%p;
        return ret;
    }
}base,res;
int main()
{
    int ans=0;
    n=rd();
    k=rd();
    for (int i=1;i<=k;i++)
    {
        f[i]=rd();
        if (f[i]<0) f[i]+=p;
    }
    for (int i=0;i<k;i++)
    {
        h[i]=rd();
        if (h[i]<0) h[i]+=p;
    }
    res.init();
    base.a[1]=1;
    for (n-=k-1;n;n>>=1,base=base*base)
        if (n&1) res=res*base;
    for (int i=k;i<2*k-1;i++)
        for (int j=1;j<=k;j++)
            h[i]=(h[i]+(LL)h[i-j]*f[j])%p;
    for (int i=0;i<k;i++)
        ans=(ans+(LL)res.a[i]*h[k+i-1])%p;
    printf("%d\n",ans);
}