CP-Algorithms Library

This documentation is automatically generated by competitive-verifier/competitive-verifier

View the Project on GitHub cp-algorithms/cp-algorithms-aux

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

Depends on

Required by

Verified with

Code

#ifndef CP_ALGO_MATH_POLY_IMPL_DIV_HPP
#define CP_ALGO_MATH_POLY_IMPL_DIV_HPP
#include "../../fft.hpp"
#include "../../common.hpp"
#include <cassert>
// operations related to polynomial division
namespace cp_algo::math::poly::impl {
    auto divmod_slow(auto const& p, auto const& q) {
        auto R = p;
        auto D = decltype(p){};
        auto q_lead_inv = q.lead().inv();
        while(R.deg() >= q.deg()) {
            D.a.push_back(R.lead() * q_lead_inv);
            if(D.lead() != 0) {
                for(size_t i = 1; i <= q.a.size(); i++) {
                    R.a[R.a.size() - i] -= D.lead() * q.a[q.a.size() - i];
                }
            }
            R.a.pop_back();
        }
        std::ranges::reverse(D.a);
        R.normalize();
        return std::array{D, R};
    }
    template<typename poly>
    auto divmod_hint(poly const& p, poly const& q, poly const& qri) {
        assert(!q.is_zero());
        int d = p.deg() - q.deg();
        if(std::min(d, q.deg()) < magic) {
            return divmod_slow(p, q);
        }
        poly D;
        if(d >= 0) {
            D = (p.reversed().mod_xk(d + 1) * qri.mod_xk(d + 1)).mod_xk(d + 1).reversed(d + 1);
        }
        return std::array{D, p - D * q};
    }
    auto divmod(auto const& p, auto const& q) {
        assert(!q.is_zero());
        int d = p.deg() - q.deg();
        if(std::min(d, q.deg()) < magic) {
            return divmod_slow(p, q);
        }
        return divmod_hint(p, q, q.reversed().inv(d + 1));
    }

    template<typename poly>
    poly powmod_hint(poly const& p, int64_t k, poly const& md, poly const& mdri) {
        return bpow(p, k, poly(1), [&](auto const& p, auto const& q){
            return divmod_hint(p * q, md, mdri)[1];
        });
    }
    template<typename poly>
    auto powmod(poly const& p, int64_t k, poly const& md) {
        int d = md.deg();
        if(p == poly::xk(1) && false) { // does it actually speed anything up?..
            if(k < md.deg()) {
                return poly::xk(k);
            } else {
                auto mdr = md.reversed();
                return (mdr.inv(k - md.deg() + 1, md.deg()) * mdr).reversed(md.deg());
            }
        }
        if(md == poly::xk(d)) {
            return p.pow(k, d);
        }
        if(md == poly::xk(d) - poly(1)) {
            return p.powmod_circular(k, d);
        }
        return powmod_hint(p, k, md, md.reversed().inv(md.deg() + 1));
    }
    template<typename poly>
    poly& inv_inplace(poly& q, int64_t k, size_t n) {
        using poly_t = std::decay_t<poly>;
        using base = poly_t::base;
        if(k <= std::max<int64_t>(n, size(q.a))) {
            return q.inv_inplace(k + n).div_xk_inplace(k);
        }
        if(k % 2) {
            return inv_inplace(q, k - 1, n + 1).div_xk_inplace(1);
        }
        auto [q0, q1] = q.bisect();
        auto qq = q0 * q0 - (q1 * q1).mul_xk_inplace(1);
        inv_inplace(qq, k / 2 - q.deg() / 2, (n + 1) / 2 + q.deg() / 2);
        size_t N = fft::com_size(size(q0.a), size(qq.a));
        auto q0f = fft::dft<base>(q0.a, N);
        auto q1f = fft::dft<base>(q1.a, N);
        auto qqf = fft::dft<base>(qq.a, N);
        size_t M = q0.deg() + (n + 1) / 2;
        std::vector<base> A(M), B(M);
        q0f.mul(qqf, A, M);
        q1f.mul_inplace(qqf, B, M);
        q.a.resize(n + 1);
        for(size_t i = 0; i < n; i += 2) {
            q.a[i] = A[q0.deg() + i / 2];
            q.a[i + 1] = -B[q0.deg() + i / 2];
        }
        q.a.pop_back();
        q.normalize();
        return q;
    }
    template<typename poly>
    poly& inv_inplace(poly& p, size_t n) {
        using poly_t = std::decay_t<poly>;
        using base = poly_t::base;
        if(n == 1) {
            return p = base(1) / p[0];
        }
        // Q(-x) = P0(x^2) + xP1(x^2)
        auto [q0, q1] = p.bisect(n);
        
        size_t N = fft::com_size(size(q0.a), (n + 1) / 2);
        
        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 qq = poly_t(q0f * q0f) - poly_t(q1f * q1f).mul_xk_inplace(1);

        inv_inplace(qq, (n + 1) / 2);
        auto qqf = fft::dft<base>(qq.a, N);
        
        std::vector<base> A((n + 1) / 2), B((n + 1) / 2);
        q0f.mul(qqf, A, (n + 1) / 2);
        q1f.mul_inplace(qqf, B, (n + 1) / 2);
        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();
        p.normalize();
        return p;
    }
}
#endif // CP_ALGO_MATH_POLY_IMPL_DIV_HPP

