Merge pull request #335 from hitonanode/f2-determinant

F2 determinant

Author: hitonanode

linear_algebra_matrix/linalg_bitset.hpp

@@ -9,7 +9,7 @@
 // Complexity: O(nm + nm rank(M) / 64)
 // Verified: abc276_h (2000 x 8000)
 template <int Wmax>
-std::vector<std::bitset<Wmax>> gauss_jordan(int W, std::vector<std::bitset<Wmax>> M) {
+std::vector<std::bitset<Wmax>> f2_gauss_jordan(int W, std::vector<std::bitset<Wmax>> M) {
     assert(W <= Wmax);
     int H = M.size(), c = 0;
     for (int h = 0; h < H and c < W; ++h, ++c) {
@@ -33,7 +33,7 @@ std::vector<std::bitset<Wmax>> gauss_jordan(int W, std::vector<std::bitset<Wmax>
 }
 
 // Rank of Gauss-Jordan eliminated matrix
-template <int Wmax> int rank_gauss_jordan(int W, const std::vector<std::bitset<Wmax>> &M) {
+template <int Wmax> int f2_rank_gauss_jordan(int W, const std::vector<std::bitset<Wmax>> &M) {
     assert(W <= Wmax);
     for (int h = (int)M.size() - 1; h >= 0; h--) {
         int j = 0;
@@ -43,27 +43,54 @@ template <int Wmax> int rank_gauss_jordan(int W, const std::vector<std::bitset<W
     return 0;
 }
 
+// determinant of F2 matrix.
+// Return 0 if the matrix is singular, otherwise return 1.
+// Complexity: O(W^3 / 64)
+template <int Wmax> int f2_determinant(const std::vector<std::bitset<Wmax>> &M) {
+    const int H = M.size();
+    if (H > Wmax) return 0;
+
+    auto tmp = M;
+    for (int h = 0; h < H; ++h) {
+        int piv = -1;
+        for (int j = h; j < H; ++j) {
+            if (tmp.at(j).test(h)) {
+                piv = j;
+                break;
+            }
+        }
+        if (piv == -1) return 0; // singular
+
+        if (piv != h) std::swap(tmp.at(piv), tmp.at(h));
+        for (int hh = h + 1; hh < H; ++hh) {
+            if (tmp.at(hh).test(h)) tmp.at(hh) ^= tmp.at(h);
+        }
+    }
+
+    return 1; // nonsingular
+}
+
 template <int W1, int W2>
 std::vector<std::bitset<W2>>
-matmul(const std::vector<std::bitset<W1>> &A, const std::vector<std::bitset<W2>> &B) {
+f2_matmul(const std::vector<std::bitset<W1>> &A, const std::vector<std::bitset<W2>> &B) {
     int H = A.size(), K = B.size();
     std::vector<std::bitset<W2>> C(H);
     for (int i = 0; i < H; i++) {
         for (int j = 0; j < K; j++) {
-            if (A[i][j]) C[i] ^= B[j];
+            if (A.at(i).test(j)) C.at(i) ^= B.at(j);
         }
     }
     return C;
 }
 
 template <int Wmax>
-std::vector<std::bitset<Wmax>> matpower(std::vector<std::bitset<Wmax>> X, long long n) {
+std::vector<std::bitset<Wmax>> f2_matpower(std::vector<std::bitset<Wmax>> X, long long n) {
     int D = X.size();
     std::vector<std::bitset<Wmax>> ret(D);
     for (int i = 0; i < D; i++) ret[i][i] = 1;
     while (n) {
-        if (n & 1) ret = matmul<Wmax, Wmax>(ret, X);
-        X = matmul<Wmax, Wmax>(X, X), n >>= 1;
+        if (n & 1) ret = f2_matmul<Wmax, Wmax>(ret, X);
+        X = f2_matmul<Wmax, Wmax>(X, X), n >>= 1;
     }
     return ret;
 }
@@ -74,7 +101,7 @@ std::vector<std::bitset<Wmax>> matpower(std::vector<std::bitset<Wmax>> X, long l
 // Complexity: O(HW + HW rank(A) / 64 + W^2 len(freedoms))
 template <int Wmax, class Vec>
 std::tuple<bool, std::bitset<Wmax>, std::vector<std::bitset<Wmax>>>
-system_of_linear_equations(std::vector<std::bitset<Wmax>> A, Vec b, int W) {
+f2_system_of_linear_equations(std::vector<std::bitset<Wmax>> A, Vec b, int W) {
     int H = A.size();
     assert(W <= Wmax);
     assert(A.size() == b.size());
@@ -84,7 +111,7 @@ system_of_linear_equations(std::vector<std::bitset<Wmax>> A, Vec b, int W) {
         for (int j = 0; j < W; ++j) M[i][j] = A[i][j];
         M[i][W] = b[i];
     }
-    M = gauss_jordan<Wmax + 1>(W + 1, M);
+    M = f2_gauss_jordan<Wmax + 1>(W + 1, M);
     std::vector<int> ss(W, -1);
     std::vector<int> ss_nonneg_js;
     for (int i = 0; i < H; i++) {

linear_algebra_matrix/linalg_bitset.md

@@ -15,11 +15,15 @@ vector<bitset<Wmax>> A;
 vector<bool> b;
 
 // Solve Ax = b (x: F_2^W)
-auto [feasible, x0, freedoms] = system_of_linear_equations<Wmax, vector<bool>>(A, b, W);
+auto [feasible, x0, freedoms] = f2_system_of_linear_equations<Wmax, vector<bool>>(A, b, W);
+
+// Calc determinant (or check whether A is regular)
+int det = f2_determinant<dim>(mat);
 ```
 
 ## 問題例
 
 - [AOJ 2624: Graph Automata Player](https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=2624)
 - [No.1421 国勢調査 (Hard) - yukicoder](https://yukicoder.me/problems/no/1421)
 - [AtCoder Beginner Contest 276 Ex - Construct a Matrix](https://atcoder.jp/contests/abc276/tasks/abc276_h) $2000 \times 8000$ 行列の線型方程式を解く．
+- [Library Checker: Determinant of Matrix (Mod 2)](https://judge.yosupo.jp/problem/matrix_det_mod_2)

linear_algebra_matrix/test/linalg_bitset.test.cpp

@@ -26,9 +26,9 @@ int main() {
     }
 
     cin >> T;
-    A = matpower<Wmax>(A, T);
+    A = f2_matpower<Wmax>(A, T);
     for (int i = 0; i < N; i++) A[i][N] = v[i];
-    A = gauss_jordan<Wmax>(N + 1, A);
+    A = f2_gauss_jordan<Wmax>(N + 1, A);
 
     for (int i = 0; i < N; i++) {
         if (A[i].count() == 1 and A[i][N]) {

linear_algebra_matrix/test/linalg_bitset.yuki1421.test.cpp

@@ -59,7 +59,7 @@ int main() {
             }
             A.emplace_back(a);
         }
-        auto [ok, solution, freedoms] = system_of_linear_equations<Wmax, vector<bool>>(A, b, N);
+        auto [ok, solution, freedoms] = f2_system_of_linear_equations<Wmax, vector<bool>>(A, b, N);
         if (!ok) bad();
         for (int i = 0; i < N; ++i) ret[i] += int(solution[i]) << d;
     }

linear_algebra_matrix/test/linalg_bitset_det.test.cpp

@@ -0,0 +1,28 @@
+#define PROBLEM "https://judge.yosupo.jp/problem/matrix_det_mod_2"
+#include "../linalg_bitset.hpp"
+
+#include <bitset>
+#include <iostream>
+#include <string>
+#include <vector>
+
+using namespace std;
+
+int main() {
+    cin.tie(nullptr);
+    ios::sync_with_stdio(false);
+
+    constexpr int dim = 1 << 12;
+    using BS = bitset<dim>;
+
+    int N;
+    cin >> N;
+    vector<BS> mat(N);
+    for (auto &v : mat) {
+        string s;
+        cin >> s;
+        v = BS(s);
+    }
+
+    cout << f2_determinant<dim>(mat) << '\n';
+}

linear_algebra_matrix/test/linalg_bitset_mul.test.cpp

@@ -0,0 +1,43 @@
+#define PROBLEM "https://judge.yosupo.jp/problem/matrix_product_mod_2"
+#include "../linalg_bitset.hpp"
+
+#include <algorithm>
+#include <bitset>
+#include <iostream>
+#include <string>
+#include <vector>
+using namespace std;
+
+int main() {
+    cin.tie(nullptr);
+    ios::sync_with_stdio(false);
+
+    constexpr int dim = 1 << 12;
+
+    using BS = bitset<dim>;
+
+    int N, M, K;
+    cin >> N >> M >> K;
+
+    vector<BS> A(N), B(M);
+    for (auto &v : A) {
+        string s;
+        cin >> s;
+        std::reverse(s.begin(), s.end()); // bitset に文字列 s を与えると s[0] が最高位になるので反転
+        v = BS(s);
+    }
+
+    for (auto &v : B) {
+        string s;
+        cin >> s;
+        std::reverse(s.begin(), s.end());
+        v = BS(s);
+    }
+
+    auto C = f2_matmul<dim, dim>(A, B);
+    for (const auto &v : C) {
+        auto tmp = v.to_string();
+        std::reverse(tmp.begin(), tmp.end());
+        cout << tmp.substr(0, K) << '\n';
+    }
+}