Skip to content

Commit 94056a4

Browse files
hitonanodehitonanode
committed
Determinant of Matrix (arbitrary mod)
1 parent 21595d8 commit 94056a4

6 files changed

+99
-7
lines changed

‎linear_algebra_matrix/characteristic_poly.hpp

Lines changed: 3 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1,12 +1,11 @@
11
#pragma once
2-
#include "hessenberg_reduction.hpp"
2+
#include "matrix.hpp"
33
#include <vector>
44

5-
// Characteristic polynomial of matrix M (|xI - M|)
5+
// Characteristic polynomial of upper Hessenberg matrix M (|xI - M|)
66
// Complexity: O(n^3)
77
// R. Rehman, I. C. Ipsen, "La Budde's Method for Computing Characteristic Polynomials," 2011.
8-
template <class Tp> std::vector<Tp> characteristic_poly(matrix<Tp> M) {
9-
hessenberg_reduction(M);
8+
template <class Tp> std::vector<Tp> characteristic_poly_of_hessenberg(matrix<Tp> &M) {
109
const int N = M.height();
1110
// p[i + 1] = (Characteristic polynomial of i-th leading principal minor)
1211
std::vector<std::vector<Tp>> p(N + 1);

‎linear_algebra_matrix/characteristic_poly.md

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -3,11 +3,11 @@ title: Characteristic polynomial （行列の特性多項式）
33
documentation_of: ./characteristic_poly.hpp
44
---
55

6-
体上の $n$ 次正方行列 $M$ の特性多項式 $p_M(x) = \det (xI - M)$ を $O(n^3)$ で求める決定的アルゴリズム．
6+
Upper Hessenberg な体上の $n$ 次正方行列 $M$ の特性多項式 $p_M(x) = \det (xI - M)$ を $O(n^3)$ で求める決定的アルゴリズム．
77

88
## やっていること
99

10-
入力として与えられた正方行列 $M$ に (upper) Hessenberg reduction を施した後，これに対して $M$ の首座小行列の特性多項式を再帰的に求めていく（La Budde's method, [1]）．
10+
入力として与えられた (upper) Hessenberg な正方行列 $M$ に対して $M$ の首座小行列の特性多項式を再帰的に求めていく（La Budde's method, [1]）．
1111

1212
## References
1313

‎linear_algebra_matrix/hessenberg_reduction.hpp

Lines changed: 41 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -24,3 +24,44 @@ template <class Tp> void hessenberg_reduction(matrix<Tp> &M) {
2424
}
2525
}
2626
}
27+
28+
template <class Ring> void ring_hessenberg_reduction(matrix<Ring> &M) {
29+
assert(M.height() == M.width());
30+
const int N = M.height();
31+
for (int r = 0; r < N - 2; r++) {
32+
int piv = matrix<Ring>::choose_pivot(M, r + 1, r);
33+
if (piv < 0) continue;
34+
for (int i = 0; i < N; i++) std::swap(M[r + 1][i], M[piv][i]);
35+
for (int i = 0; i < N; i++) std::swap(M[i][r + 1], M[i][piv]);
36+
37+
for (int i = r + 2; i < N; i++) {
38+
if (M[i][r] == Ring()) continue;
39+
Ring a = M[r + 1][r], b = M[i][r], m00 = 1, m01 = 0, m10 = 0, m11 = 1;
40+
while (a != Ring() and b != Ring()) {
41+
if (a.val() > b.val()) {
42+
auto d = a.val() / b.val();
43+
a -= b * d, m00 -= m10 * d, m01 -= m11 * d;
44+
} else {
45+
auto d = b.val() / a.val();
46+
b -= a * d, m10 -= m00 * d, m11 -= m01 * d;
47+
}
48+
}
49+
if (a == Ring()) std::swap(a, b), std::swap(m00, m10), std::swap(m01, m11);
50+
51+
for (int j = 0; j < N; j++) {
52+
Ring anew = M[r + 1][j] * m00 + M[i][j] * m01;
53+
Ring bnew = M[r + 1][j] * m10 + M[i][j] * m11;
54+
M[r + 1][j] = anew;
55+
M[i][j] = bnew;
56+
}
57+
assert(M[i][r] == 0);
58+
59+
for (int j = 0; j < N; j++) {
60+
Ring anew = M[j][r + 1] * m11 - M[j][i] * m10;
61+
Ring bnew = -M[j][r + 1] * m01 + M[j][i] * m00;
62+
M[j][r + 1] = anew;
63+
M[j][i] = bnew;
64+
}
65+
}
66+
}
67+
}

