Improve inv

Author: adamant-pwn

cp-algo/math/fft.hpp

@@ -7,6 +7,7 @@
 #include <cassert>
 #include <ranges>
 #include <vector>
+#include <deque>
 #include <bit>
 
 namespace cp_algo::math::fft {
@@ -249,8 +250,8 @@ namespace cp_algo::math::fft {
         void mul(auto &&B, auto& res, size_t k) {
             mul(B.A, B.B, res, k);
         }
-        std::vector<base> operator *= (auto &&B) {
-            std::vector<base> res(2 * A.size());
+        std::deque<base> operator *= (auto &&B) {
+            std::deque<base> res(2 * A.size());
             mul(B.A, B.B, res, size(res));
             return res;
         }

cp-algo/math/poly.hpp

@@ -538,7 +538,7 @@ namespace cp_algo::math {
         }
 
         poly_t x2() { // P(x) -> P(x^2)
-            std::vector<T> res(2 * a.size());
+            std::deque<T> res(2 * a.size());
             for(size_t i = 0; i < a.size(); i++) {
                 res[2 * i] = a[i];
             }
@@ -583,8 +583,12 @@ namespace cp_algo::math {
         }
 
         // inverse series mod x^n
+        poly_t& inv_inplace(size_t n) {
+            poly::impl::inv_inplace(*this, n);
+            return *this;
+        }
         poly_t inv(size_t n) const {
-            return poly::impl::inv(*this, n);
+            return poly_t(*this).inv_inplace(n);
         }
         // [x^k]..[x^{k+n-1}] of inv()
         // supports negative k if k+n >= 0

cp-algo/math/poly/impl/div.hpp

@@ -3,6 +3,7 @@
 #include "../../fft.hpp"
 #include "../../common.hpp"
 #include <cassert>
+#include <deque>
 // operations related to polynomial division
 namespace cp_algo::math::poly::impl {
     auto divmod_slow(auto const& p, auto const& q) {
@@ -90,27 +91,34 @@ namespace cp_algo::math::poly::impl {
         ).div_xk(2*q0.deg()).mod_xk(n);
     }
     template<typename poly>
-    poly inv(poly const& p, size_t n) {
+    void inv_inplace(poly&& p, size_t n) {
+        using poly_t = std::decay_t<poly>;
+        using base = poly_t::base;
         if(n == 1) {
-            return poly(1) / p[0];
+            p = base(1) / p[0];
+            return;
         }
         // Q(-x) = P0(x^2) + xP1(x^2)
         auto [q0, q1] = p.bisect(n);
-        q1 *= -1;
 
         int N = fft::com_size((n + 1) / 2, (n + 1) / 2);
 
-        auto q0f = fft::dft<typename poly::base>(q0.a, N);
-        auto q1f = fft::dft<typename poly::base>(q1.a, N);
+        auto q0f = fft::dft<base>(q0.a, N);
+        auto q1f = fft::dft<base>(q1.a, N);
 
         // Q(x)*Q(-x) = Q0(x^2)^2 - x^2 Q1(x^2)^2
-        auto qqf = fft::dft<typename poly::base>(inv(
-            poly(q0f * q0f) - poly(q1f * q1f).mul_xk(1)
-        , (n + 1) / 2).a, N);
+        auto qq = poly_t(q0f * q0f) - poly_t(q1f * q1f).mul_xk(1);
+        inv_inplace(qq, (n + 1) / 2);
+        auto qqf = fft::dft<base>(qq.a, N);
 
-        return (
-            poly(q0f * qqf).x2() + poly(q1f * qqf).x2().mul_xk(1)
-        ).mod_xk(n);
+        auto A = q0f * qqf;
+        auto B = q1f * qqf;
+        p.a.resize(n + 1);
+        for(size_t i = 0; i < n; i += 2) {
+            p.a[i] = A[i / 2];
+            p.a[i + 1] = -B[i / 2];
+        }
+        p.a.pop_back();
     }
 }
 #endif // CP_ALGO_MATH_POLY_IMPL_DIV_HPP

verify/poly/inv.test.cpp

@@ -16,7 +16,7 @@ void solve() {
     cin >> n;
     vector<base> a(n);
     copy_n(istream_iterator<base>(cin), n, begin(a));
-    polyn(a).inv(n).print(n);
+    polyn(a).inv_inplace(n).print(n);
 }
 
 signed main() {