library

This documentation is automatically generated by online-judge-tools/verification-helper

---

Project maintained by tko919 Hosted on GitHub Pages — Theme by mattgraham

Black Box Linear Algebra
(Math/bbla.hpp)

• View this file on GitHub
• Last update: 2025-07-02 00:02:32+09:00
• Include: #include "Math/bbla.hpp"

Depends on

• Berlekamp Massey Algorithm (FPS/berlekampmassey.hpp)
• Random (Utility/random.hpp)

Verified with

• Verify/LC_sparse_matrix_det.test.cpp
• Verify/YUKI_1112.test.cpp
• Verify/YUKI_310.test.cpp

Code

#pragma once

#include "FPS/berlekampmassey.hpp"

#include "Utility/random.hpp"


template <typename T> Poly<T> RandPoly(int n) {
    Poly<T> ret(n);
    for (auto &x : ret)
        x = Random::get(1, T::get_mod() - 1);
    return ret;
}
template <typename T> struct SparseMatrix {
    vector<T> base;
    vector<map<int, T>> extra;
    SparseMatrix(int n, T v = 0) : base(n, v), extra(n) {}
    int size() const {
        return base.size();
    }
    inline void add(int i, int j, T x) {
        extra[i][j] += x;
    }
    friend Poly<T> operator*(const SparseMatrix<T> &A, const Poly<T> &b) {
        int n = A.size();
        Poly<T> ret(n);
        T sum;
        for (auto &v : b)
            sum += v;
        rep(i, 0, n) {
            T add = sum;
            for (auto &[j, v] : A.extra[i]) {
                ret[i] += v * b[j];
                add -= b[j];
            }
            ret[i] += add * A.base[i];
        }
        return ret;
    }
    void mul(int i, T x) {
        base[i] *= x;
        for (auto &[_, v] : extra[i])
            v *= x;
    }
};

template <typename T> Poly<T> MinPolyforVector(const vector<Poly<T>> &b) {
    int n = b.size(), m = b[0].size();
    Poly<T> base = RandPoly<T>(m), a(n);
    rep(i, 0, n) rep(j, 0, m) a[i] += base[j] * b[i][j];
    return Poly<T>(BerlekampMassey(a));
}
template <typename T> Poly<T> MinPolyforMatrix(const SparseMatrix<T> &A) {
    int n = A.size();
    Poly<T> base = RandPoly<T>(n);
    vector<Poly<T>> b(n * 2 + 1);
    rep(i, 0, n * 2 + 1) b[i] = base, base = A * base;
    return MinPolyforVector(b);
}
template <typename T>
Poly<T> FastPow(const SparseMatrix<T> &A, Poly<T> b, ll t) {
    int n = A.size();
    auto mp = MinPolyforMatrix(A).rev();
    Poly<T> cs({T(1)}), base({T(0), T(1)});
    while (t) {
        if (t & 1) {
            cs *= base;
            cs %= mp;
        }
        base = base.square();
        base %= mp;
        t >>= 1;
    }
    Poly<T> ret(n);
    for (auto &c : cs)
        ret += b * c, b = A * b;
    return ret;
}
template <typename T> T FastDet(const SparseMatrix<T> &A) {
    int n = A.size();
    for (;;) {
        Poly<T> d = RandPoly<T>(n);
        SparseMatrix<T> AD = A;
        rep(i, 0, n) AD.mul(i, d[i]);
        auto mp = MinPolyforMatrix(AD);
        if (mp.back() == 0)
            return 0;
        if (int(mp.size()) != n + 1)
            continue;
        T ret = mp.back(), base = 1;
        if (n & 1)
            ret = -ret;
        for (auto &v : d)
            base *= v;
        return ret / base;
    }
}

/**
 * @brief Black Box Linear Algebra
 */

#line 2 "Math/bbla.hpp"

#line 2 "FPS/berlekampmassey.hpp"

template <typename T> vector<T> BerlekampMassey(vector<T> &a) {
    int n = a.size();
    T d = 1;
    vector<T> b(1), c(1);
    b[0] = c[0] = 1;
    rep(j, 1, n + 1) {
        int l = c.size(), m = b.size();
        T x = 0;
        rep(i, 0, l) x += c[i] * a[j - l + i];
        b.push_back(0);
        m++;
        if (x == 0)
            continue;
        T coeff = -x / d;
        if (l < m) {
            auto tmp = c;
            c.insert(c.begin(), m - l, 0);
            rep(i, 0, m) c[m - 1 - i] += coeff * b[m - 1 - i];
            b = tmp;
            d = x;
        } else
            rep(i, 0, m) c[l - 1 - i] += coeff * b[m - 1 - i];
    }
    reverse(ALL(c));
    return c;
}

/**
 * @brief Berlekamp Massey Algorithm
 */
#line 2 "Utility/random.hpp"

