【BZOJ4161】Shlw loves matrixI

題目描述

給定數列 {hn}前k項，其後每一項知足
hn = a1h(n-1) + a2h(n-2) + ... + ak*h(n-k)
其中 a1,a2...ak 爲給定數列。請計算 h(n)，並將結果對 1000000007 取模輸出。

解法

一個顯然的思路就是矩陣乘法，但這樣的話顯然超時。
實際上，咱們還能夠繼續對這個矩陣乘法進行優化。

首先，因爲這是常係數齊次線性遞推式，簡單來講就是：
\[f_i=\sum_{j=1}^k a_i*f_{i-j}\]
而後，咱們須要引進特徵多項式這個概念。
對於一個矩陣\(A\)，它的特徵多項式是\(f(x)=|Ix-A|\)
把行列式展開以後得，\(f(x)=|Ix-A|=x^K-\sum_{i=1}^K a_i*x^{K-i}\)
由Cayley-hamilton定理，那麼咱們知道\(f(A)=0\)
而後就能知道一個關鍵的式子：
\[A^K=\sum_{i=1}^K a_i*A^{K-i}\]

而後因爲\(A^i\)都能表示成\(A^1,A^2,....,A^K\)的線性組合，
因此如今矩陣乘法直接使用\(O(K^2)\)的卷積便可。

當獲得了\(A^n=\sum_{i=1}^K c_i*A^i\)以後，咱們乘上題目給的向量h。
那麼就會有\(\sum_{i=1}^K c_i*A^i*h_K\)，即\(Ans=\sum_{i=1}^K c_i*h_{K+i}\)。

複雜度就被優化爲\(O(K^2 log n)\)

Code

#include<bits/stdc++.h>
#define ll long long
#define fo(i,x,y) for(int i=x;i<=y;i++)
#define fd(i,x,y) for(int i=x;i>=y;i--)
using namespace std;

const int maxn=4007,mo=1e9+7;

int n,K;
int a[maxn],h[maxn],f[maxn];
int Ans;

void Init(){
    scanf("%d%d",&n,&K);
    fo(i,1,K) scanf("%d",&a[i]);
    memcpy(f,a,sizeof a);
    fo(i,1,K) scanf("%d",&h[i]);
}

#define PLUS(x,y) (x)=((x)+(y))%mo
void mult(int *a,int *b){
    static int c[maxn];
    memset(c,0,sizeof c);
    fo(i,1,K)
        fo(j,1,K)
            PLUS(c[i+j],1ll*a[i]*b[j]);
    fd(i,2*K,K+1)
        fo(j,1,K)
            PLUS(c[i-j],1ll*c[i]*f[j]);
    memcpy(a,c,sizeof c);
}

void qPower(int x){
    bool bz=false;
    static int b[maxn];
    b[1]=1;
    while (x){
        if (x&1){
            if (bz) mult(a,b);
            else{
                bz=true;
                memcpy(a,b,sizeof b);
            }
        }
        mult(b,b);
        x>>=1;
    }
}
void Solve(){
    n++;
    fo(i,K+1,2*K) fo(j,1,K) PLUS(h[i],1ll*h[i-j]*a[j]);
    if (n<=2*K){
        Ans=h[n];
        return;
    }
    qPower(n-K);
    Ans=0;
    fo(i,1,K) PLUS(Ans,1ll*a[i]*h[i+K]);
}

void Print(){
    Ans=(Ans+mo)%mo;
    printf("%d\n",Ans);
}

int main(){
    freopen("1.in","r",stdin);
    freopen("1.out","w",stdout);
    Init();
    Solve();
    Print();
    return 0;
}