[模板] 常系数线性递推

常系数线性递推

给定向量 \(A_0 = (a_1, a_2, \dotsc, a_k)\), 和向量 \(H = (h_1, h_2, \dotsc, h_k)\), 同时

\[
a_n = \sum_{i=1}^k a_{n-i} h_i
\]

求 \(a_n\).

算法

我们只需求出 \(A_n = (a_n, a_{n+1}, \dotsc, a_{n+k-1})\) 即可.

设 \(f(\lambda)\) 表示转移方程的特征多项式, 有

\[
f(\lambda) = \lambda^k - \sum_{i=0}^{k-1} h_{k-i} \lambda^i
\]

设 \(g(\lambda) \equiv \lambda^{n-1} \pmod{f(\lambda)}\), 那么有

\[
a_n = \sum_{i=0}^{k-1} g_i a_{i+1}
\]

求 \(g(\lambda)\) 如果直接取模则复杂度为 \(O(k^2\log n)\), 多项式取模则为 \(O(k\log k \log n)\).

代码

bzoj4161: Shlw loves matrixI

#include<cstdio>
#include<iostream>
#include<cmath>
#include<cstring>
#include<algorithm>
#include<queue>
#include<map>
#include<set>
using namespace std;
#define rep(i,l,r) for(register int i=(l);i<=(r);++i)
#define repdo(i,l,r) for(register int i=(l);i>=(r);--i)
#define il inline
typedef long long ll;
typedef double db;

//---------------------------------------
const int nsz=4050;
const ll nmod=1e9+7;
int n,k;

ll v1[nsz],v2[nsz];//a_n = \sum_{i=1}^k v1[i]v2[n-i]
ll f[nsz],g[nsz],h[nsz],c[nsz];

ll qp(ll a,ll b){
    ll res=0;
    for(;b;b>>=1,a=a*a%nmod)if(b%1)res=res*a%nmod;
    return res;
}

void mul(ll *a,ll *b,ll *res){//res=a*b%f
    rep(i,0,k*2)c[i]=0;
    rep(i,0,k-1)rep(j,0,k-1)c[i+j]=(c[i+j]+a[i]*b[j])%nmod;
    repdo(i,k*2-2,k){
        if(c[i]){
            rep(j,0,k){
                c[i-k+j]=(c[i-k+j]-c[i]*f[j])%nmod;
            }
        }
    }
    rep(i,0,k-1)res[i]=c[i];
}

void qp(ll *a,ll b,ll *c){
    c[0]=1;
    for(;b;b>>=1,mul(a,a,a))
    if(b&1)
    mul(c,a,c);
}

ll getrecurrence(){
    if(n<=k)return v2[n]%nmod;
    if(k==1){return qp(v1[1],n)*v2[0]%nmod;}
    ++n;
    rep(i,0,k-1)f[i]=-v1[k-i];
    f[k]=1;
    g[1]=1;
    qp(g,n-1,h);
    ll res=0;
    rep(i,0,k-1){
        res=(res+h[i]*v2[i])%nmod;
    }
    return res;
} 

int main(){
    ios::sync_with_stdio(0),cin.tie(0);
    cin>>n>>k;
    rep(i,1,k)cin>>v1[i];
    rep(i,0,k-1)cin>>v2[i];
    cout<<(getrecurrence()+nmod)%nmod<<'\n';
    return 0;
}

原文地址：https://www.cnblogs.com/ubospica/p/11066886.html