#line 1 "cp-algo/math/poly/impl/div.hpp"


#line 1 "cp-algo/math/fft.hpp"


#line 1 "cp-algo/number_theory/modint.hpp"


#line 1 "cp-algo/math/common.hpp"


#include <functional>
#include <cstdint>
namespace cp_algo::math {
#ifdef CP_ALGO_MAXN
    const int maxn = CP_ALGO_MAXN;
#else
    const int maxn = 1 << 19;
#endif
    const int magic = 64; // threshold for sizes to run the naive algo

    auto bpow(auto const& x, auto n, auto const& one, auto op) {
        if(n == 0) {
            return one;
        } else {
            auto t = bpow(x, n / 2, one, op);
            t = op(t, t);
            if(n % 2) {
                t = op(t, x);
            }
            return t;
        }
    }
    auto bpow(auto x, auto n, auto ans) {
        return bpow(x, n, ans, std::multiplies{});
    }
    template<typename T>
    T bpow(T const& x, auto n) {
        return bpow(x, n, T(1));
    }
}

#line 4 "cp-algo/number_theory/modint.hpp"
#include <iostream>
#include <cassert>
namespace cp_algo::math {
    inline constexpr auto inv2(auto x) {
        assert(x % 2);
        std::make_unsigned_t<decltype(x)> y = 1;
        while(y * x != 1) {
            y *= 2 - x * y;
        }
        return y;
    }

