Improve Frobenius normal form computation

Author: adamant-pwn

cp-algo/linalg/frobenius.hpp

@@ -39,14 +39,9 @@ namespace cp_algo::linalg {
                 if(full_rec.mod_xk(start) != polyn()) {
                     auto charp = full_rec.div_xk(start);
                     auto x = basis_init[start];
-                    start = 0;
-                    for(auto &rec: charps) {
-                        polyn cur_rec = full_rec.substr(start, rec.deg());
-                        auto shift = cur_rec / charp;
-                        for(int j = 0; j <= shift.deg(); j++) {
-                            x.add_scaled(basis_init[start + j], shift[j]);
-                        }
-                        start += rec.deg();
+                    auto shift = full_rec / charp;
+                    for(int j = 0; j < shift.deg(); j++) {
+                        x.add_scaled(basis_init[j], shift[j]);
                     }
                     basis.resize(start);
                     basis_init.resize(start);
@@ -75,13 +70,12 @@ namespace cp_algo::linalg {
     }
 
     template<typename base>
-    auto frobenius_pow(matrix<base> A, uint64_t k) {
-        using polyn = math::poly_t<base>;
+    auto with_frobenius(matrix<base> const& A, auto &&callback) {
         auto [T, Tinv, charps] = frobenius_form<full>(A);
         std::vector<matrix<base>> blocks;
         for(auto charp: charps) {
             matrix<base> block(charp.deg());
-            auto xk = polyn::xk(1).powmod(k, charp);
+            auto xk = callback(charp);
             for(size_t i = 0; i < block.n(); i++) {
                 std::ranges::copy(xk.a, begin(block[i]));
                 xk = xk.mul_xk(1) % charp;
@@ -91,5 +85,12 @@ namespace cp_algo::linalg {
         auto S = matrix<base>::block_diagonal(blocks);
         return Tinv * S * T;
     }
+
+    template<typename base>
+    auto frobenius_pow(matrix<base> const& A, uint64_t k) {
+        return with_frobenius(A, [k](auto const& charp) {
+            return math::poly_t<base>::xk(1).powmod(k, charp);
+        });
+    }
 };
 #endif // CP_ALGO_LINALG_FROBENIUS_HPP

cp-algo/linalg/vector.hpp

@@ -56,8 +56,10 @@ namespace cp_algo::linalg {
         }
 
         virtual void add_scaled(vec const& b, base scale, size_t i = 0) {
-            for(; i < size(*this); i++) {
-                (*this)[i] += scale * b[i];
+            if(scale != base(0)) {
+                for(; i < size(*this); i++) {
+                    (*this)[i] += scale * b[i];
+                }
             }
         }
         virtual vec const& normalize() {
@@ -125,14 +127,16 @@ namespace cp_algo::linalg {
         using Base::Base;
 
         void add_scaled(vec const& b, base scale, size_t i = 0) override {
-            for(; i < size(*this); i++) {
-                (*this)[i].add_unsafe(scale.getr() * b[i].getr());
-            }
-            if(++counter == 8) {
-                for(auto &it: *this) {
-                    it.pseudonormalize();
+            if(scale != base(0)) {
+                for(; i < size(*this); i++) {
+                    (*this)[i].add_unsafe(scale.getr() * b[i].getr());
+                }
+                if(++counter == 8) {
+                    for(auto &it: *this) {
+                        it.pseudonormalize();
+                    }
+                    counter = 0;
                 }
-                counter = 0;
             }
         }
         vec const& normalize() override {

verify/linalg/characteristic.test.cpp

@@ -1,8 +1,8 @@
 // @brief Characteristic Polynomial
 #define PROBLEM "https://judge.yosupo.jp/problem/characteristic_polynomial"
 #pragma GCC optimize("Ofast,unroll-loops")
-#pragma GCC target("tune=native")
-#define CP_ALGO_MAXN 1 << 12
+//#pragma GCC target("tune=native")
+#define CP_ALGO_MAXN 1 << 10
 #include "cp-algo/linalg/frobenius.hpp"
 #include <bits/stdc++.h>
 

verify/linalg/pow_fast.test.cpp

@@ -2,7 +2,7 @@
 #define PROBLEM "https://judge.yosupo.jp/problem/pow_of_matrix"
 #pragma GCC optimize("Ofast,unroll-loops")
 #pragma GCC target("tune=native")
-#define CP_ALGO_MAXN 1 << 12
+#define CP_ALGO_MAXN 1 << 10
 #include "cp-algo/linalg/frobenius.hpp"
 #include <bits/stdc++.h>