Linear Recurrence (Berlekamp Massey, K-th term)

BerlekampMassey.cpp

vector<ll> BerlekampMassey(vector<ll> s) {
    int n = s.size(), L = 0, m = 0;
    vector<ll> C(n), B(n), T;
    C[0] = B[0] = 1;

    ll b = 1;
    for (int i = 0; i < n; i++) {
        m++;
        ll d = s[i] % MOD;
        for (int j = 1; j < L + 1; j++)
            d = (d + C[j] * s[i - j]) % MOD;
        if (!d)
            continue;
        T = C;
        ll coef = d * binExp(b, MOD - 2) % MOD;
        for (int j = m; j < n; j++)
            C[j] = (C[j] - coef * B[j - m]) % MOD;
        if (2 * L > i)
            continue;
        L = i + 1 - L;
        B = T, b = d, m = 0;
    }
    C.resize(L + 1);
    C.erase(C.begin());
    for (auto &x : C)
        x = (MOD - x) % MOD;
    return C;
}

Kth-term.cpp

typedef vector<ll> Poly;
ll linearRec(Poly S, Poly tr, ll k) {
    int n = S.size();
    auto combine = [&](Poly a, Poly b) {
        Poly res(n * 2 + 1);
        for (int i = 0; i < n + 1; i++)
            for (int j = 0; j < n + 1; j++)
                res[i + j] = (res[i + j] + a[i] * b[j]) % MOD;
        for (int i = 2 * n; i > n; --i)
            for (int j = 0; j < n; j++)
                res[i - 1 - j] = (res[i - 1 - j] + res[i] * tr[j]) % MOD;
        res.resize(n + 1);
        return res;
    };
    Poly pol(n + 1), e(pol);
    pol[0] = e[1] = 1;
    for (++k; k; k /= 2) {
        if (k % 2)
            pol = combine(pol, e);
        e = combine(e, e);
    }
    ll res = 0;
    for (int i = 0; i < n; i++)
        res = (res + pol[i + 1] * S[i]) % MOD;
    return res;
}

Berlekamp Massey recovers any n-order linear recurrence from the first >= 2n terms of the sequence. Useful for guessing linear recurrences after bruteforcing the first terms. Should work on any field, but numerical stability for floats is not guaranteed.

Usage: BerlekampMassey({0, 1, 1, 3, 5, 11}) // {1, 2}

Kth-term finds the k-th term of an n-order linear recurrence. Useful together with Berlekamp Massey. Faster than matrix exponentiation: O(n^2 log k)

Usage: linearRec({0,1}, {1,1}, k) // k-th Fibonacci number

Combined usage:

vector<ll> x; // some brute forced linear-recurrence S
vector<ll> tr = BerlekampMassey(x);
x.resize(tr.size());
linearRec(x, tr, k) // Finds k-th term of the linear recurrence represented by x.