    template<typename modint, typename _Int>
    struct modint_base {
        using Int = _Int;
        using UInt = std::make_unsigned_t<Int>;
        static constexpr size_t bits = sizeof(Int) * 8;
        using Int2 = std::conditional_t<bits <= 32, int64_t, __int128_t>;
        using UInt2 = std::conditional_t<bits <= 32, uint64_t, __uint128_t>;
        static Int mod() {
            return modint::mod();
        }
        static UInt imod() {
            return modint::imod();
        }
        static UInt2 pw128() {
            return modint::pw128();
        }
        static UInt m_reduce(UInt2 ab) {
            if(mod() % 2 == 0) [[unlikely]] {
                return UInt(ab % mod());
            } else {
                UInt2 m = (UInt)ab * imod();
                return UInt((ab + m * mod()) >> bits);
            }
        }
        static UInt m_transform(UInt a) {
            if(mod() % 2 == 0) [[unlikely]] {
                return a;
            } else {
                return m_reduce(a * pw128());
            }
        }
        modint_base(): r(0) {}
        modint_base(Int2 rr): r(UInt(rr % mod())) {
            r = std::min(r, r + mod());
            r = m_transform(r);
        }
        modint inv() const {
            return bpow(to_modint(), mod() - 2);
        }
        modint operator - () const {
            modint neg;
            neg.r = std::min(-r, 2 * mod() - r);
            return neg;
        }
        modint& operator /= (const modint &t) {
            return to_modint() *= t.inv();
        }
        modint& operator *= (const modint &t) {
            r = m_reduce((UInt2)r * t.r);
            return to_modint();
        }
        modint& operator += (const modint &t) {
            r += t.r; r = std::min(r, r - 2 * mod());
            return to_modint();
        }
        modint& operator -= (const modint &t) {
            r -= t.r; r = std::min(r, r + 2 * mod());
            return to_modint();
        }
        modint operator + (const modint &t) const {return modint(to_modint()) += t;}
        modint operator - (const modint &t) const {return modint(to_modint()) -= t;}
        modint operator * (const modint &t) const {return modint(to_modint()) *= t;}
        modint operator / (const modint &t) const {return modint(to_modint()) /= t;}
        // Why <=> doesn't work?..
        auto operator == (const modint_base &t) const {return getr() == t.getr();}
        auto operator != (const modint_base &t) const {return getr() != t.getr();}
        auto operator <= (const modint_base &t) const {return getr() <= t.getr();}
        auto operator >= (const modint_base &t) const {return getr() >= t.getr();}
        auto operator < (const modint_base &t) const {return getr() < t.getr();}
        auto operator > (const modint_base &t) const {return getr() > t.getr();}
        Int rem() const {
            UInt R = getr();
            return 2 * R > (UInt)mod() ? R - mod() : R;
        }

        // Only use if you really know what you're doing!
        UInt modmod() const {return (UInt)8 * mod() * mod();};
        void add_unsafe(UInt t) {r += t;}
        void pseudonormalize() {r = std::min(r, r - modmod());}
        modint const& normalize() {
            if(r >= (UInt)mod()) {
                r %= mod();
            }
            return to_modint();
        }
        void setr(UInt rr) {r = m_transform(rr);}
        UInt getr() const {
            UInt res = m_reduce(r);
            return std::min(res, res - mod());
        }
        void setr_direct(UInt rr) {r = rr;}
        UInt getr_direct() const {return r;}
    private:
        UInt r;
        modint& to_modint() {return static_cast<modint&>(*this);}
        modint const& to_modint() const {return static_cast<modint const&>(*this);}
    };
    template<typename modint>
    concept modint_type = std::is_base_of_v<modint_base<modint, typename modint::Int>, modint>;
    template<modint_type modint>
    std::istream& operator >> (std::istream &in, modint &x) {
        typename modint::UInt r;
        auto &res = in >> r;
        x.setr(r);
        return res;
    }
    template<modint_type modint>
    std::ostream& operator << (std::ostream &out, modint const& x) {
        return out << x.getr();
    }

    template<auto m>
    struct modint: modint_base<modint<m>, decltype(m)> {
        using Base = modint_base<modint<m>, decltype(m)>;
        using Base::Base;
        static constexpr Base::UInt im = m % 2 ? inv2(-m) : 0;
        static constexpr Base::UInt r2 = (typename Base::UInt2)(-1) % m + 1;
        static constexpr Base::Int mod() {return m;}
        static constexpr Base::UInt imod() {return im;}
        static constexpr Base::UInt2 pw128() {return r2;}
    };

    template<typename Int = int64_t>
    struct dynamic_modint: modint_base<dynamic_modint<Int>, Int> {
        using Base = modint_base<dynamic_modint<Int>, Int>;
        using Base::Base;
        static Int mod() {return m;}
        static Base::UInt imod() {return im;}
        static Base::UInt2 pw128() {return r2;}
        static void switch_mod(Int nm) {
            m = nm;
            im = m % 2 ? inv2(-m) : 0;
            r2 = static_cast<Base::UInt>(static_cast<Base::UInt2>(-1) % m + 1);
        }