namespace Random {
mt19937_64 randgen(chrono::steady_clock::now().time_since_epoch().count());
using u64 = unsigned long long;
u64 get() {
    return randgen();
}
template <typename T> T get(T L) { // [0,L]

    return get() % (L + 1);
}
template <typename T> T get(T L, T R) { // [L,R]

    return get(R - L) + L;
}
double uniform() {
    return double(get(1000000000)) / 1000000000;
}
string str(int n) {
    string ret;
    rep(i, 0, n) ret += get('a', 'z');
    return ret;
}
template <typename Iter> void shuffle(Iter first, Iter last) {
    if (first == last)
        return;
    int len = 1;
    for (auto it = first + 1; it != last; it++) {
        len++;
        int j = get(0, len - 1);
        if (j != len - 1)
            iter_swap(it, first + j);
    }
}
template <typename T> vector<T> select(int n, T L, T R) { // [L,R]

    if (n * 2 >= R - L + 1) {
        vector<T> ret(R - L + 1);
        iota(ALL(ret), L);
        shuffle(ALL(ret));
        ret.resize(n);
        return ret;
    } else {
        unordered_set<T> used;
        vector<T> ret;
        while (SZ(used) < n) {
            T x = get(L, R);
            if (!used.count(x)) {
                used.insert(x);
                ret.push_back(x);
            }
        }
        return ret;
    }
}

void relabel(int n, vector<pair<int, int>> &es) {
    shuffle(ALL(es));
    vector<int> ord(n);
    iota(ALL(ord), 0);
    shuffle(ALL(ord));
    for (auto &[u, v] : es)
        u = ord[u], v = ord[v];
}
template <bool directed, bool multi, bool self>
vector<pair<int, int>> genGraph(int n, int m) {
    vector<pair<int, int>> cand, es;
    rep(u, 0, n) rep(v, 0, n) {
        if (!self and u == v)
            continue;
        if (!directed and u > v)
            continue;
        cand.push_back({u, v});
    }
    if (m == -1)
        m = get(SZ(cand));
    // chmin(m, SZ(cand));

    vector<int> ord;
    if (multi)
        rep(_, 0, m) ord.push_back(get(SZ(cand) - 1));
    else {
        ord = select(m, 0, SZ(cand) - 1);
    }
    for (auto &i : ord)
        es.push_back(cand[i]);
    relabel(n, es);
    return es;
}
vector<pair<int, int>> genTree(int n) {
    vector<pair<int, int>> es;
    rep(i, 1, n) es.push_back({get(i - 1), i});
    relabel(n, es);
    return es;
}
}; // namespace Random


/**
 * @brief Random
 */
#line 5 "Math/bbla.hpp"

template <typename T> Poly<T> RandPoly(int n) {
    Poly<T> ret(n);
    for (auto &x : ret)
        x = Random::get(1, T::get_mod() - 1);
    return ret;
}
template <typename T> struct SparseMatrix {
    vector<T> base;
    vector<map<int, T>> extra;
    SparseMatrix(int n, T v = 0) : base(n, v), extra(n) {}
    int size() const {
        return base.size();
    }
    inline void add(int i, int j, T x) {
        extra[i][j] += x;
    }
    friend Poly<T> operator*(const SparseMatrix<T> &A, const Poly<T> &b) {
        int n = A.size();
        Poly<T> ret(n);
        T sum;
        for (auto &v : b)
            sum += v;
        rep(i, 0, n) {
            T add = sum;
            for (auto &[j, v] : A.extra[i]) {
                ret[i] += v * b[j];
                add -= b[j];
            }
            ret[i] += add * A.base[i];
        }
        return ret;
    }
    void mul(int i, T x) {
        base[i] *= x;
        for (auto &[_, v] : extra[i])
            v *= x;
    }
};

template <typename T> Poly<T> MinPolyforVector(const vector<Poly<T>> &b) {
    int n = b.size(), m = b[0].size();
    Poly<T> base = RandPoly<T>(m), a(n);
    rep(i, 0, n) rep(j, 0, m) a[i] += base[j] * b[i][j];
    return Poly<T>(BerlekampMassey(a));
}
template <typename T> Poly<T> MinPolyforMatrix(const SparseMatrix<T> &A) {
    int n = A.size();
    Poly<T> base = RandPoly<T>(n);
    vector<Poly<T>> b(n * 2 + 1);
    rep(i, 0, n * 2 + 1) b[i] = base, base = A * base;
    return MinPolyforVector(b);
}
template <typename T>
Poly<T> FastPow(const SparseMatrix<T> &A, Poly<T> b, ll t) {
    int n = A.size();
    auto mp = MinPolyforMatrix(A).rev();
    Poly<T> cs({T(1)}), base({T(0), T(1)});
    while (t) {
        if (t & 1) {
            cs *= base;
            cs %= mp;
        }
        base = base.square();
        base %= mp;
        t >>= 1;
    }
    Poly<T> ret(n);
    for (auto &c : cs)
        ret += b * c, b = A * b;
    return ret;
}
template <typename T> T FastDet(const SparseMatrix<T> &A) {
    int n = A.size();
    for (;;) {
        Poly<T> d = RandPoly<T>(n);
        SparseMatrix<T> AD = A;
        rep(i, 0, n) AD.mul(i, d[i]);
        auto mp = MinPolyforMatrix(AD);
        if (mp.back() == 0)
            return 0;
        if (int(mp.size()) != n + 1)
            continue;
        T ret = mp.back(), base = 1;
        if (n & 1)
            ret = -ret;
        for (auto &v : d)
            base *= v;
        return ret / base;
    }
}

/**
 * @brief Black Box Linear Algebra
 */

Back to top page