bzoj 4161: Shlw loves matrixI 常系数线性齐次递推

题意

给定数列 {hn}前k项，其后每一项满足
hn = a1*h(n-1) + a2*h(n-2) + … + ak*h(n-k)
其中 a1,a2…ak 为给定数列。请计算 h(n)，并将结果对 1000000007 取模输出。
对于 100% 数据，满足 n <= 10^9;k <= 2000; abs(hi)<=10^9; abs(ai)<=10^9

分析

搞了一下午中午搞了出来。
大概就是先把状态转移矩阵$A$画出来，不难发现$A$的特征多项式就是$x^k-\sum\limits_{i=1}^ka_i*x^{k-i}$。
设$f(x)=x^k-\sum\limits_{i=1}^ka_i*x^{k-i}$，根据Cayley-Hamilton定理可以得到$f(A)=0$。
那么$A^n=A^n\bmod f(A)$。
于是我们可以先快速幂求出$x^n\bmod f(x)$的值，设为$\sum\limits_{i=0}^{k-1}c_i*x^i$然后把$A$带到$x$里面，我们要求的是$A^nH=\sum\limits_{i=0}^{k-1}c_i*A^i*H$的最后一个元素。
就等于$\sum\limits_{i=0}^{k-1}c_i*h_{i+k}$。
那么我们可以先$O(k^2)$把$h_{k...2k}$求出来然后直接算即可。

代码

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<algorithm>
using namespace std;

typedef long long LL;

const int N=4005;
const int MOD=1000000007;

int n,k,a[N],c[N],mo[N],tmp[N],ans[N],b[N];

int read()
{
    int x=0,f=1;char ch=getchar();
    while (ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
    while (ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();}
    return x*f;
}

void mul(int *a,int *b,int *c)
{
    for (int i=0;i<=k*2-2;i++) tmp[i]=0;
    for (int i=0;i<k;i++)
        for (int j=0;j<k;j++)
            tmp[i+j]+=(LL)a[i]*b[j]%MOD,tmp[i+j]-=tmp[i+j]>=MOD?MOD:0;
    for (int i=k*2-2;i>=k;i--)
        for (int j=k-1;j>=0;j--)
            tmp[i-k+j]+=MOD-(LL)tmp[i]*mo[j]%MOD,tmp[i-k+j]-=tmp[i-k+j]>=MOD?MOD:0;
    for (int i=0;i<k;i++) c[i]=tmp[i];
}

int main()
{
    n=read();k=read();
    for (int i=1;i<=k;i++) c[i]=read(),c[i]+=c[i]<0?MOD:0;
    for (int i=1;i<=k;i++) a[i]=read(),a[i]+=a[i]<0?MOD:0;
    for (int i=1;i<=k;i++) mo[k-i]=MOD-c[i];
    mo[k]=1;b[1]=1;ans[0]=1;
    for (n-=k-1;n;n>>=1)
    {
        if (n&1) mul(ans,b,ans);
        mul(b,b,b);
    }
    for (int i=k+1;i<=k*2-1;i++)
        for (int j=1;j<=k;j++)
            a[i]+=(LL)a[i-j]*c[j]%MOD,a[i]-=a[i]>=MOD?MOD:0;
    int tot=0;
    for (int i=0;i<=k-1;i++) tot+=(LL)ans[i]*a[k+i]%MOD,tot-=tot>=MOD?MOD:0;
    printf("%d",tot);
    return 0;
}