Berlekamp–Massey 算法求最短递推式 & Bostan-Mori 算法数列远端求值 / ARIS0_0 - 20

ARIS0_0 · 2026-04-28 11:21:49 · 题解

本文将介绍 Berlekamp-Massey 算法求最短递推式和 Bostan-Mori 算法数列远端求值。

Berlekamp-Massey 算法

给长为 n 的序列 a_0, a_1, \dots, a_{n - 1}，求长度为 m 的序列 f_0, f_1, \dots, f_{m - 1}，使得 \displaystyle \forall i \ge m, a_i = \sum_{j = 0}^{m - 1} f_ja_{i - j - 1}。

考虑增量构造：按照 i 从小到大的顺序依次考虑，如果 a_i 满足当前线性递推式，则不对 f 做任何改变；否则考虑修改 f。初始时，显然有 m = 0, f = \{\}。

假设我们考虑到 a_i 时，原来的 f 已经不满足当前线性递推式，我们需要构造新的 m_{\text{new}} 和 f_{\text{new}}。

令 \displaystyle d_i = a_i - \sum_{j = 0}^{m - 1} f_ja_{i - j - 1}。

如果 i 是第一个需要对递推系数进行修改的位置

说明 \forall j < i, a_j = 0，那么直接令 m_{\text{new}} = i + 1, f_0 = f_1 = \cdots = f_{m_{\text{new}} - 1} = 0 即可。显然这是一个合法的最短线性递推式。注意那个 i + 1 是因为我们是 0-index。

如果 i 不是第一个需要对递推系数进行修改的位置

设上一个需要对递推系数进行修改的位置为 k。

