行列（高速版）

math/matrix.hpp · GitHub で見る

全体サイト（256 KiB向け） · 高速実装

軽量実装（64 KiB）はこちら

#include "lib/math/matrix.hpp"
#include "lib/math/mint.hpp"
#include <cassert>

int main() {
	Matrix<mint> a = { {1, 2}, {3, 4} };
	auto det = a.det();
	auto inverse = a.inverse();
	auto power = a.pow(10);
	auto solution = a.solve({5, 11});
	auto polynomial = a.characteristic_polynomial();
	assert(det == mint(-2));
	assert(inverse && a * *inverse == Matrix<mint>::identity(2));
	assert(power.rows() == 2 && power.cols() == 2);
	assert((solution && solution->particular == std::vector<mint>{1, 2}));
	assert((polynomial == std::vector<mint>{-2, -5, 1}));
}

Matrix<T> は行優先の連続領域に保存する可変サイズの行列です。添字は 0 始まり、区間は半開区間です。bool は使えません。

操作	内容・返り値	計算量
Matrix<T>() / Matrix<T>(n)	0 行 0 列／0 初期化した正方行列	空: \(O(1)\)、正方: \(O(n^2)\)
Matrix<T>(n,m,value = T(0))	n 行 m 列を同じ値で初期化	\(O(nm)\)
Matrix<T>{ {...},{...} } / Matrix<T>(vector<vector<T>>)	係数列から構築。各行の長さは同じにする	\(O(nm)\)
rows() / cols() / empty()	行数／列数／行数または列数が 0 か	\(O(1)\)
a[i][j] / a(i,j)	要素の参照。const 版も利用可	\(O(1)\)
a[i] / data()	行の span／行優先の先頭ポインタ	\(O(1)\)
Matrix<T>::identity(n)	単位行列	\(O(n^2)\)
fill(x)	全要素を同じ値で上書き	\(O(nm)\)
swap(b)	行列を交換	\(O(1)\)
transpose()	転置した行列を返す	\(O(nm)\)
submatrix(row_l,row_r,col_l,col_r)	指定した行・列区間をコピー	切り出す要素数に比例
-a / a + b / a - b	符号反転／同じ形の行列の加減算	\(O(nm)\)
a * b	行列積。a.cols() == b.rows()	\(O(nkm)\)
a * x / x * a / a / x	スカラー倍／スカラー除算。除算は体の要素で非零のみ	\(O(nm)\)
a += b / a -= b / a *= b	行列演算の結果を代入	加減算: \(O(nm)\)、積: \(O(nkm)\)
a *= x / a /= x	スカラー演算の結果を代入	\(O(nm)\)
a == b / a != b	形と全要素を比較	\(O(nm)\)
a * vector<T>	行列と列ベクトルの積。ベクトル長は列数	\(O(nm)\)
pow(k)	正方行列の非負整数乗。k は uint64_t、0 乗は単位行列	固定の素数 modint で分解成功時: \(O(n^3+n^2\log k)\)、その他・最悪: \(O(n^3\log k)\)
pow(k,normalize)	単位行列・入力行列・各積の結果に normalize(Matrix&) を適用して累乗	\(O(n^3\log k)\) ＋ normalize の処理時間
trace()	正方行列の対角成分の和	\(O(n)\)
det(eps = default_epsilon)	行列式。特異なら 0、0 行 0 列なら 1	\(O(n^3)\)
inverse(eps = default_epsilon)	逆行列を optional<Matrix<T>> で返す。特異なら nullopt	\(O(n^3)\)
rank(eps = default_epsilon)	行列の rank	\(O(nm\min(n,m))\)
row_echelon(eps = default_epsilon)	行階段形とピボットの列番号の組。各ピボットは 1	\(O(nm\min(n,m))\)
rref(eps = default_epsilon)	簡約行階段形とピボットの列番号の組	\(O(nm\min(n,m))\)
solve(b,eps = default_epsilon)	\(Ax=b\) の解。解なしなら nullopt	\(O(nm\min(n,m)+m^2)\)
kernel(eps = default_epsilon)	\(Ax=0\) の解空間の基底。各要素が列ベクトル	\(O(nm\min(n,m)+m^2)\)
characteristic_polynomial(eps = default_epsilon)	\(\det(xI-A)\) の係数を定数項から返す	\(O(n^3)\)
adjugate(eps = default_epsilon)	余因子行列の転置。特異行列にも対応	\(O(n^3)\)

solve は optional<LinearSolution<T>> を返します。b の長さは行数とします。

フィールド	内容
particular	解を 1 つ返す。自由変数は 0
basis	同次方程式の解空間の基底。自由変数の列番号順

全ての解は particular と basis の線形結合の和です。ピボットの列番号も昇順です。行階段形などを返す操作を含め、元の行列は変更しません。

加減算・乗算・累乗には環の要素を使えます。除算・逆行列・rank・掃き出し・方程式・核・特性多項式・余因子行列には、非零要素で割れる体を使います。素数 mod の modint や Fraction<> に対応します。通常の整数型をこれらの除算を伴う操作には使えません。

符号付き整数型の det は例外として Bareiss 法で厳密に計算します。途中の積・差・商は __int128_t、結果は T に収まる必要があります。整数の他の演算も途中の値が T に収まる範囲で使います。

