cp-library
競技プログラミング用 C++ ライブラリ
行列（軽量版）
math/matrix_light.hpp · GitHub で見る
64 KiBサイト · 軽量実装
この項目の軽量実装は LightMatrix<T> です。以下の使用例・API一覧は共有ヘッダーの説明です。
LightMatrix<T> は 行列 と同じ操作を使える軽量実装です。通常版と同じ共通コードを使い、数値型専用の行列積と Frobenius 分解は読み込みません。通常版の Matrix<T> と同じソースで併用できます。
#include "lib/math/matrix_light.hpp"
#include "lib/fps/fps_light.hpp"
using M = LightMatrix<light_fps>;
M a{ {light_fps{1}, light_fps{0, 1}},
{light_fps{0, 1}, light_fps{1}} };
int degree = 8;
auto b = a.pow(1000000000000000000ULL, [degree](M &m) {
for (int i = 0; i < m.rows(); ++i)
for (auto &f : m[i])
if (f.size() > std::size_t(degree)) f.resize(degree);
});
この例は \(x^8\) を法として FPS を要素とする行列を累乗します。pow(k,normalize) は単位行列・入力行列と、各積の結果に normalize を呼びます。FPS の次数は積のたびに打ち切るため、大きな指数も扱えます。通常版でも Matrix<fps> と同じ書き方ができます。
型 T に関する条件や、累乗以外の API は 通常版 を参照してください。固定法 modint の行列で速度を重視する場合は通常版を使います。通常版の行列単独も 64 KiB 以内に収まります。
行列積 \(n\times k\) と \(k\times m\) は \(O(nkm)\)、正方行列の累乗は \(O(n^3\log k)\)、追加空間は \(O(nk+km+nm)\) です。要素の演算・コピーを \(O(1)\) とします。FPS を要素にするときは、その演算の計算量を掛けます。normalize の処理時間は別途必要です。
依存ライブラリ
コード
#ifndef MATRIX_LIGHT_HPP
#define MATRIX_LIGHT_HPP
#include "../components/matrix_core.hpp"
namespace matrix_detail {
template <class T>
struct PlainMatrixArithmetic {
static constexpr bool modular = false;
template <class F>
static auto field(std::vector<T> &a, F &&f) {
static_assert(!integer<T>);
return f.template operator()<Field<T>>(a);
}
static std::vector<T> multiply(const std::vector<T> &a, const std::vector<T> &b, int n, int inner, int m) {
return basic_multiply(a, b, n, inner, m);
}
static std::optional<std::vector<T>> power(const std::vector<T> &, int, std::uint64_t) { return std::nullopt; }
};
}
template <class T>
using LightMatrix = Matrix<T, matrix_detail::PlainMatrixArithmetic<T>>;
#endif
#line 1 "math/matrix_light.hpp"
#line 1 "components/matrix_core.hpp"
#line 1 "components/matrix_base_arithmetic.hpp"
#line 1 "components/field_base.hpp"
#include <cmath>
#include <type_traits>
namespace field_detail {
template <class T>
struct Field {
using V = T;
static T add(T a, T b) { return a + b; }
static T mul(T a, T b) { return a * b; }
static T neg(T a) { return -a; }
static T inv(T a) { return T(1) / a; }
static bool zero(T a, long double eps) {
if constexpr (std::is_floating_point_v<T>) return std::abs(a) <= eps;
else { (void)eps; return a == T(0); }
}
static bool better(T a, T b) {
if constexpr (std::is_floating_point_v<T>) return std::abs(a) > std::abs(b);
else { (void)a; (void)b; return false; }
}
static T dot(const T *a, const T *b, int n) {
T result = 0;
for (int i = 0; i < n; ++i) result += a[i] * b[i];
return result;
}
static void scale(T *a, int n, T x) { for (int i = 0; i < n; ++i) a[i] *= x; }
static void axpy(T *a, const T *b, int n, T x) { for (int i = 0; i < n; ++i) a[i] += x * b[i]; }
};
}
#line 5 "components/matrix_base_arithmetic.hpp"
#include <algorithm>
#include <cassert>
#include <cstddef>
#include <cstdint>
#include <limits>
#line 11 "components/matrix_base_arithmetic.hpp"
#include <utility>
#include <vector>
namespace matrix_detail {
template <class T>
inline constexpr bool integer = std::is_integral_v<T> || std::is_same_v<T, __int128_t> || std::is_same_v<T, __uint128_t>;
using field_detail::Field;
template <class Ops>
struct Elimination {
std::vector<int> pivots;
typename Ops::V determinant = 1;
};
template <class Ops, bool exact = false>
Elimination<Ops> eliminate(std::vector<typename Ops::V> &a, int h, int stride, int columns, bool reduced, long double eps) {
using V = typename Ops::V;
Elimination<Ops> result;
result.pivots.reserve(std::min(h, columns));
int rank = 0;
for (int col = 0; col < columns && rank < h; ++col) {
int p = -1;
for (int i = rank; i < h; ++i) {
V x = a[std::size_t(i) * stride + col];
if (!Ops::zero(x, eps) && (p < 0 || Ops::better(x, a[std::size_t(p) * stride + col]))) p = i;
}
if (p < 0) continue;
V *row = a.data() + std::size_t(rank) * stride;
if (p != rank) {
std::swap_ranges(row, row + stride, a.data() + std::size_t(p) * stride);
result.determinant = Ops::neg(result.determinant);
}
result.determinant = Ops::mul(result.determinant, row[col]);
Ops::scale(row + col + 1, stride - col - 1, Ops::inv(row[col]));
row[col] = V(1);
for (int i = rank + 1; i < h; ++i) {
V *dst = a.data() + std::size_t(i) * stride;
V x = dst[col];
if (!Ops::zero(x, eps)) Ops::axpy(dst + col + 1, row + col + 1, stride - col - 1, Ops::neg(x));
dst[col] = V(0);
}
result.pivots.push_back(col);
++rank;
}
bool full = exact && rank == columns;
if (reduced) for (int k = rank - 1; k >= 0; --k) {
int col = result.pivots[k], start = full ? columns : col + 1;
const V *row = a.data() + std::size_t(k) * stride;
for (int i = 0; i < k; ++i) {
V *dst = a.data() + std::size_t(i) * stride;
V x = dst[col];
if (!Ops::zero(x, eps) && (!full || start < stride)) Ops::axpy(dst + start, row + start, stride - start, Ops::neg(x));
dst[col] = V(0);
}
}
return result;
}
template <class T>
[[gnu::noinline]] T integer_determinant(const std::vector<T> &a, int n) {
static_assert(std::numeric_limits<T>::is_signed);
using W = __int128_t;
std::vector<W> work(a.begin(), a.end());
W previous = 1, sign = 1;
if (!n) return T(1);
for (int k = 0; k + 1 < n; ++k) {
int pivot = k;
while (pivot < n && work[std::size_t(pivot) * n + k] == 0) ++pivot;
if (pivot == n) return T(0);
if (pivot != k) { for (int j = 0; j < n; ++j) std::swap(work[std::size_t(k) * n + j], work[std::size_t(pivot) * n + j]); sign = -sign; }
W p = work[std::size_t(k) * n + k];
for (int i = k + 1; i < n; ++i) {
for (int j = k + 1; j < n; ++j) {
W value = work[std::size_t(i) * n + j] * p - work[std::size_t(i) * n + k] * work[std::size_t(k) * n + j];
assert(value % previous == 0);
work[std::size_t(i) * n + j] = value / previous;
}
work[std::size_t(i) * n + k] = 0;
}
previous = p;
}
W result = sign * work.back();
assert(result >= std::numeric_limits<T>::lowest() && result <= std::numeric_limits<T>::max());
return T(result);
}
template <class T>
std::vector<T> basic_multiply(const std::vector<T> &a, const std::vector<T> &b, int n, int inner, int m) {
std::vector<T> c(std::size_t(n) * m, T(0));
for (int i = 0; i < n; ++i) for (int k = 0; k < inner; ++k) {
const T &x = a[std::size_t(i) * inner + k];
for (int j = 0; j < m; ++j) c[std::size_t(i) * m + j] += x * b[std::size_t(k) * m + j];
}
return c;
}
}
#line 9 "components/matrix_core.hpp"
#include <initializer_list>
#line 11 "components/matrix_core.hpp"
#include <optional>
#include <span>
#line 16 "components/matrix_core.hpp"
template <class T>
struct LinearSolution {
std::vector<T> particular;
std::vector<std::vector<T>> basis;
};
namespace matrix_detail {
template <class T> struct FastMatrixArithmetic;
}
template <class T, class Arithmetic = matrix_detail::FastMatrixArithmetic<T>>
class [[nodiscard, gnu::warn_unused]] Matrix {
static_assert(!std::is_same_v<T, bool>);
int h = 0, w = 0;
std::vector<T> a;
std::pair<T, std::vector<int>> reduce(int columns, bool reduced, long double eps) {
assert(eps >= 0);
return Arithmetic::field(a, [&]<class Ops>(auto &work) {
auto res = [&] {
if constexpr (requires { Arithmetic::template eliminate<Ops>(work, h, w, columns, reduced, eps); })
return Arithmetic::template eliminate<Ops>(work, h, w, columns, reduced, eps);
else return matrix_detail::eliminate<Ops, Arithmetic::modular>(work, h, w, columns, reduced, eps);
}();
return std::pair{T(res.determinant), std::move(res.pivots)};
});
}
static std::vector<std::vector<T>> nullspace(const Matrix &r, const std::vector<int> &piv, int vars) {
std::vector<bool> pivot(vars);
for (int p : piv) pivot[p] = true;
std::vector<std::vector<T>> basis;
basis.reserve(vars - piv.size());
for (int col = 0; col < vars; ++col) if (!pivot[col]) {
std::vector<T> v(vars, T(0));
v[col] = T{1};
for (int i = 0; i < int(piv.size()); ++i) v[piv[i]] = -r(i, col);
basis.push_back(std::move(v));
}
return basis;
}
Matrix augment() const {
static_assert(!matrix_detail::integer<T>);
assert(h == w && w <= std::numeric_limits<int>::max() / 2);
Matrix work(h, w * 2);
for (int i = 0; i < h; ++i) {
std::copy((*this)[i].begin(), (*this)[i].end(), work[i].begin());
work(i, w + i) = T{1};
}
return work;
}
std::pair<Matrix, std::vector<int>> echelon(bool reduced, long double eps) const {
Matrix work = *this;
auto piv = work.reduce(w, reduced, eps).second;
return {std::move(work), std::move(piv)};
}
public:
using value_type = T;
static constexpr long double default_epsilon = std::is_floating_point_v<T> ? 1e-12L : 0;
Matrix() = default;
explicit Matrix(int n) : Matrix(n, n) {}
Matrix(int rows, int columns, T value = T(0)) : h(rows), w(columns) {
assert(rows >= 0 && columns >= 0);
a.assign(std::size_t(h) * w, value);
}
Matrix(std::initializer_list<std::initializer_list<T>> values) : Matrix(int(values.size()), values.size() ? int(values.begin()->size()) : 0) {
int i = 0;
for (const auto &row : values) {
assert(int(row.size()) == w);
std::copy(row.begin(), row.end(), (*this)[i++].begin());
}
}
[[nodiscard]] explicit Matrix(const std::vector<std::vector<T>> &values) : Matrix(int(values.size()), values.empty() ? 0 : int(values[0].size())) {
for (int i = 0; i < h; ++i) {
assert(int(values[i].size()) == w);
std::copy(values[i].begin(), values[i].end(), (*this)[i].begin());
}
}
Matrix(const Matrix &) = default;
Matrix &operator=(const Matrix &) = default;
Matrix(Matrix &&other) noexcept { swap(other); }
Matrix &operator=(Matrix &&other) noexcept {
if (this != &other) { Matrix next(std::move(other)); swap(next); }
return *this;
}
void swap(Matrix &other) noexcept { std::swap(h, other.h); std::swap(w, other.w); a.swap(other.a); }
[[nodiscard]] int rows() const { return h; }
[[nodiscard]] int cols() const { return w; }
[[nodiscard]] bool empty() const { return !h || !w; }
[[nodiscard]] T *data() { return a.data(); }
[[nodiscard]] const T *data() const { return a.data(); }
[[nodiscard]] std::span<T> operator[](int i) { assert(0 <= i && i < h); return {w ? a.data() + std::size_t(i) * w : a.data(), std::size_t(w)}; }
[[nodiscard]] std::span<const T> operator[](int i) const { assert(0 <= i && i < h); return {w ? a.data() + std::size_t(i) * w : a.data(), std::size_t(w)}; }
[[nodiscard]] T &operator()(int i, int j) { assert(0 <= i && i < h && 0 <= j && j < w); return a[std::size_t(i) * w + j]; }
[[nodiscard]] const T &operator()(int i, int j) const { assert(0 <= i && i < h && 0 <= j && j < w); return a[std::size_t(i) * w + j]; }
[[nodiscard]] bool operator==(const Matrix &) const = default;
static Matrix identity(int n) {
Matrix res(n);
for (int i = 0; i < n; ++i) res(i, i) = T{1};
return res;
}
void fill(T value) { std::fill(a.begin(), a.end(), value); }
Matrix transpose() const {
Matrix res(w, h);
for (int i = 0; i < h; i += 32) for (int j = 0; j < w; j += 32)
for (int x = i; x < std::min(h, i + 32); ++x) for (int y = j; y < std::min(w, j + 32); ++y) res(y, x) = (*this)(x, y);
return res;
}
Matrix submatrix(int row_l, int row_r, int col_l, int col_r) const {
assert(0 <= row_l && row_l <= row_r && row_r <= h && 0 <= col_l && col_l <= col_r && col_r <= w);
Matrix res(row_r - row_l, col_r - col_l);
for (int i = row_l; i < row_r; ++i) std::copy((*this)[i].begin() + col_l, (*this)[i].begin() + col_r, res[i - row_l].begin());
return res;
}
Matrix operator-() const { Matrix res = *this; for (auto &x : res.a) x = -x; return res; }
Matrix &operator+=(const Matrix &b) { assert(h == b.h && w == b.w); for (std::size_t i = 0; i < a.size(); ++i) a[i] += b.a[i]; return *this; }
Matrix &operator-=(const Matrix &b) { assert(h == b.h && w == b.w); for (std::size_t i = 0; i < a.size(); ++i) a[i] -= b.a[i]; return *this; }
Matrix &operator*=(const Matrix &b) { return *this = *this * b; }
Matrix &operator*=(T x) { for (auto &v : a) v *= x; return *this; }
Matrix &operator/=(T x) { static_assert(!matrix_detail::integer<T>); assert(x != T(0)); return *this *= T{1} / x; }
friend Matrix operator+(Matrix a, const Matrix &b) { a += b; return a; }
friend Matrix operator-(Matrix a, const Matrix &b) { a -= b; return a; }
friend Matrix operator*(const Matrix &a, const Matrix &b) {
assert(a.w == b.h);
Matrix res;
res.h = a.h; res.w = b.w;
res.a = Arithmetic::multiply(a.a, b.a, a.h, a.w, b.w);
return res;
}
friend Matrix operator*(Matrix a, T x) { a *= x; return a; }
friend Matrix operator*(T x, Matrix a) { a *= x; return a; }
friend Matrix operator/(Matrix a, T x) { a /= x; return a; }
[[nodiscard]] std::vector<T> operator*(const std::vector<T> &v) const {
assert(int(v.size()) == w);
return Arithmetic::multiply(a, v, h, w, 1);
}
Matrix pow(std::uint64_t exponent) const {
assert(h == w);
if (auto power = Arithmetic::power(a, h, exponent)) {
Matrix res(h);
res.a = std::move(*power);
return res;
}
return pow(exponent, [](Matrix &) {});
}
template <class Normalize>
Matrix pow(std::uint64_t exponent, Normalize normalize) const {
assert(h == w);
Matrix res = identity(h), x = *this;
normalize(res);
normalize(x);
while (exponent) {
if (exponent & 1) {
res *= x;
normalize(res);
}
exponent >>= 1;
if (exponent) {
x *= x;
normalize(x);
}
}
return res;
}
[[nodiscard]] T trace() const { assert(h == w); T res = T(0); for (int i = 0; i < h; ++i) res += (*this)(i, i); return res; }
[[nodiscard]] T det(long double eps = default_epsilon) const {
assert(h == w);
if constexpr (matrix_detail::integer<T>) {
return matrix_detail::integer_determinant(a, h);
} else {
Matrix work = *this;
auto [value, piv] = work.reduce(w, false, eps);
return int(piv.size()) == h ? value : T(0);
}
}
[[nodiscard]] int rank(long double eps = default_epsilon) const { Matrix work = *this; return int(work.reduce(w, false, eps).second.size()); }
[[nodiscard]] std::pair<Matrix, std::vector<int>> row_echelon(long double eps = default_epsilon) const { return echelon(false, eps); }
[[nodiscard]] std::pair<Matrix, std::vector<int>> rref(long double eps = default_epsilon) const { return echelon(true, eps); }
[[nodiscard]] std::optional<Matrix> inverse(long double eps = default_epsilon) const {
assert(h == w && eps >= 0);
if constexpr (requires { Arithmetic::inverse(a, h, false); }) if (auto fast = Arithmetic::inverse(a, h, false)) {
if (fast->second.empty()) return std::nullopt;
Matrix res(h);
res.a = std::move(fast->second);
return res;
}
auto work = augment();
auto piv = work.reduce(w, true, eps).second;
if (int(piv.size()) < h) return std::nullopt;
return work.submatrix(0, h, w, w * 2);
}
[[nodiscard]] std::vector<std::vector<T>> kernel(long double eps = default_epsilon) const {
auto [work, piv] = rref(eps);
return nullspace(work, piv, w);
}
[[nodiscard]] std::optional<LinearSolution<T>> solve(const std::vector<T> &b, long double eps = default_epsilon) const {
assert(int(b.size()) == h && w < std::numeric_limits<int>::max());
Matrix work(h, w + 1);
for (int i = 0; i < h; ++i) {
std::copy((*this)[i].begin(), (*this)[i].end(), work[i].begin());
work(i, w) = b[i];
}
auto piv = work.reduce(w, true, eps).second;
for (int i = int(piv.size()); i < h; ++i) if (!matrix_detail::Field<T>::zero(work(i, w), eps)) return std::nullopt;
LinearSolution<T> res{std::vector<T>(w, T(0)), nullspace(work, piv, w)};
for (int i = 0; i < int(piv.size()); ++i) res.particular[piv[i]] = work(i, w);
return res;
}
[[nodiscard]] std::vector<T> characteristic_polynomial(long double eps = default_epsilon) const {
assert(h == w && eps >= 0);
if constexpr (requires { Arithmetic::characteristic(a, h); })
if (auto res = Arithmetic::characteristic(a, h)) return std::move(*res);
auto work = a;
return Arithmetic::field(work, [&]<class Ops>(auto &b) {
using V = typename Ops::V;
auto at = [&](int i, int j) -> V& { return b[std::size_t(i) * w + j]; };
std::vector<V> factors;
for (int k = 1; k < h; ++k) {
int pivot = -1;
for (int i = k; i < h; ++i) if (!Ops::zero(at(i, k - 1), eps)
&& (pivot < 0 || Ops::better(at(i, k - 1), at(pivot, k - 1)))) pivot = i;
if (pivot < 0) continue;
if (pivot != k) {
for (int j = 0; j < w; ++j) std::swap(at(k, j), at(pivot, j));
for (int i = 0; i < h; ++i) std::swap(at(i, k), at(i, pivot));
}
V inverse = Ops::inv(at(k, k - 1));
if constexpr (requires { Arithmetic::template hessenberg<Ops>(b, h, k, inverse, factors); })
if (Arithmetic::template hessenberg<Ops>(b, h, k, inverse, factors)) continue;
for (int i = k + 1; i < h; ++i) {
V factor = Ops::mul(at(i, k - 1), inverse);
if (Ops::zero(factor, eps)) continue;
Ops::axpy(b.data() + std::size_t(i) * w + k, b.data() + std::size_t(k) * w + k, w - k, Ops::neg(factor));
at(i, k - 1) = V(0);
for (int j = 0; j < h; ++j) at(j, k) = Ops::add(at(j, k), Ops::mul(factor, at(j, i)));
}
}
std::vector<std::vector<V>> p(h + 1);
p[0] = {V(1)};
for (int n = 1; n <= h; ++n) {
p[n].assign(n + 1, V(0));
std::copy(p[n - 1].begin(), p[n - 1].end(), p[n].begin() + 1);
Ops::axpy(p[n].data(), p[n - 1].data(), n, Ops::neg(at(n - 1, n - 1)));
V product = 1;
for (int i = n - 2; i >= 0; --i) {
product = Ops::mul(product, at(i + 1, i));
V factor = Ops::mul(product, at(i, n - 1));
Ops::axpy(p[n].data(), p[i].data(), i + 1, Ops::neg(factor));
}
}
std::vector<T> res(h + 1);
for (int i = 0; i <= h; ++i) res[i] = T(p[h][i]);
return res;
});
}
Matrix adjugate(long double eps = default_epsilon) const {
assert(h == w && eps >= 0);
if constexpr (requires { Arithmetic::inverse(a, h, true); }) if (auto fast = Arithmetic::inverse(a, h, true)) {
Matrix res(h);
res.a = std::move(fast->second);
return res;
}
auto work = augment();
auto [det, piv] = work.reduce(w, true, eps);
if (int(piv.size()) == h) return work.submatrix(0, h, w, w * 2) * det;
Matrix res(h);
if (int(piv.size()) + 1 < h) return res;
auto right = nullspace(work, piv, w)[0];
int col = 0;
while (std::find(piv.begin(), piv.end(), col) != piv.end()) ++col;
if ((h - 1 - col) & 1) det = -det;
for (int i = 0; i < h; ++i) for (int j = 0; j < w; ++j) res(i, j) = det * right[i] * work(h - 1, w + j);
return res;
}
};
#line 5 "math/matrix_light.hpp"
namespace matrix_detail {
template <class T>
struct PlainMatrixArithmetic {
static constexpr bool modular = false;
template <class F>
static auto field(std::vector<T> &a, F &&f) {
static_assert(!integer<T>);
return f.template operator()<Field<T>>(a);
}
static std::vector<T> multiply(const std::vector<T> &a, const std::vector<T> &b, int n, int inner, int m) {
return basic_multiply(a, b, n, inner, m);
}
static std::optional<std::vector<T>> power(const std::vector<T> &, int, std::uint64_t) { return std::nullopt; }
};
}
template <class T>
using LightMatrix = Matrix<T, matrix_detail::PlainMatrixArithmetic<T>>;