library
This documentation is automatically generated by online-judge-tools/verification-helper
Black Box Linear Algebra
(Math/bbla.hpp)
- View this file on GitHub
- Last update: 2024-04-26 03:18:17+09:00
- Include:
#include "Math/bbla.hpp"
Depends on
Verified with
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)).rev();
}
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];
}
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 simple> vector<pair<int, int>> genGraph(int n) {
vector<pair<int, int>> cand, es;
rep(u, 0, n) rep(v, 0, n) {
if (simple and u == v)
continue;
if (!directed and u > v)
continue;
cand.push_back({u, v});
}
int m = get(SZ(cand));
vector<int> ord;
if (simple)
ord = select(m, 0, SZ(cand) - 1);
else {
rep(_, 0, m) ord.push_back(get(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)).rev();
}
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
*/