浮動小数点では絶対値が最大のピボットを選び、絶対値が eps 以下なら 0 と判定します。既定は 1e-12L、それ以外の型は 0 です。eps >= 0 とし、必要に応じて各操作の引数で変更します。浮動小数点の結果には丸め誤差を含みます。

大きな累乗では、上記の modint の法が素数なら Frobenius 分解を使います。基底の独立性と各ブロックの関係式を確認し、確認できない場合は繰り返し二乗法に切り替えます。合成数 mod・整数・浮動小数点でも pow を使えます。

FPS も要素型にでき、Matrix<fps> の単位行列の対角成分は定数 FPS の {1} になります。\(x^D\) を法として計算するときは、pow(k,normalize) のコールバックで各要素を先頭 D 項に打ち切ってください。この累乗は繰り返し二乗法を使います。例は 軽量版の行列 にあります。FPS と行列を組み合わせて 64 KiB 制限で使う場合は LightMatrix<light_fps> を使います。

\(n\times m\) 行列について、構築・加減算・スカラー演算・転置は \(O(nm)\)、行列積 \(n\times k\) と \(k\times m\) は \(O(nkm)\)、行列とベクトルの積は \(O(nm)\)、trace は \(O(n)\) です。正方行列の累乗は固定の素数 modint で分解に成功すれば \(O(n^3+n^2\log k)\)、最悪や他の型で \(O(n^3\log k)\)、行列式・逆行列・特性多項式・余因子行列は \(O(n^3)\) です。長方行列の rank・掃き出し・方程式・核は、基底の出力も含めて \(O(nm\min(n,m)+m^2)\) です。rank・掃き出しだけなら \(O(nm\min(n,m))\) です。追加空間は行列式・rank・掃き出しが \(O(nm)\)、方程式・核は \(O(nm+m^2)\)、正方行列の各操作は \(O(n^2)\) です。行列積は入力の変換と出力を含めて \(O(nk+km+nm)\) です。要素の演算・コピーを \(O(1)\) とします。

依存ライブラリ

• 体の内部演算
• 体の基本演算
• 行列の内部演算
• 行列の基本演算
• 行列のブロック逆行列・余因子行列
• 行列の共通クラス
• 固定法の行列消去
• 行列累乗・特性多項式の内部演算

検証コード

• test/matrix_adjugate.test.cpp
• test/matrix_characteristic_polynomial.test.cpp
• test/matrix_det.test.cpp
• test/matrix_inverse.test.cpp
• test/matrix_pow.test.cpp
• test/matrix_product.test.cpp
• test/matrix_rank.test.cpp
• test/matrix_solve.test.cpp

コード

#ifndef MATRIX_HPP
#define MATRIX_HPP

#include "../components/matrix_arithmetic.hpp"
#include "../components/matrix_power.hpp"
#include "../components/matrix_elimination_modular.hpp"

namespace matrix_detail {
template <class T>
std::optional<std::pair<T, std::vector<T>>> block_inverse(const std::vector<T> &a, int n, bool adjugate);

template <class T>
struct FastMatrixArithmetic {
	static constexpr bool modular = small_modint<T>;
	template <class Ops>
	static auto eliminate(std::vector<typename Ops::V> &a, int h, int stride, int columns, bool reduced, long double eps) {
		if constexpr (modular) if (h >= 32 && columns >= 32) return ModularElimination<Ops>::run(a, h, stride, columns, reduced);
		return matrix_detail::eliminate<Ops, modular>(a, h, stride, columns, reduced, eps);
	}
	template <class Ops>
	static bool hessenberg(std::vector<typename Ops::V> &a, int n, int k, typename Ops::V inverse, std::vector<typename Ops::V> &factors) {
		if constexpr (modular) if (n >= 16) {
			factors.resize(n);
			ModularElimination<Ops>::hessenberg(a, n, k, inverse, factors);
			return true;
		}
		return false;
	}

	template <class F>
	static auto field(std::vector<T> &a, F &&f) { return matrix_detail::field(a, std::forward<F>(f)); }
	static std::vector<T> multiply(const std::vector<T> &a, const std::vector<T> &b, int n, int inner, int m) {
		return matrix_detail::multiply(a, b, n, inner, m);
	}
	static std::optional<std::pair<T, std::vector<T>>> inverse(const std::vector<T> &a, int n, bool adjugate) {
		if constexpr (modular) return block_inverse(a, n, adjugate);
		return std::nullopt;
	}
	static std::optional<std::vector<T>> characteristic(const std::vector<T> &a, int n) {
		if constexpr (small_modint<T>) if constexpr (prime_modulus(T::mod())) {
			if (n >= 32) return FrobeniusPower<T>(a, n).characteristic();
		}
		return std::nullopt;
	}
	static std::optional<std::vector<T>> power(const std::vector<T> &a, int n, std::uint64_t exponent) {
		if constexpr (small_modint<T>) if constexpr (prime_modulus(T::mod())) {
			if (n >= 16 && exponent >= 256) return FrobeniusPower<T>(a, n).run(exponent);
		}
		return std::nullopt;
	}
};
}

#include "../components/matrix_core.hpp"
#include "../components/matrix_block_inverse.hpp"

#endif

Back to top page