Factor Frobenius basis out from the test

Author: adamant-pwn

cp-algo/linalg/all

@@ -2,4 +2,5 @@
 #define CP_ALGO_LINALG_ALL
 #include "vector.hpp"
 #include "matrix.hpp"
+#include "frobenius.hpp"
 #endif // CP_ALGO_LINALG_ALL

cp-algo/linalg/frobenius.hpp

@@ -0,0 +1,66 @@
+#ifndef CP_ALGO_LINALG_FROBENIUS_HPP
+#define CP_ALGO_LINALG_FROBENIUS_HPP
+#include "matrix.hpp"
+#include "../algebra/poly.hpp"
+#include <vector>
+namespace cp_algo::linalg {    
+    template<bool reduce = false>
+    auto frobenius_basis(auto const& A) {
+        using matrix = std::decay_t<decltype(A)>;
+        using base = matrix::base;
+        using polyn = algebra::poly_t<base>;
+        assert(A.n() == A.m());
+        size_t n = A.n();
+        struct krylov {
+            std::vector<vec<base>> basis; 
+            polyn rec;
+        };
+        std::vector<krylov> blocks;
+        std::vector<vec<base>> reduced;
+        while(size(reduced) < n) {
+            auto generate_block = [&](auto x) {
+                krylov block;
+                while(true) {
+                    vec<base> y = x | vec<base>::ei(n + 1, size(reduced));
+                    for(auto &it: reduced) {
+                        y.reduce_by(it);
+                    }
+                    y.normalize();
+                    if(vec<base>(y[std::slice(0, n, 1)]) == vec<base>(n)) {
+                        block.rec = std::vector<base>(
+                            begin(y) + n + size(reduced) - size(block.basis),
+                            begin(y) + n + size(reduced) + 1
+                        );
+                        return std::pair{block, vec<base>(y[std::slice(n, n, 1)])};
+                    } else {
+                        block.basis.push_back(x);
+                        reduced.push_back(y);
+                        x = A.apply(x);
+                    }
+                }
+            };
+            auto [block, full_rec] = generate_block(vec<base>::random(n));
+            if constexpr (reduce) {
+                if(vec<base>(full_rec[std::slice(0, size(reduced), 1)]) != vec<base>(size(reduced))) {
+                    auto x = block.basis[0];
+                    size_t start = 0;
+                    for(auto &[basis, rec]: blocks) {
+                        polyn cur_rec = std::vector<base>(
+                            begin(full_rec) + start, begin(full_rec) + start + rec.deg()
+                        );
+                        auto shift = cur_rec / block.rec;
+                        for(int j = 0; j <= shift.deg(); j++) {
+                            x.add_scaled(basis[j], shift[j]);
+                        }
+                        start += rec.deg();
+                    }
+                    reduced.erase(begin(reduced) + start, end(reduced));
+                    tie(block, full_rec) = generate_block(x.normalize());
+                }
+            }
+            blocks.push_back(block);
+        }
+        return blocks;
+    }
+};
+#endif // CP_ALGO_LINALG_FROBENIUS_HPP

cp-algo/linalg/matrix.hpp

@@ -9,8 +9,9 @@
 #include <vector>
 #include <array>
 namespace cp_algo::linalg {
-    template<typename base>
-    struct matrix: valarray_base<matrix<base>, vec<base>> {
+    template<typename base_t>
+    struct matrix: valarray_base<matrix<base_t>, vec<base_t>> {
+        using base = base_t;
         using Base = valarray_base<matrix<base>, vec<base>>;
         using Base::Base;
 

verify/algebra/matrix/characteristic.test.cpp

@@ -2,8 +2,7 @@
 #define PROBLEM "https://judge.yosupo.jp/problem/characteristic_polynomial"
 #pragma GCC optimize("Ofast,unroll-loops")
 #pragma GCC target("avx2,tune=native")
-#include "cp-algo/algebra/poly.hpp"
-#include "cp-algo/linalg/matrix.hpp"
+#include "cp-algo/linalg/frobenius.hpp"
 #include <bits/stdc++.h>
 
 using namespace std;
@@ -14,62 +13,6 @@ const int mod = 998244353;
 using base = modular<mod>;
 using polyn = poly_t<base>;
 
-template<bool reduce = false>
-auto frobenius_basis(auto const& A) {
-    assert(A.n() == A.m());
-    size_t n = A.n();
-    struct krylov {
-        vector<vec<base>> basis; 
-        poly_t<base> rec;
-    };
-    vector<krylov> blocks;
-    vector<vec<base>> reduced;
-    while(size(reduced) < n) {
-        auto generate_block = [&](auto x) {
-            krylov block;
-            while(true) {
-                vec<base> y = x | vec<base>::ei(n + 1, size(reduced));
-                for(auto &it: reduced) {
-                    y.reduce_by(it);
-                }
-                y.normalize();
-                if(vec<base>(y[slice(0, n, 1)]) == vec<base>(n)) {
-                    block.rec = vector<base>(
-                        begin(y) + n + size(reduced) - size(block.basis),
-                        begin(y) + n + size(reduced) + 1
-                    );
-                    return std::pair{block, vec<base>(y[slice(n, n, 1)])};
-                } else {
-                    block.basis.push_back(x);
-                    reduced.push_back(y);
-                    x = A.apply(x);
-                }
-            }
-        };
-        auto [block, full_rec] = generate_block(vec<base>::random(n));
-        if constexpr (reduce) {
-            if(vec<base>(full_rec[slice(0, size(reduced), 1)]) != vec<base>(size(reduced))) {
-                auto x = block.basis[0];
-                size_t start = 0;
-                for(auto &[basis, rec]: blocks) {
-                    polyn cur_rec = vector<base>(
-                        begin(full_rec) + start, begin(full_rec) + start + rec.deg()
-                    );
-                    auto shift = cur_rec / block.rec;
-                    for(int j = 0; j <= shift.deg(); j++) {
-                        x.add_scaled(basis[j], shift[j]);
-                    }
-                    start += rec.deg();
-                }
-                reduced.erase(begin(reduced) + start, end(reduced));
-                tie(block, full_rec) = generate_block(x.normalize());
-            }
-        }
-        blocks.push_back(block);
-    }
-    return blocks;
-}
-
 void solve() {
     size_t n;
     cin >> n;