[模板] Berlekamp–Massey 算法
一、题目
二、解法
首先说一下这东西能干什么：在你遇到无法下手的计数题时，不妨默写下 \(\tt BM\) 算法，然后丢前几项进去，他就会给你返回一个线性递推式，然后我们依照信仰使用这个递推式即可。
考虑数列 \(\{a_1,a_2...a_n\}\)，设 \(\{a_1,a_2...a_i\}\) 的最短线性递推为 \(R_i\)，\(a_{i+1}\) 与 \(R_i\) 推算出的下一个数的差值为 \(d_i\)，特别地，\(R_0=\{\}\)
现在考虑如果我们已经知道了 \(R_1,R_2...R_{i-1}\)，如何求出 \(R_i\)，我们分类讨论一下：
- 如果 \(d_{i-1}=0\)，那么 \(R_{i}=R_{i-1}\)
- 如果 \(R_{i-1}=\{\}\)，那么令 \(R_{i}=\{0,0...0\}\)（共有 \(i\) 个 \(0\)），显然可以让除了这一位之外的数都为 \(0\)
- 否则 \(R_i\) 需要在 \(R_{i-1}\) 的基础上修正，考虑用下面的方法修正：
设线性递推 \(R'=\{r'_1,r'_2...r'_k\}\) 满足 \(\sum_{j=1}^kr'_j\cdot a_{w-j}=0\) 对所有 \(k<w<i\) 都成立，并且在 \(w=i\) 时，得到的值是 \(d_{i-1}\)，那么 \(R_i=R_{i-1}+R'\)
考虑构造 \(R'=\{0,0...0,\frac{d_{i-1}}{d_{o}},-\frac{d_{i-1}}{d_o}\times R_o\}\)，其中 \(0\) 的个数是 \(i-o-2\)，在前面补 \(0\) 的原因是我们要还原在 \(o\) 处的情况，考虑带入 \(i\) 会得到 \(d_o\) 的值，我们通过添加系数可以把他改成 \(d_{i-1}\)，其他地方的取值就是 \(0\)
为了是递推式阶数最小，我们可以让 \(o\) 为 \(R_o\) 的阶数\(-o\) 最小的位置 \(o\)，时间复杂度 \(O(n^2)\)
至于后面的快速递推，可以使用 常系数齐次线性递推，这里我们取模和乘法都暴力做就可以通过本题，时间复杂度 \(O(n^2\log m)\)
#include <cstdio>
#include <vector>
#include <iostream>
using namespace std;
const int M = 10005;
const int MOD = 998244353;
#define int long long
#define pb push_back
int read()
{
int x=0,f=1;char c;
while((c=getchar())<'0' || c>'9') {if(c=='-') f=-1;}
while(c>='0' && c<='9') {x=(x<<3)+(x<<1)+(c^48);c=getchar();}
return x*f;
}
int n,m,a[M],f[M],g[M],tmp[M],p[M];
int qkpow(int a,int b)
{
int r=1;
while(b>0)
{
if(b&1) r=r*a%MOD;
a=a*a%MOD;
b>>=1;
}
return r;
}
void add(int &x,int y) {x=(x+y)%MOD;}
void BM(int n,int *a,vector<int> &r)
{
r.clear();vector<int> ls;int w=0,d=0;
for(int i=1;i<=n;i++)
{
int tmp=0;
for(int j=0;j<r.size();j++)
tmp=(tmp+a[i-1-j]*r[j])%MOD;
if((a[i]-tmp)%MOD==0) continue;
if(!w)
{
w=i;d=a[i]-tmp;
for(int j=i;j;j--) r.pb(0);
continue;
}
vector<int> nw=r;
int mul=(a[i]-tmp)*qkpow(d,MOD-2)%MOD;
if(r.size()<ls.size()+i-w)
r.resize(ls.size()+i-w);
add(r[i-w-1],mul);
for(int j=0;j<ls.size();j++)
add(r[i-w+j],-mul*ls[j]%MOD);
if((int)nw.size()-i<(int)ls.size()-w)
ls=nw,w=i,d=a[i]-tmp;
}
}
int calc(vector<int> &r)
{
if(m<n) return a[m];
int k=r.size(),ans=0;p[0]=-1;
for(int i=1;i<=k;i++) p[i]=r[i-1];
f[0]=1;(k>1)?g[1]=1:g[0]=p[1];
auto mul = [&] (int *a,int *b,int *c)
{
for(int i=0;i<2*k;i++) tmp[i]=0;
for(int i=0;i<k;i++)
for(int j=0;j<k;j++)
add(tmp[i+j],a[i]*b[j]%MOD);
for(int i=2*k;i>=k;i--) if(tmp[i])
for(int j=k;j>=0;j--)
add(tmp[i-j],tmp[i]*p[j]%MOD);
for(int i=0;i<k;i++) c[i]=tmp[i];
return 0;
};
while(m)
{
if(m&1) mul(f,g,f);
mul(g,g,g);
m>>=1;
}
for(int i=0;i<k;i++)
ans=(ans+a[i+1]*f[i])%MOD;
return ans;
}
signed main()
{
n=read();m=read();
for(int i=1;i<=n;i++) a[i]=read();
vector<int> r;BM(n,a,r);
for(auto x:r) printf("%lld ",(x+MOD)%MOD);
puts("");
printf("%lld\n",(calc(r)+MOD)%MOD);
}
