Tutte matrix matching test

Author: adamant-pwn

cp-algo/linalg/matrix.hpp

@@ -9,6 +9,7 @@
 #include <vector>
 #include <array>
 namespace cp_algo::linalg {
+    enum gauss_mode {normal, reverse};
     template<typename base_t>
     struct matrix: valarray_base<matrix<base_t>, vec<base_t>> {
         using base = base_t;
@@ -127,17 +128,28 @@ namespace cp_algo::linalg {
             }
             return *this;
         }
-
-        enum gauss_mode {normal, reverse};
+        template<gauss_mode mode = normal>
+        void eliminate(size_t i, size_t k) {
+            auto kinv = base(1) / row(i).normalize()[k];
+            for(size_t j = (mode == normal) * i; j < n(); j++) {
+                if(j != i) {
+                    row(j).add_scaled(row(i), -row(j).normalize(k) * kinv);
+                }
+            }
+        }
+        template<gauss_mode mode = normal>
+        void eliminate(size_t i) {
+            row(i).normalize();
+            for(size_t j = (mode == normal) * i; j < n(); j++) {
+                if(j != i) {
+                    row(j).reduce_by(row(i));
+                }
+            }
+        }
         template<gauss_mode mode = normal>
         matrix& gauss() {
             for(size_t i = 0; i < n(); i++) {
-                row(i).normalize();
-                for(size_t j = (mode == normal) * i; j < n(); j++) {
-                    if(j != i) {
-                        row(j).reduce_by(row(i));
-                    }
-                }
+                eliminate<mode>(i);
             }
             return normalize();
         }

verify/linalg/tutte.test.cpp

@@ -0,0 +1,64 @@
+// @brief Matching on General Graph (Tutte matrix)
+#define PROBLEM "https://judge.yosupo.jp/problem/general_matching"
+#pragma GCC optimize("Ofast,unroll-loops")
+#pragma GCC target("tune=native")
+#include "cp-algo/linalg/matrix.hpp"
+#include "cp-algo/random/rng.hpp"
+#include <bits/stdc++.h>
+
+using namespace std;
+using namespace cp_algo::math;
+using namespace cp_algo::linalg;
+using namespace cp_algo::random;
+
+const int mod = 998244353;
+using base = modint<mod>;
+
+void solve() {
+    int n, m;
+    cin >> n >> m;
+    matrix<base> T(n);
+    for(int i = 0; i < m; i++) {
+        int u, v;
+        cin >> u >> v;
+        base x = rng();
+        T[u][v] += x;
+        T[v][u] -= x;
+    }
+    auto [pivots, free] = matrix(T).echelonize();
+    matrix<base> B(size(pivots));
+    for(int i = 0; i < ssize(pivots); i++) {
+        for(int j = 0; j < ssize(pivots); j++) {
+            B[i][j] = T[pivots[i]][pivots[j]];
+        }
+    }
+    B = *B.inv();
+    vector<pair<int, int>> ans;
+    for(size_t i = 0; i < size(pivots); i++) {
+        for(size_t j = 0; j < size(pivots); j++) {
+            if(T[pivots[i]][pivots[j]] != 0 && B[i][j] != 0) {
+                ans.emplace_back(pivots[i], pivots[j]);
+                B.eliminate<gauss_mode::reverse>(i, j);
+                B.eliminate<gauss_mode::reverse>(j, i);
+                B.normalize();
+                B[i] *= 0;
+                B[j] *= 0;
+            }
+        }
+    }
+    cout << size(ans) << "\n";
+    for(auto [u, v]: ans) {
+        cout << u << ' ' << v << "\n";
+    }
+}
+
+signed main() {
+    //freopen("input.txt", "r", stdin);
+    ios::sync_with_stdio(0);
+    cin.tie(0);
+    int t = 1;
+    // cin >> t;
+    while(t--) {
+        solve();
+    }
+}