[高端操作]常系数线性递推式

对一个常系数线性递推式$$f_n = sum_{i=1}^k a_i imes f_{n-i}$$

矩阵快速幂需要 $O(k^3logn)$ ?

这篇文章将教您在至多为 $O(k^2logn + k^4)$ 时间内搞这个式子

有什么用嘛，可能对我这种省选注定退役的人没啥用，但对于 NOI 及以上比赛想 AK 的同学们可能有用

1.矩阵的特征多项式

矩阵 $A$ 的特征多项式是一个关于 $lambda$ 的函数 $f(lambda)=det(lambda I - A)$，其中 $I$ 为单位矩阵，$det()$ 为行列式

2.Caylay-Hamilton 定理

就一句话，$f(A) = 0$，证明略

3.求矩阵特征多项式的方法

　　1) 消一消　

　　考虑根据定义，需要求行列式，消一下就可以了，复杂度 $O(k^5)$ 或 $O(k^4)$

　　2)插一插

　　考虑到特征多项式是一个 $k$ 次的多项式，我们可以把 $lambda = 0,1,2, cdots ,k$ 时的点值求出来，然后拉格朗日插值，复杂度 $O(k^4)$

4.有啥用

根据小学知识

$$被除数 ÷ 除数 = 商 cdots cdots 余数$$

我们发现，一个次数大于等于 $k$ 的矩阵幂，可以把它除以 $f(A)$

之后我们发现：因为除数等于 $0$ ，所以$余数 = 被除数$

然后我们发现，余数的次数不超过 $k$

于是我们可以用类似快速幂的做法倍增，把 $n$ 次的矩阵幂变成一个不超过 $k$ 次的多项式

其中每一步如果用快速的多项式取模 (FFT) ，是 $O(klogk)$ 的，如果暴力，是 $O(k^2)$ 的

这一步的复杂度是 $O(klogklogn)$ 或者 $O(k^2logn)$

然后就做完了。。。看上去很简单，但我们可以数数，一道完整的题，您要写多少代码

1.多项式相关操作（FFT，求逆，取膜）

2.构造矩阵（dfs / 根据题目情况而定）

3.构造特征多项式（消元 / 拉格朗日插值）

4.倍增的单步转移

5.倍增的全过程

嗯...不是很长？

附一个拉格朗日插值的代码吧

给一个长度为 $n$ 的多项式函数（注意，是长度，所以次数为 $n-1$），求它在 k 处的点值

#include<bits/stdc++.h>
#define LL long long 
using namespace std;
inline int read()
{
    int x=0,f=1; char ch=getchar();
    for(;!isdigit(ch);ch=getchar())if(ch=='-')f=-1;
    for(;isdigit(ch);ch=getchar())x=ch-'0'+(x<<3)+(x<<1);
    return x*f;
}
inline void write(int x)
{
    if(x<0)putchar('-'),x=-x;
    int t=10,len=1; while(t<=x)t*=10,len++;
    while(len--)t/=10,putchar(x/t+48),x%=t;
    putchar(' '); return ;
}
const int maxn = 2010,mod = 998244353;
int n,k,x[maxn],y[maxn];
inline int skr(int x,int t)
{
    int res = 1;
    while(t)
    {
        if(t & 1)res = 1LL * res * x % mod;
        x = 1LL * x * x % mod;
        t = t >> 1;
    }return res;
}
int main()
{
    n = read() - 1,k = read();
    for(int i=0;i<=n;i++)x[i] = read(),y[i] = read();
    int ans = 0;
    for(int i=0;i<=n;i++)
    {
        int t1 = 1,t2 = 1;
        for(int j=0;j<=n;j++)
        {
            if(i != j)
            {
                t1 = 1LL * t1 * (k - x[j]) % mod;
                t2 = 1LL * t2 * (x[i] - x[j]) % mod;
            }
        }
        ans = ((LL)ans + (LL)y[i] * t1 % mod * skr(t2,mod - 2) % mod) % mod;
    }
    ans = ((ans % mod) + mod) % mod;
    cout<<ans;
}