        // Wrapper for temp switching
        auto static with_mod(Int tmp, auto callback) {
            struct scoped {
                Int prev = mod();
                ~scoped() {switch_mod(prev);}
            } _;
            switch_mod(tmp);
            return callback();
        }
    private:
        static thread_local Int m;
        static thread_local Base::UInt im, r2;
    };
    template<typename Int>
    Int thread_local dynamic_modint<Int>::m = 1;
    template<typename Int>
    dynamic_modint<Int>::Base::UInt thread_local dynamic_modint<Int>::im = -1;
    template<typename Int>
    dynamic_modint<Int>::Base::UInt thread_local dynamic_modint<Int>::r2 = 0;
}

#line 1 "cp-algo/math/cvector.hpp"


#include <algorithm>
#line 5 "cp-algo/math/cvector.hpp"
#include <complex>
#include <vector>
#include <ranges>
namespace cp_algo::math::fft {
    using ftype = double;
    static constexpr size_t bytes = 32;
    static constexpr size_t flen = bytes / sizeof(ftype);
    using point = std::complex<ftype>;
    using vftype [[gnu::vector_size(bytes)]] = ftype;
    using vpoint = std::complex<vftype>;

#define WITH_IV(...)                             \
  [&]<size_t ... i>(std::index_sequence<i...>) { \
      return __VA_ARGS__;                        \
  }(std::make_index_sequence<flen>());

    template<typename ft>
    constexpr ft to_ft(auto x) {
        return ft{} + x;
    }
    template<typename pt>
    constexpr pt to_pt(point r) {
        using ft = std::conditional_t<std::is_same_v<point, pt>, ftype, vftype>;
        return {to_ft<ft>(r.real()), to_ft<ft>(r.imag())};
    }
    struct cvector {
        static constexpr size_t pre_roots = 1 << 17;
        std::vector<vftype> x, y;
        cvector(size_t n) {
            n = std::max(flen, std::bit_ceil(n));
            x.resize(n / flen);
            y.resize(n / flen);
        }
        template<class pt = point>
        void set(size_t k, pt t) {
            if constexpr(std::is_same_v<pt, point>) {
                x[k / flen][k % flen] = real(t);
                y[k / flen][k % flen] = imag(t);
            } else {
                x[k / flen] = real(t);
                y[k / flen] = imag(t);
            }
        }
        template<class pt = point>
        pt get(size_t k) const {
            if constexpr(std::is_same_v<pt, point>) {
                return {x[k / flen][k % flen], y[k / flen][k % flen]};
            } else {
                return {x[k / flen], y[k / flen]};
            }
        }
        vpoint vget(size_t k) const {
            return get<vpoint>(k);
        }

        size_t size() const {
            return flen * std::size(x);
        }
        void dot(cvector const& t) {
            size_t n = size();
            for(size_t k = 0; k < n; k += flen) {
                set(k, get<vpoint>(k) * t.get<vpoint>(k));
            }
        }
        static const cvector roots;
        template<class pt = point>
        static pt root(size_t n, size_t k) {
            if(n < pre_roots) {
                return roots.get<pt>(n + k);
            } else {
                auto arg = std::numbers::pi / ftype(n);
                if constexpr(std::is_same_v<pt, point>) {
                    return {cos(ftype(k) * arg), sin(ftype(k) * arg)};
                } else {
                    return WITH_IV(pt{vftype{cos(ftype(k + i) * arg)...},
                                      vftype{sin(ftype(k + i) * arg)...}});
                }
            }
        }
        template<class pt = point>
        static void exec_on_roots(size_t n, size_t m, auto &&callback) {
            size_t step = sizeof(pt) / sizeof(point);
            pt cur;
            pt arg = to_pt<pt>(root<point>(n, step));
            for(size_t i = 0; i < m; i += step) {
                if(i % 64 == 0 || n < pre_roots) {
                    cur = root<pt>(n, i);
                } else {
                    cur *= arg;
                }
                callback(i, cur);
            }
        }