‎linear_algebra_matrix/hessenberg_reduction.md

Lines changed: 21 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -8,3 +8,24 @@ documentation_of: ./hessenberg_reduction.hpp
88
## やっていること
99

1010
(Upper) Hessenberg reduction とは，行列に相似変換を施すことでその対角成分より2つ以上左下側の成分を全てゼロにするというもので，このような変換は特に Householder 変換の組合せによって可能である．相似変換で特性多項式は不変なため，本ライブラリでは特性多項式の導出などに応用されている．
11+
12+
## 使用方法
13+
14+
`matrix<T>` に対して upper Hessenberg reduction を行う関数は以下のように使用する．
15+
### `T` が逆元がとれるデータ構造の場合
16+
17+
```cpp
18+
matrix<T> mat(N, N);
19+
hessenberg_reduction(mat);
20+
```
21+
22+
### `T` が逆元がとれないがユークリッドの互除法が可能なデータ構造の場合（例：合成数 modint）
23+
24+
```cpp
25+
matrix<T> mat(N, N);
26+
ring_hessenberg_reduction(mat);
27+
```
28+
29+
## 問題例
30+
31+
- [Library Checker: Determinant of Matrix (arbitrary mod)](https://judge.yosupo.jp/problem/matrix_det_arbitrary_mod): 任意 mod のケースでも特性多項式が $O(n^3)$ で計算可能である．

‎linear_algebra_matrix/test/characteristic_poly.test.cpp

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,7 @@
11
#define PROBLEM "https://judge.yosupo.jp/problem/characteristic_polynomial"
22
#include "../characteristic_poly.hpp"
33
#include "../../modint.hpp"
4+
#include "../hessenberg_reduction.hpp"
45
#include "../matrix.hpp"
56
#include <algorithm>
67
#include <iostream>
@@ -13,7 +14,8 @@ int main() {
1314
using mint = ModInt<998244353>;
1415
matrix<mint> M(N, N);
1516
cin >> M;
16-
auto poly = characteristic_poly<mint>(M);
17+
hessenberg_reduction(M);
18+
auto poly = characteristic_poly_of_hessenberg<mint>(M);
1719
for (auto x : poly) cout << x << ' ';
1820
cout << '\n';
1921
}

‎linear_algebra_matrix/test/determinant_arbitrary_mod.test.cpp

Lines changed: 29 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,29 @@
1+
#define PROBLEM "https://judge.yosupo.jp/problem/matrix_det_arbitrary_mod"
2+
#if __cplusplus < 201402L
3+
#define IGNORE
4+
#endif
5+
#include "../characteristic_poly.hpp"
6+
#include "../hessenberg_reduction.hpp"
7+
#include "../matrix.hpp"
8+
#include <atcoder/modint>
9+
#include <iostream>
10+
using namespace std;
11+
12+
int main() {
13+
cin.tie(nullptr), ios::sync_with_stdio(false);
14+
int N, M;
15+
cin >> N >> M;
16+
using mint = atcoder::modint;
17+
mint::set_mod(M);
18+
matrix<mint> mat(N, N);
19+
for (int i = 0; i < N; ++i) {
20+
for (int j = 0; j < N; ++j) {
21+
int v;
22+
cin >> v;
23+
mat[i][j] = v;
24+
}
25+
}
26+
mat = -mat;
27+
ring_hessenberg_reduction(mat);
28+
cout << characteristic_poly_of_hessenberg(mat)[0].val() << '\n';
29+
}

0 commit comments

Comments
 (0)