题目连接

对于 $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;
}
分类: 文章

Remmina

No puzzle that couldn't be solved.

0 条评论

发表评论

电子邮件地址不会被公开。 必填项已用*标注