        void ifft() {
            size_t n = size();
            for(size_t i = 1; i < n; i *= 2) {
                for(size_t j = 0; j < n; j += 2 * i) {
                    auto butterfly = [&]<class pt>(size_t k, pt rt) {
                        k += j;
                        auto t = get<pt>(k + i) * conj(rt);
                        set(k + i, get<pt>(k) - t);
                        set(k, get<pt>(k) + t);
                    };
                    if(2 * i <= flen) {
                        exec_on_roots(i, i, butterfly);
                    } else {
                        exec_on_roots<vpoint>(i, i, butterfly);
                    }
                }
            }
            for(size_t k = 0; k < n; k += flen) {
                set(k, get<vpoint>(k) /= to_pt<vpoint>(ftype(n)));
            }
        }
        void fft() {
            size_t n = size();
            for(size_t i = n / 2; i >= 1; i /= 2) {
                for(size_t j = 0; j < n; j += 2 * i) {
                    auto butterfly = [&]<class pt>(size_t k, pt rt) {
                        k += j;
                        auto A = get<pt>(k) + get<pt>(k + i);
                        auto B = get<pt>(k) - get<pt>(k + i);
                        set(k, A);
                        set(k + i, B * rt);
                    };
                    if(2 * i <= flen) {
                        exec_on_roots(i, i, butterfly);
                    } else {
                        exec_on_roots<vpoint>(i, i, butterfly);
                    }
                }
            }
        }
    };
    const cvector cvector::roots = []() {
        cvector res(pre_roots);
        for(size_t n = 1; n < res.size(); n *= 2) {
            auto base = std::polar(1., std::numbers::pi / ftype(n));
            point cur = 1;
            for(size_t k = 0; k < n; k++) {
                if((k & 15) == 0) {
                    cur = std::polar(1., std::numbers::pi * ftype(k) / ftype(n));
                }
                res.set(n + k, cur);
                cur *= base;
            }
        }
        return res;
    }();

    template<typename base>
    struct dft {
        cvector A;
        
        dft(std::vector<base> const& a, size_t n): A(n) {
            for(size_t i = 0; i < std::min(n, a.size()); i++) {
                A.set(i, a[i]);
            }
            if(n) {
                A.fft();
            }
        }

        std::vector<base> operator *= (dft const& B) {
            assert(A.size() == B.A.size());
            size_t n = A.size();
            if(!n) {
                return std::vector<base>();
            }
            A.dot(B.A);
            A.ifft();
            std::vector<base> res(n);
            for(size_t k = 0; k < n; k++) {
                res[k] = A.get(k);
            }
            return res;
        }

        auto operator * (dft const& B) const {
            return dft(*this) *= B;
        }

        point operator [](int i) const {return A.get(i);}
    };
}

#line 5 "cp-algo/math/fft.hpp"
namespace cp_algo::math::fft {
    template<modint_type base>
    struct dft<base> {
        int split;
        cvector A, B;
        
        dft(auto const& a, size_t n): A(n), B(n) {
            split = int(std::sqrt(base::mod())) + 1;
            cvector::exec_on_roots(2 * n, size(a), [&](size_t i, point rt) {
                size_t ti = std::min(i, i - n);
                auto rem = std::remainder(a[i].rem(), split);
                auto quo = (ftype(a[i].rem()) - rem) / split;
                A.set(ti, A.get(ti) + rem * rt);
                B.set(ti, B.get(ti) + quo * rt);
    
            });
            if(n) {
                A.fft();
                B.fft();
            }
        }