考虑构造一个序列 g = \{g_0, g_1, \dots, g_{m' - 1}\} 使得 \displaystyle \forall l \in [m', i), \sum_{j = 0}^{m' - 1} g_ja_{l - j - 1} = 0，且 \displaystyle \sum_{j = 0}^{m' - 1} g_ja_{i - j - 1} = d_i，那么可以令 m_{\text{new}} = \max(m, m'), f_{\text{new}} = f + g。

设在考虑完 a_{k - 1} 之后序列的最小线性递推系数序列为 h，最小线性递推式长度为 m_h，那么 g 可以这样构造：g_0 = g_1 = \cdots = g_{i - k - 2} = 0, g_{i - k - 1} = \dfrac{d_i}{d_k}, g_{i - k} = -\dfrac{d_i}{d_k}h_0, g_{i - k + 1} = -\dfrac{d_i}{d_k}h_1, g_{i - k + 2} = -\dfrac{d_i}{d_k}h_2, \dots, g_{i - k + m_h - 1} = -\dfrac{d_i}{d_k}h_{m_h - 1}，这样构造 m' = i - k + m_h。

验证此构造是否对于 l \in [m', i) 成立：

由于 \displaystyle \forall l \in [0, k), a_l = \sum_{j = 0}^{m_h - 1} h_ja_{l - j - 1}，且 \forall l \in [m', i), l - i + k < k，所以有上式等于 0。

再验证此构造是否在 i 处成立：

也成立。

令 m_{\text{new}} = \max(m, m'), f_{\text{new}} = f + g，得到新的递推式。可见，Berlekamp-Massey 算法的时间复杂度为 O(nm)。

一个实现细节：如果 f_{new} 的长度等于 f 的长度，那么直接把 f 更新为 f_{new} 即可，我们需要“假装”上一个需要对递推系数进行修改的位置仍然是 k 而不是 i。否则，如果有连续的三个位置 i - 2, i - 1, i，它们的最短递推式长度相同为 m 但系数不同，那么在考虑第三个位置 i 的时候，它对应的 k 就是 i - 1，对应的 k - 1 就是 i - 2，算出来的最短递推式长度就是 1 + m，其中 1 是头插的 \dfrac{d_i}{d_k}，而 m 是 -\dfrac{d_i}{d_k}h，就错了。

常系数齐次线性递推数列表示为有理函数

给定序列 a = \{a_0, a_1, a_2, \dots\} 满足 \displaystyle \forall i \ge m, a_i = \sum_{j = 0}^{m - 1} c_ja_{i - j - 1}，把 a 表示为 \dfrac{p(x)}{q(x)} 的形式，其中 p(x) 和 q(x) 是两个有限次多项式。

考虑 q(x) = 1 - c_0x - c_1x^2 - c_2x^3 - \dots - c_{m - 1}x^m，那么对于 p(x) = q(x)A(x) 和一个 i \in [m, +\infin)，有 \displaystyle p_i = \sum_{j = 0}^{m} q_ja_{i - j} = a_i - \sum_{j = 0}^{m - 1}c_ja_{i - j - 1} = 0，也就是说 p(x) 的次数最高是 m - 1 的。

所以，构造 \displaystyle q(x) = 1 - c_0x - c_1x^2 - c_2x^3 - \dots - c_{m - 1}x^m, p(x) = \left(q(x)\sum_{i = 0}^{m - 1} a_ix^i\right) \bmod x^m 是合法的。

Bostan-Mori 算法

给多项式 p(x) 和 q(x)，求 [x^n]\dfrac{p(x)}{q(x)}，其中 \deg p \le 10^5, \deg q \le 10^5, n \le 10^{18}。

考虑把 \dfrac{p(x)}{q(x)} 的上下同乘 q(-x)，变成 \dfrac{p(x)q(-x)}{q(x)q(-x)}。

这时，分母是偶函数，可以表示为 q_{\text{new}}(x^2) = q(x)q(-x)；又由于偶函数的逆也是偶函数，所以第 n 项的系数，就由分子决定了。

考虑把分子奇偶项拆开，表示为 p(x)q(-x) = p_0(x^2) + xp_1(x^2)。如果 n 是奇数，那么令 p_{\text{new}}(x) = p_1(x)；否则令 p_{\text{new}}(x) = p_0(x)，问题转化为求 \left[x^{\left\lfloor\frac{n}{2}\right\rfloor}\right]\dfrac{p_{\text{new}}(x)}{q_{\text{new}}(x)}。

这样，每次 n 减半，q 的长度不变，p 的长度每次向 q 趋近，如果 \deg p 和 \deg q 差距不大的话，复杂度可以表示为 O(\deg q \log \deg q \log n)。

Bostan-Mori 数列远端求值

给定序列 a = \{a_0, a_1, a_2, \dots\} 满足 \displaystyle \forall i \ge m, a_i = \sum_{j = 0}^{m - 1} c_ja_{i - j - 1}，其中 m \le 5000，和一个正整数 n \le 10^{18}，求 a_n。

把 a 的生成函数表示为有理函数，然后做 Bostan-Mori 即可。时间复杂度 O(m \log m \log n)。

#include <cstdio>
#include <vector>
#define bit_ceil(x) ((x) <= 1 ? 1 : (1 << (32 - __builtin_clz((x) - 1))))
typedef std::vector<int> poly;
namespace FastIn
{
    const int S = 4e6 + 7;
    char *p1, *p2, buf[S];
    #define gc() (p1 == p2 && (p2 = (p1 = buf) + fread(buf, 1, S, stdin), p1 == p2) ? EOF : *p1++)
    inline void read(int &x)
    {
        x = 0;
        char c = gc();
        while (c < '0' || c > '9')
            c = gc();
        for (; c >= '0' && c <= '9'; c = gc())
            x = (x << 3) + (x << 1) + (c & 15);
    }
    #undef gc
}
using FastIn::read;
const int mod = 998244353;
inline int M(const int &x) { return x >= mod ? x - mod : x; }
inline int qpow(int a, int b)
{
    int res = 1;
    for (; b; a = (long long)(a) * a % mod, b >>= 1)
        (b & 1) && (res = (long long)(res) * a % mod);
    return res;
}
inline void NTT(poly &a, const int &rvs)
{
    static std::vector<int> rev;
    int n = a.size(); (int)(rev.size()) < n && (rev.resize(n), true);
    for (int i = 0; i < n; ++i)
        if (i < (rev[i] = (rev[i >> 1] >> 1) | ((i & 1) * (n >> 1))))
            std::swap(a[i], a[rev[i]]);
    for (int len = 1, W; len < n; len <<= 1)
    {
        W = qpow(rvs == 1 ? 3 : 332748118, (mod - 1) / (len << 1));
        for (int i = 0; i < n; i += (len << 1))
            for (int j = i, f0, f1, w = 1; j < i + len; ++j, w = (long long)(w) * W % mod)
                f0 = a[j], f1 = (long long)(w) * a[j + len] % mod, a[j] = M(f0 + f1), a[j + len] = M(f0 + mod - f1);
    }
    if (rvs == -1)
        for (int i = 0, invn = qpow(n, mod - 2); i < n; ++i)
            a[i] = (long long)(a[i]) * invn % mod;
}
inline poly operator*(const poly &A, const poly &B)
{
    int n = A.size() + B.size() - 1, len = bit_ceil(n);
    poly a = A, b = B;
    a.resize(len), b.resize(len), NTT(a, 1), NTT(b, 1);
    for (int i = 0; i < len; ++i)
        a[i] = (long long)(a[i]) * b[i] % mod;
    return NTT(a, -1), a.resize(n), a;
}
inline poly operator*=(poly &a, const poly &b) { return a = a * b; }
inline poly operator+(const poly &A, const poly &B)
{
    poly a = A, b = B;
    a.resize(std::max(A.size(), B.size())), b.resize(std::max(A.size(), B.size()));
    for (int i = 0; i < (int)(a.size()); ++i)
        a[i] = M(a[i] + b[i]);
    return a;
}
namespace BM
{
    int n, a[10007], d[10007], p[10007];
    poly f[10007];
    inline void solve()
    {
        (d[0] = a[0]) ? (f[0] = (poly){0}, p[0] = 0) : (p[0] = -1);
        for (int i = 1, k; i < n; ++i)
        {
            if (~p[i - 1])
                for (int j = 0; j < (int)(f[p[i - 1]].size()); ++j)
                    d[i] = (d[i] + (long long)(f[p[i - 1]][j]) * a[i - j - 1]) % mod;
            if (!(d[i] = M(a[i] + mod - d[i])))
            {
                p[i] = p[i - 1];
                continue;
            }
            else if (p[i - 1] == -1)
            {
                f[p[i] = i].resize(i + 1);
                continue;
            }
            int c = (long long)(d[i]) * qpow(d[k = p[i - 1]], mod - 2) % mod;
            poly g;
            for (int j = 0; j < i - k - 1; ++j)
                g.push_back(0);
            g.push_back(c);
            if ((~(k - 1)) && (~(p[k - 1])))
                for (const int &j : f[p[k - 1]])
                    g.push_back((long long)(mod - c) * j % mod);
            poly tmp = f[k] + g;
            if (f[k].size() == tmp.size())
                f[p[i] = k] = tmp;
            else
                f[p[i] = i] = tmp;
        }
    }
}
inline int BostanMori(poly p, poly q, int k)
{
    while (k)
    {
        int n = p.size(), m = q.size();
        poly tmp = q;
        for (int i = 1; i < m; i += 2)
            tmp[i] = M(mod - q[i]);
        p *= tmp, q *= tmp, n = (p.size() + ((k & 1) ^ 1)) >> 1;
        for (int i = 0; i < n; ++i)
            p[i] = p[(i << 1) + (k & 1)];
        for (int i = 0; i < m; ++i)
            q[i] = q[i << 1];
        p.resize(n), q.resize(m), k >>= 1;
    }
    return (long long)(p[0]) * qpow(q[0], mod - 2) % mod;
}
int n, m, ans, a[10007];
poly f;
int main()
{
    read(n), read(m);
    for (int i = 0; i < n; ++i)
        read(a[i]), BM::a[i] = a[i];
    BM::n = n, BM::solve();
    if (BM::p[BM::n - 1] == -1)
        return putchar('0'), putchar('\n'), putchar('0'), 0;
    poly p, q = (poly){1};
    for (int i = 0; i < (int)(BM::f[BM::p[BM::n - 1]].size()); ++i)
        p.push_back(a[i]);
    for (const int &i : BM::f[BM::p[BM::n - 1]])
        q.push_back(M(mod - i));
    p *= q, p.resize(BM::f[BM::p[BM::n - 1]].size());
    for (const int &i : BM::f[BM::p[BM::n - 1]])
        printf("%d ", i);
    printf("\n%d", BostanMori(p, q, m));
    return 0;
}
// 不曾走过炫暇盛宴的绝顶，唯能步入终幕遗信的黄昏
// 20