Black Box Linear Algebra
ベクトルと行列の掛け算に\(T(n)\)時間かかるとします.
最小多項式
\(O(n^2 + n T(n))\)
#include "../berlekamp_massey.hpp"
#include <vector>
template<class F>
std::vector<F> find_minimal_polynomial(const std::vector<F>& a) {
std::vector<F> c = berlekamp_massey(a);
c.insert(c.begin(), -F(1));
return c;
}
template<class F, class NonZeroRandGen>
std::vector<F> find_minimal_polynomial_from_vector(int n, const std::vector<std::vector<F>>& a, NonZeroRandGen rnd) {
std::vector<F> u(n);
for(int i = 0; i < n; i++) u[i] = rnd();
std::vector<F> b(a.size(), F(0));
for(int i = 0; i < a.size(); i++) {
for(int j = 0; j < n; j++) {
b[i] += a[i][j] * u[j];
}
}
return find_minimal_polynomial(b);
}
template<class F, class NonZeroRandGen>
std::vector<F> find_minimal_polynomial_from_dense_matrix_pow_b(const std::vector<std::vector<F>>& a, std::vector<F> b, NonZeroRandGen rnd) {
int n = a.size();
std::vector<F> bf;
std::vector<F> u(n);
for(int i = 0; i < n; i++) u[i] = rnd();
std::vector<F> c(n * 2);
for(int i = 0; i < 2 * n; i++) {
for(int j = 0; j < n; j++) {
c[i] += b[j] * u[j];
}
if(i + 1 < 2 * n) {
bf = b;
for(int j = 0; j < n; j++) {
b[j] = F(0);
for(int k = 0; k < n; k++) {
b[j] += a[j][k] * bf[k];
}
}
}
}
return find_minimal_polynomial(c);
}
// fast for dense matrix
template<class F, class NonZeroRandGen>
std::vector<F> find_minimal_polynomial_from_dense_matrix_pow(const std::vector<std::vector<F>>& a, NonZeroRandGen rnd) {
int n = a.size();
std::vector<F> b(n);
for(int i = 0; i < n; i++) b[i] = rnd();
return find_minimal_polynomial_from_dense_matrix_pow_b(a, std::move(b), rnd);
}
#include <tuple>
template<class F, class NonZeroRandGen>
std::vector<F> find_minimal_polynomial_from_sparse_matrix_pow_b(const std::vector<std::tuple<int, int, F>>& a, std::vector<F> b, int n, NonZeroRandGen rnd) {
std::vector<F> bf;
std::vector<F> u(n);
for(int i = 0; i < n; i++) u[i] = rnd();
std::vector<F> c(n * 2);
for(int i = 0; i < 2 * n; i++) {
for(int j = 0; j < n; j++) {
c[i] += b[j] * u[j];
}
if(i + 1 < 2 * n) {
bf = b;
for(int j = 0; j < n; j++) {
b[j] = F(0);
}
for(auto& [j, k, v]: a) {
b[j] += v * bf[k];
}
}
}
return find_minimal_polynomial(c);
}
template<class F, class NonZeroRandGen>
std::vector<F> find_minimal_polynomial_from_sparse_matrix_pow(const std::vector<std::tuple<int, int, F>>& a, int n, NonZeroRandGen rnd) {
std::vector<F> b(n);
for(int i = 0; i < n; i++) b[i] = rnd();
return find_minimal_polynomial_from_sparse_matrix_pow_b(a, std::move(b), n, rnd);
}
行列式
\(O(n^2 + n T(n))\)
#include "minimal_polynomial.hpp"
#include <vector>
#include <tuple>
template<class F, class NonZeroRandGen>
F fast_determinant_dense(std::vector<std::vector<F>> a, NonZeroRandGen rng) {
int n = a.size();
std::vector<F> d(n);
F d_det(1);
for(int i = 0; i < n; i++) {
d[i] = rng();
d_det *= d[i];
}
for(int i = 0; i < n; i++) {
for(int j = 0; j < n; j++) {
a[i][j] = a[i][j] * d[j];
}
}
auto c = find_minimal_polynomial_from_dense_matrix_pow2(a, rng);
F det = c.back() / c.front();
det = n & 1 ? -det : det;
return det / d_det;
}
template<class F, class NonZeroRandGen>
F fast_determinant_sparse(std::vector<std::tuple<int, int, F>> a, int n, NonZeroRandGen rng) {
std::vector<F> d(n);
F d_det(1);
for(int i = 0; i < n; i++) {
d[i] = rng();
d_det *= d[i];
}
for(auto& [i, j, v]: a) {
v *= d[j];
}
auto c = find_minimal_polynomial_from_sparse_matrix_pow2(a, n, rng);
F det = c.back() / c.front();
det = n & 1 ? -det : det;
return det / d_det;
}
行列累乗
\(A^k b\)を求める. fast_kitamasa
を用いて\(O(N^3 + N \log N \log k)\)
#include "../fast_kitamasa.hpp"
#include <vector>
#include <tuple>
template<class F, class FM, class NonZeroRandGen>
std::vector<F> bbla_dense_matrix_pow(const std::vector<std::vector<F>>& a, std::vector<F> b, long long r, NonZeroRandGen rnd) {
int n = a.size();
auto c = find_minimal_polynomial_from_dense_matrix_pow_b(a, b, rnd);
auto d = fast_kitamasa<F, FM>(std::move(c), r);
std::vector<F> ans(n);
std::vector<F> bf;
for(int i = 0; i < d.size(); i++) {
for(int j = 0; j < n; j++) {
ans[j] += d[d.size() - 1 - i] * b[j];
}
if(i + 1 < d.size()) {
bf = b;
for(int j = 0; j < n; j++) {
b[j] = F(0);
for(int k = 0; k < n; k++) {
b[j] += a[j][k] * bf[k];
}
}
}
}
return ans;
}
template<class F, class FM, class NonZeroRandGen>
std::vector<F> bbla_sparse_matrix_pow(const std::vector<std::tuple<int, int, F>>& a, std::vector<F> b, long long r, NonZeroRandGen rnd) {
int n = b.size();
auto c = find_minimal_polynomial_from_sparse_matrix_pow_b(a, b, n, rnd);
auto d = fast_kitamasa<F, FM>(std::move(c), r);
std::vector<F> ans(n);
std::vector<F> bf;
for(int i = 0; i < d.size(); i++) {
for(int j = 0; j < n; j++) {
ans[j] += d[d.size() - 1 - i] * b[j];
}
if(i + 1 < d.size()) {
bf = b;
for(int j = 0; j < n; j++) {
b[j] = F(0);
}
for(auto& [j, k, v]: a) {
b[j] += v * bf[k];
}
}
}
return ans;
}