        void mul(auto &&C, auto const& D, auto &res, size_t k) {
            assert(A.size() == C.size());
            size_t n = A.size();
            if(!n) {
                res = {};
                return;
            }
            for(size_t i = 0; i < n; i += flen) {
                auto tmp = A.vget(i) * D.vget(i) + B.vget(i) * C.vget(i);
                A.set(i, A.vget(i) * C.vget(i));
                B.set(i, B.vget(i) * D.vget(i));
                C.set(i, tmp);
            }
            A.ifft();
            B.ifft();
            C.ifft();
            auto splitsplit = (base(split) * split).rem();
            cvector::exec_on_roots(2 * n, std::min(n, k), [&](size_t i, point rt) {
                rt = conj(rt);
                auto Ai = A.get(i) * rt;
                auto Bi = B.get(i) * rt;
                auto Ci = C.get(i) * rt;
                int64_t A0 = llround(real(Ai));
                int64_t A1 = llround(real(Ci));
                int64_t A2 = llround(real(Bi));
                res[i] = A0 + A1 * split + A2 * splitsplit;
                if(n + i >= k) {
                    return;
                }
                int64_t B0 = llround(imag(Ai));
                int64_t B1 = llround(imag(Ci));
                int64_t B2 = llround(imag(Bi));
                res[n + i] = B0 + B1 * split + B2 * splitsplit;
            });
        }
        void mul_inplace(auto &&B, auto& res, size_t k) {
            mul(B.A, B.B, res, k);
        }
        void mul(auto const& B, auto& res, size_t k) {
            mul(cvector(B.A), B.B, res, k);
        }
        std::vector<base> operator *= (dft &B) {
            std::vector<base> res(2 * A.size());
            mul_inplace(B, res, size(res));
            return res;
        }
        std::vector<base> operator *= (dft const& B) {
            std::vector<base> res(2 * A.size());
            mul(B, res, size(res));
            return res;
        }
        auto operator * (dft const& B) const {
            return dft(*this) *= B;
        }
        
        point operator [](int i) const {return A.get(i);}
    };
    
    void mul_slow(auto &a, auto const& b, size_t k) {
        if(empty(a) || empty(b)) {
            a.clear();
        } else {
            size_t n = std::min(k, size(a));
            size_t m = std::min(k, size(b));
            a.resize(k);
            for(int j = int(k - 1); j >= 0; j--) {
                a[j] *= b[0];
                for(int i = std::max(j - (int)n, 0) + 1; i < std::min(j + 1, (int)m); i++) {
                    a[j] += a[j - i] * b[i];
                }
            }
        }
    }
    size_t com_size(size_t as, size_t bs) {
        if(!as || !bs) {
            return 0;
        }
        return std::max(flen, std::bit_ceil(as + bs - 1) / 2);
    }
    void mul_truncate(auto &a, auto const& b, size_t k) {
        using base = std::decay_t<decltype(a[0])>;
        if(std::min({k, size(a), size(b)}) < magic) {
            mul_slow(a, b, k);
            return;
        }
        auto n = std::max(flen, std::bit_ceil(
            std::min(k, size(a)) + std::min(k, size(b)) - 1
        ) / 2);
        a.resize(k);
        auto A = dft<base>(a, n);
        if(&a == &b) {
            A.mul(A, a, k);
        } else {
            A.mul_inplace(dft<base>(b | std::views::take(k), n), a, k);
        }
    }
    void mul(auto &a, auto const& b) {
        if(size(a)) {
            mul_truncate(a, b, size(a) + size(b) - 1);
        }
    }
}

#line 6 "cp-algo/math/poly/impl/div.hpp"
// operations related to polynomial division
namespace cp_algo::math::poly::impl {
    auto divmod_slow(auto const& p, auto const& q) {
        auto R = p;
        auto D = decltype(p){};
        auto q_lead_inv = q.lead().inv();
        while(R.deg() >= q.deg()) {
            D.a.push_back(R.lead() * q_lead_inv);
            if(D.lead() != 0) {
                for(size_t i = 1; i <= q.a.size(); i++) {
                    R.a[R.a.size() - i] -= D.lead() * q.a[q.a.size() - i];
                }
            }
            R.a.pop_back();
        }
        std::ranges::reverse(D.a);
        R.normalize();
        return std::array{D, R};
    }
    template<typename poly>
    auto divmod_hint(poly const& p, poly const& q, poly const& qri) {
        assert(!q.is_zero());
        int d = p.deg() - q.deg();
        if(std::min(d, q.deg()) < magic) {
            return divmod_slow(p, q);
        }
        poly D;
        if(d >= 0) {
            D = (p.reversed().mod_xk(d + 1) * qri.mod_xk(d + 1)).mod_xk(d + 1).reversed(d + 1);
        }
        return std::array{D, p - D * q};
    }
    auto divmod(auto const& p, auto const& q) {
        assert(!q.is_zero());
        int d = p.deg() - q.deg();
        if(std::min(d, q.deg()) < magic) {
            return divmod_slow(p, q);
        }
        return divmod_hint(p, q, q.reversed().inv(d + 1));
    }

