BM 算法可以根据一个数列计算其递推式（需要满足数列长度大于等于二倍递推式长度），复杂度 $O(n ^ 2)$（$n$为给定数列长度）

（比如给定 $\{ 1, 1, 2, 3 \}$BM 算法就能算出 $\{ 1, 1 \}$，当然为了保险起见最好能算多少项就算多少项）

#include <bits/stdc++.h>

#define NS (3005)
#define MOD (1000000007)
#define Inv(a) (qpow((a), MOD - 2))
#define pls(a, b) ((a) + (b) < MOD ? (a) + (b) : (a) + (b) - MOD)
#define mns(a, b) ((a) - (b) < 0 ? (a) - (b) + MOD : (a) - (b))
#define mul(a, b) (1ll * (a) * (b) % MOD)

using namespace std;

template<typename _Tp> inline void IN(_Tp& dig)
{
char c; bool flag = 0; dig = 0;
while (c = getchar(), !isdigit(c)) if (c == '-') flag = 1;
while (isdigit(c)) dig = dig * 10 + c - '0', c = getchar();
if (flag) dig = -dig;
}

int qpow(int a, int b)
{
int res = 1;
while (b)
{
if (b & 1) res = mul(res, a);
a = mul(a, a), b >>= 1;
}
return res;
}

int n, A[NS];

int f[NS], fs;

void BM()
{
static int g[NS], gs, t[NS], ts;
memset(f, 0, sizeof(f)), memset(g, 0, sizeof(g)), fs = gs = 0;
int dt, ld, p, k;
for (int i = 1; i <= n; i += 1)
{
dt = A[i];
for (int j = 1; j <= fs; j += 1) dt = mns(dt, mul(f[j], A[i - j]));
if (!dt) continue;
if (!fs) {fs = p = i, ld = dt; continue;}
k = mns(0, mul(dt, Inv(ld))), ts = max(fs, i - p + gs);
memset(t, 0, sizeof(int) * ts), t[i - p] = mns(0, k);
for (int j = 1; j <= gs; j += 1) t[i - p + j] = mul(g[j], k);
for (int j = 1; j <= fs; j += 1) t[j] = pls(t[j], f[j]);
if (i - p + gs >= fs)
gs = fs, memmove(g + 1, f + 1, sizeof(int) * gs), p = i, ld = dt;
fs = ts, memmove(f + 1, t + 1, sizeof(int) * fs);
}
}

int main(int argc, char const* argv[])
{
IN(n);
for (int i = 1; i <= n; i += 1) IN(A[i]);
BM(), printf("%d\n", fs);
for (int i = 1; i <= fs; i += 1) printf("%d ", f[i]);
putchar(10);
return 0;
}


#include <bits/stdc++.h>

#define TS (100)
#define MS (50)
#define MOD (65521)
#define Inv(a) (qpow((a), MOD - 2))
#define pls(a, b) ((a) + (b) < MOD ? (a) + (b) : (a) + (b) - MOD)
#define mns(a, b) ((a) - (b) < 0 ? (a) - (b) + MOD : (a) - (b))
#define mul(a, b) (1ll * (a) * (b) % MOD)

using namespace std;

typedef long long LL;

int k, D[TS + 5][TS + 5], A[TS + 5];

LL n;

int Gs(int p)
{
int res = 1;
for (int i = 1; i <= p; i += 1)
{
for (int j = i + 1; j <= p; j += 1)
while (D[j][i])
{
int tmp = D[i][i] / D[j][i];
for (int k = i; k <= p; k += 1)
D[i][k] = mns(D[i][k], mul(tmp, D[j][k]));
swap(D[i], D[j]), res = mns(0, res);
}
res = mul(res, D[i][i]);
}
return res;
}

void Get_Tab()
{
for (int q = 1; q <= TS; q += 1)
{
int p = k + q - 1;
for (int i = 1; i <= p; i += 1)
for (int j = 1; j <= p; j += 1)
D[i][j] = 0;
for (int i = 1; i <= p; i += 1)
for (int j = max(1, i - k); j <= min(p, i + k); j += 1)
if (i != j) D[i][j] = mns(0, 1), D[i][i]++;
A[q] = Gs(p - 1);
}
}

int qpow(int a, int b)
{
int res = 1;
while (b)
{
if (b & 1) res = mul(res, a);
a = mul(a, a), b >>= 1;
}
return res;
}

int f[TS + 5], fs;

void BM()
{
static int g[TS + 5], gs, t[TS + 5], ts;
memset(f, 0, sizeof(f)), memset(g, 0, sizeof(g)), fs = gs = 0;
int dt, ld, p, k;
for (int i = 1; i <= TS; i += 1)
{
dt = A[i];
for (int j = 1; j <= fs; j += 1) dt = mns(dt, mul(f[j], A[i - j]));
if (!dt) continue;
if (!fs) {fs = p = i, ld = dt; continue;}
k = mns(0, mul(dt, Inv(ld))), ts = max(fs, i - p + gs);
memset(t, 0, sizeof(int) * ts), t[i - p] = mns(0, k);
for (int j = 1; j <= gs; j += 1) t[i - p + j] = mul(g[j], k);
for (int j = 1; j <= fs; j += 1) t[j] = pls(t[j], f[j]);
if (i - p + gs >= fs)
gs = fs, memmove(g + 1, f + 1, sizeof(int) * gs), p = i, ld = dt;
fs = ts, memmove(f + 1, t + 1, sizeof(int) * fs);
}
}

int m, M[MS][MS], R[MS][MS];

void M_mul(int (&a)[MS][MS], int (&b)[MS][MS])
{
static int tmp[MS][MS];
memset(tmp, 0, sizeof(tmp));
for (register int k = 0; k < m; ++k)
for (register int i = 0; i < m; ++i)
for (register int j = 0; j < m; ++j)
tmp[i][j] = pls(tmp[i][j], mul(a[i][k], b[k][j]));
memmove(a, tmp, sizeof(a));
}

void M_pow(LL b)
{
for (int i = 0; i < m; i += 1) R[i][i] = 1;
while (b)
{
if (b & 1) M_mul(R, M);
M_mul(M, M), b >>= 1;
}
}

int main(int argc, char const* argv[])
{
scanf("%d%lld", &k, &n), Get_Tab(), BM(), m = fs;
for (int i = 0; i < m - 1; i += 1) M[i + 1][i] = 1;
for (int i = 0; i < m; i += 1) M[i][m - 1] = f[m - i];
M_pow(n - k);
int ans = 0;
for (int i = 0; i < m; i += 1) ans = pls(ans, mul(A[i + 1], R[i][0]));
printf("%d\n", ans);
return 0;
}


#### Remmina

No puzzle that couldn't be solved.