对于 $n$的范围如此巨大得不正常的题目我们通常会考虑用矩阵快速幂解决
能用矩阵快速幂解决的话需要满足状态转移方程是个 “齐次线性常系数递推”
何谓 “齐次线性常系数递推” 呢？
比如这一题，当 $k$已知，设 $f[i]$表示有 $i$个点时的答案
若存在一个数列 $p$，满足对于 $\forall f[i](i > |p|)$都有 $f[i] = \sum _ {j = 1} ^ {|p|} p[j] \times f[i – j]$，则称 $p$为 $f$的递推式，$f$也满足齐次线性常系数递推关系
比如等比数列 $\{ 1, 2, 4, 8, 16 … \}$的递推式就是 $\{ 2 \}$，斐波那契数列的递推式就是 $\{1, 1\}$
此题正解是状态压缩然后手推递推式，但是似乎非常麻烦
有没有什么办法轻松获取递推式呢？
答案是肯定的，Berlekamp-Massey 算法就是用来计算递推式的
BM 算法可以根据一个数列计算其递推式（需要满足数列长度大于等于二倍递推式长度），复杂度 $O(n ^ 2)$（$n$为给定数列长度）
（比如给定 $\{ 1, 1, 2, 3 \}$BM 算法就能算出 $\{ 1, 1 \}$，当然为了保险起见最好能算多少项就算多少项）
关于 BM 的理解我认为 Boshi 大佬写的挺好的：https://www.mina.moe/archives/7500
也许以后有空我会发一篇讲 BM 算法的文章吧（咕咕咕警告.jpg）
另外也可以参考 fjzzq 大佬的博客 https://www.cnblogs.com/zzqsblog/p/6877339.html（注意他的文章中被删除线划去的部分是错误的，即最近的不一定最优）
我个人认为自己的 BM 写的还是挺好的 QwQ，这是我第四个版本的 BM 算法代码了：
#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;
}
就这一题而言只需要用矩阵树定理先算出 $f[1], f[2], f[3] …$的值，然后把这个数列丢给 BM，BM 算出其递推式后转成矩阵，快速幂跑出答案即可，是不是非常无脑，哈哈哈~
本题中递推式大概是 $50$项的样子（实测不到 $50$项），因此要打一个长度为 $100$的表，即 $f[1], f[2], f[3], …, f[100]$。
高斯消元求行列式的值复杂度是 $O(n ^ 3)$的，因此打表复杂度 $O(n ^ 4)$，看起来挺悬的其实跑得还挺快的
如果是在线打表怕 TLE 的话可以先打好表写代码里 QvQ
代码：
#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;
}
0 条评论