    template<typename poly>
    poly powmod_hint(poly const& p, int64_t k, poly const& md, poly const& mdri) {
        return bpow(p, k, poly(1), [&](auto const& p, auto const& q){
            return divmod_hint(p * q, md, mdri)[1];
        });
    }
    template<typename poly>
    auto powmod(poly const& p, int64_t k, poly const& md) {
        int d = md.deg();
        if(p == poly::xk(1) && false) { // does it actually speed anything up?..
            if(k < md.deg()) {
                return poly::xk(k);
            } else {
                auto mdr = md.reversed();
                return (mdr.inv(k - md.deg() + 1, md.deg()) * mdr).reversed(md.deg());
            }
        }
        if(md == poly::xk(d)) {
            return p.pow(k, d);
        }
        if(md == poly::xk(d) - poly(1)) {
            return p.powmod_circular(k, d);
        }
        return powmod_hint(p, k, md, md.reversed().inv(md.deg() + 1));
    }
    template<typename poly>
    poly& inv_inplace(poly& q, int64_t k, size_t n) {
        using poly_t = std::decay_t<poly>;
        using base = poly_t::base;
        if(k <= std::max<int64_t>(n, size(q.a))) {
            return q.inv_inplace(k + n).div_xk_inplace(k);
        }
        if(k % 2) {
            return inv_inplace(q, k - 1, n + 1).div_xk_inplace(1);
        }
        auto [q0, q1] = q.bisect();
        auto qq = q0 * q0 - (q1 * q1).mul_xk_inplace(1);
        inv_inplace(qq, k / 2 - q.deg() / 2, (n + 1) / 2 + q.deg() / 2);
        size_t N = fft::com_size(size(q0.a), size(qq.a));
        auto q0f = fft::dft<base>(q0.a, N);
        auto q1f = fft::dft<base>(q1.a, N);
        auto qqf = fft::dft<base>(qq.a, N);
        size_t M = q0.deg() + (n + 1) / 2;
        std::vector<base> A(M), B(M);
        q0f.mul(qqf, A, M);
        q1f.mul_inplace(qqf, B, M);
        q.a.resize(n + 1);
        for(size_t i = 0; i < n; i += 2) {
            q.a[i] = A[q0.deg() + i / 2];
            q.a[i + 1] = -B[q0.deg() + i / 2];
        }
        q.a.pop_back();
        q.normalize();
        return q;
    }
    template<typename poly>
    poly& inv_inplace(poly& p, size_t n) {
        using poly_t = std::decay_t<poly>;
        using base = poly_t::base;
        if(n == 1) {
            return p = base(1) / p[0];
        }
        // Q(-x) = P0(x^2) + xP1(x^2)
        auto [q0, q1] = p.bisect(n);
        
        size_t N = fft::com_size(size(q0.a), (n + 1) / 2);
        
        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 qq = poly_t(q0f * q0f) - poly_t(q1f * q1f).mul_xk_inplace(1);

        inv_inplace(qq, (n + 1) / 2);
        auto qqf = fft::dft<base>(qq.a, N);
        
        std::vector<base> A((n + 1) / 2), B((n + 1) / 2);
        q0f.mul(qqf, A, (n + 1) / 2);
        q1f.mul_inplace(qqf, B, (n + 1) / 2);
        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();
        p.normalize();
        return p;
    }
}
Back to top page