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: Convolution mod $10^9+7$ (verify/poly/convolution107.test.cpp)

Depends on

Code

// @brief Convolution mod $10^9+7$
#define PROBLEM "https://judge.yosupo.jp/problem/convolution_mod_1000000007"
#pragma GCC optimize("Ofast,unroll-loops")
#include "cp-algo/math/fft.hpp"
#include <bits/stdc++.h>

using namespace std;
using namespace cp_algo::math;

const int mod = 1e9 + 7;
using base = modint<mod>;

void solve() {
    int n, m;
    cin >> n >> m;
    vector<base> a(n), b(m);
    copy_n(istream_iterator<base>(cin), n, begin(a));
    copy_n(istream_iterator<base>(cin), m, begin(b));
    fft::mul(a, b);
    ranges::copy(a, ostream_iterator<base>(cout, " "));
}

signed main() {
    //freopen("input.txt", "r", stdin);
    ios::sync_with_stdio(0);
    cin.tie(0);
    int t = 1;
    while(t--) {
        solve();
    }
}

#line 1 "verify/poly/convolution107.test.cpp"
// @brief Convolution mod $10^9+7$
#define PROBLEM "https://judge.yosupo.jp/problem/convolution_mod_1000000007"
#pragma GCC optimize("Ofast,unroll-loops")
#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/util/checkpoint.hpp"


#line 4 "cp-algo/util/checkpoint.hpp"
#include <chrono>
#include <string>
namespace cp_algo {
    template<bool final = false>
    void checkpoint([[maybe_unused]] std::string const& msg = "") {
#ifdef CP_ALGO_CHECKPOINT
        static double last = 0;
        double now = (double)clock() / CLOCKS_PER_SEC;
        double delta = now - last;
        last = now;
        if(msg.size()) {
            std::cerr << msg << ": " << (final ? now : delta) * 1000 << " ms\n";
        }
#endif
    }
}

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


#line 1 "cp-algo/util/complex.hpp"


#line 4 "cp-algo/util/complex.hpp"
#include <cmath>
namespace cp_algo {
    // Custom implementation, since std::complex is UB on non-floating types
    template<typename T>
    struct complex {
        using value_type = T;
        T x, y;
        constexpr complex() {}
        constexpr complex(T x): x(x), y() {}
        constexpr complex(T x, T y): x(x), y(y) {}
        complex& operator *= (T t) {x *= t; y *= t; return *this;}
        complex& operator /= (T t) {x /= t; y /= t; return *this;}
        complex operator * (T t) const {return complex(*this) *= t;}
        complex operator / (T t) const {return complex(*this) /= t;}
        complex& operator += (complex t) {x += t.x; y += t.y; return *this;}
        complex& operator -= (complex t) {x -= t.x; y -= t.y; return *this;}
        complex operator * (complex t) const {return {x * t.x - y * t.y, x * t.y + y * t.x};}
        complex operator / (complex t) const {return *this * t.conj() / t.norm();}
        complex operator + (complex t) const {return complex(*this) += t;}
        complex operator - (complex t) const {return complex(*this) -= t;}
        complex& operator *= (complex t) {return *this = *this * t;}
        complex& operator /= (complex t) {return *this = *this / t;}
        complex operator - () const {return {-x, -y};}
        complex conj() const {return {x, -y};}
        T norm() const {return x * x + y * y;}
        T abs() const {return std::sqrt(norm());}
        T real() const {return x;}
        T imag() const {return y;}
        T& real() {return x;}
        T& imag() {return y;}
        static constexpr complex polar(T r, T theta) {return {r * cos(theta), r * sin(theta)};}
        auto operator <=> (complex const& t) const = default;
    };
    template<typename T>
    complex<T> operator * (auto x, complex<T> y) {return y *= x;}
    template<typename T> complex<T> conj(complex<T> x) {return x.conj();}
    template<typename T> T norm(complex<T> x) {return x.norm();}
    template<typename T> T abs(complex<T> x) {return x.abs();}
    template<typename T> T& real(complex<T> &x) {return x.real();}
    template<typename T> T& imag(complex<T> &x) {return x.imag();}
    template<typename T> T real(complex<T> const& x) {return x.real();}
    template<typename T> T imag(complex<T> const& x) {return x.imag();}
    template<typename T>
    constexpr complex<T> polar(T r, T theta) {
        return complex<T>::polar(r, theta);
    }
    template<typename T>
    std::ostream& operator << (std::ostream &out, complex<T> x) {
        return out << x.real() << ' ' << x.imag();
    }
}

#line 1 "cp-algo/util/new_big.hpp"


#include <sys/mman.h>
namespace cp_algo {
    template<typename T>
    auto new_big(size_t len) {
        auto raw = mmap(nullptr, len * sizeof(T),
            PROT_READ | PROT_WRITE,
            MAP_PRIVATE | MAP_ANONYMOUS,
            -1, 0);
        madvise(raw, len * sizeof(T), MADV_HUGEPAGE);
        madvise(raw, len * sizeof(T), MADV_POPULATE_WRITE);
        return static_cast<T*>(raw);
    }
    template<typename T>
    void delete_big(T* ptr, size_t len) {
        munmap(ptr, len * sizeof(T));
    }
}

#line 6 "cp-algo/math/cvector.hpp"
#include <experimental/simd>
#include <ranges>

namespace stdx = std::experimental;
namespace cp_algo::math::fft {
    using ftype = double;
    static constexpr size_t bytes = 32;
    static constexpr size_t flen = bytes / sizeof(ftype);
    using point = complex<ftype>;
    using vftype [[gnu::vector_size(bytes)]] = ftype;
    using vpoint = complex<vftype>;
    static constexpr vftype vz = {};
    static constexpr vpoint vi = {vz, vz + 1};

    struct cvector {
        vpoint *r;
        size_t sz;
        cvector(size_t n) {
            sz = std::max(flen, std::bit_ceil(n));
            r = new_big<vpoint>(sz / flen);
            checkpoint("cvector create");
        }
        cvector(cvector const& t) {
            sz = t.sz;
            r = new_big<vpoint>(sz / flen);
            memcpy(r, t.r, (sz / flen) * sizeof(vpoint));
            checkpoint("cvector copy");
        }
        cvector(cvector&& t) noexcept {
            sz = t.sz;
            r = std::exchange(t.r, nullptr);
        }
        ~cvector() noexcept {
            if(r) {
                delete_big(r, sz / flen);
            }
        }

        vpoint& at(size_t k) {return r[k / flen];}
        vpoint at(size_t k) const {return r[k / flen];}
        template<class pt = point>
        void set(size_t k, pt t) {
            if constexpr(std::is_same_v<pt, point>) {
                real(r[k / flen])[k % flen] = real(t);
                imag(r[k / flen])[k % flen] = imag(t);
            } else {
                at(k) = t;
            }
        }
        template<class pt = point>
        pt get(size_t k) const {
            if constexpr(std::is_same_v<pt, point>) {
                return {real(r[k / flen])[k % flen], imag(r[k / flen])[k % flen]};
            } else {
                return at(k);
            }
        }

        size_t size() const {
            return sz;
        }
        static size_t eval_arg(size_t n) {
            if(n < pre_roots) {
                return eval_args[n];
            } else {
                return eval_arg(n / 2) | (n & 1) << (std::bit_width(n) - 1);
            }
        }
        static auto root(size_t n, size_t k) {
            if(n < pre_roots) {
                return roots[n + k];
            } else {
                return polar(1., std::numbers::pi / (ftype)n * (ftype)k);
            }
        }
        static point eval_point(size_t n) {
            if(n < pre_roots) {
                return evalp[n];
            } else {
                return root(2 * std::bit_floor(n), eval_arg(n));
            }
        }
        static void exec_on_roots(size_t n, size_t m, auto &&callback) {
            point cur;
            point arg = root(n, 1);
            for(size_t i = 0; i < m; i++) {
                if(i % 32 == 0 || n < pre_roots) {
                    cur = root(n, i);
                } else {
                    cur *= arg;
                }
                callback(i, cur);
            }
        }
        template<int step = 1>
        static void exec_on_evals(size_t n, auto &&callback) {
            for(size_t i = 0; i < n; i++) {
                callback(i, eval_point(step * i));
            }
        }
        static auto dot_block(size_t k, cvector const& A, cvector const& B) {
            auto rt = eval_point(k / flen / 2);
            if(k / flen % 2) {
                rt = -rt;
            }
            auto [Ax, Ay] = A.at(k);
            auto Bv = B.at(k);
            vpoint res = vz;
            for (size_t i = 0; i < flen; i++) {
                res += vpoint(vz + Ax[i], vz + Ay[i]) * Bv;
                real(Bv) = __builtin_shufflevector(real(Bv), real(Bv), 3, 0, 1, 2);
                imag(Bv) = __builtin_shufflevector(imag(Bv), imag(Bv), 3, 0, 1, 2);
                auto x = real(Bv)[0], y = imag(Bv)[0];
                real(Bv)[0] = x * real(rt) - y * imag(rt);
                imag(Bv)[0] = x * imag(rt) + y * real(rt);
            }
            return res;
        }

        void dot(cvector const& t) {
            size_t n = this->size();
            for(size_t k = 0; k < n; k += flen) {
                set(k, dot_block(k, *this, t));
            }
            checkpoint("dot");
        }

        void ifft() {
            size_t n = size();
            for(size_t i = flen; i <= n / 2; i *= 2) {
                if (4 * i <= n) { // radix-4
                    exec_on_evals<2>(n / (4 * i), [&](size_t k, point rt) {
                        k *= 4 * i;
                        vpoint v1 = {vz + real(rt), vz - imag(rt)};
                        vpoint v2 = v1 * v1;
                        vpoint v3 = v1 * v2;
                        for(size_t j = k; j < k + i; j += flen) {
                            auto A = at(j);
                            auto B = at(j + i);
                            auto C = at(j + 2 * i);
                            auto D = at(j + 3 * i);
                            at(j) = (A + B + C + D);
                            at(j + 2 * i) = (A + B - C - D) * v2;
                            at(j +     i) = (A - B - vi * (C - D)) * v1;
                            at(j + 3 * i) = (A - B + vi * (C - D)) * v3;
                        }
                    });
                    i *= 2;
                } else { // radix-2 fallback
                    exec_on_evals(n / (2 * i), [&](size_t k, point rt) {
                        k *= 2 * i;
                        vpoint cvrt = {vz + real(rt), vz - imag(rt)};
                        for(size_t j = k; j < k + i; j += flen) {
                            auto B = at(j) - at(j + i);
                            at(j) += at(j + i);
                            at(j + i) = B * cvrt;
                        }
                    });
                }
            }
            checkpoint("ifft");
            for(size_t k = 0; k < n; k += flen) {
                set(k, get<vpoint>(k) /= vz + (ftype)(n / flen));
            }
        }
        void fft() {
            size_t n = size();
            for(size_t i = n / 2; i >= flen; i /= 2) {
                if (i / 2 >= flen) { // radix-4
                    i /= 2;
                    exec_on_evals<2>(n / (4 * i), [&](size_t k, point rt) {
                        k *= 4 * i;
                        vpoint v1 = {vz + real(rt), vz + imag(rt)};
                        vpoint v2 = v1 * v1;
                        vpoint v3 = v1 * v2;
                        for(size_t j = k; j < k + i; j += flen) {
                            auto A = at(j);
                            auto B = at(j + i) * v1;
                            auto C = at(j + 2 * i) * v2;
                            auto D = at(j + 3 * i) * v3;
                            at(j)         = (A + C) + (B + D);
                            at(j + i)     = (A + C) - (B + D);
                            at(j + 2 * i) = (A - C) + vi * (B - D);
                            at(j + 3 * i) = (A - C) - vi * (B - D);
                        }
                    });
                } else { // radix-2 fallback
                    exec_on_evals(n / (2 * i), [&](size_t k, point rt) {
                        k *= 2 * i;
                        vpoint vrt = {vz + real(rt), vz + imag(rt)};
                        for(size_t j = k; j < k + i; j += flen) {
                            auto t = at(j + i) * vrt;
                            at(j + i) = at(j) - t;
                            at(j) += t;
                        }
                    });
                }
            }
            checkpoint("fft");
        }
        static constexpr size_t pre_roots = 1 << 16;
        static constexpr std::array<point, pre_roots> roots = []() {
            std::array<point, pre_roots> res = {};
            for(size_t n = 1; n < res.size(); n *= 2) {
                for(size_t k = 0; k < n; k++) {
                    res[n + k] = polar(1., std::numbers::pi / ftype(n) * ftype(k));
                }
            }
            return res;
        }();
        static constexpr std::array<size_t, pre_roots> eval_args = []() {
            std::array<size_t, pre_roots> res = {};
            for(size_t i = 1; i < pre_roots; i++) {
                res[i] = res[i >> 1] | (i & 1) << (std::bit_width(i) - 1);
            }
            return res;
        }();
        static constexpr std::array<point, pre_roots> evalp = []() {
            std::array<point, pre_roots> res = {};
            res[0] = 1;
            for(size_t n = 1; n < pre_roots; n++) {
                res[n] = polar(1., std::numbers::pi * ftype(eval_args[n]) / ftype(2 * std::bit_floor(n)));
            }
            return res;
        }();
    };
}

#line 8 "cp-algo/math/fft.hpp"
namespace cp_algo::math::fft {
    template<modint_type base>
    struct dft {
        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, std::min(n, size(a)), [&](size_t i, auto rt) {
                auto splt = [&](size_t i) {
                    auto ai = ftype(i < size(a) ? a[i].rem() : 0);
                    auto rem = std::remainder(ai, split);
                    auto quo = (ai - rem) / split;
                    return std::pair{rem, quo};
                };
                auto [rai, qai] = splt(i);
                auto [rani, qani] = splt(n + i);
                A.set(i, point(rai, rani) * rt);
                B.set(i, point(qai, qani) * rt);
            });
            checkpoint("dft init");
            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 k = 0; k < n; k += flen) {
                auto rt = cvector::eval_point(k / flen / 2);
                if(k / flen % 2) {
                    rt = -rt;
                }
                auto [Ax, Ay] = A.at(k);
                auto [Bx, By] = B.at(k);
                vpoint AC, AD, BC, BD;
                AC = AD = BC = BD = vz;
                auto Cv = C.at(k), Dv = D.at(k);
                for (size_t i = 0; i < flen; i++) {
                    vpoint Av = {vz + Ax[i], vz + Ay[i]}, Bv = {vz + Bx[i], vz + By[i]};
                    AC += Av * Cv; AD += Av * Dv;
                    BC += Bv * Cv; BD += Bv * Dv;
                    real(Cv) = __builtin_shufflevector(real(Cv), real(Cv), 3, 0, 1, 2);
                    imag(Cv) = __builtin_shufflevector(imag(Cv), imag(Cv), 3, 0, 1, 2);
                    real(Dv) = __builtin_shufflevector(real(Dv), real(Dv), 3, 0, 1, 2);
                    imag(Dv) = __builtin_shufflevector(imag(Dv), imag(Dv), 3, 0, 1, 2);
                    auto cx = real(Cv)[0], cy = imag(Cv)[0];
                    auto dx = real(Dv)[0], dy = imag(Dv)[0];
                    real(Cv)[0] = cx * real(rt) - cy * imag(rt);
                    imag(Cv)[0] = cx * imag(rt) + cy * real(rt);
                    real(Dv)[0] = dx * real(rt) - dy * imag(rt);
                    imag(Dv)[0] = dx * imag(rt) + dy * real(rt);
                }
                A.at(k) = AC;
                C.at(k) = AD + BC;
                B.at(k) = BD;
            }
            checkpoint("dot");
            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;
            });
            checkpoint("recover mod");
        }
        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);
        auto A = dft<base>(a | std::views::take(k), n);
        a.resize(k);
        checkpoint("resize a");
        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 5 "verify/poly/convolution107.test.cpp"
#include <bits/stdc++.h>

using namespace std;
using namespace cp_algo::math;

const int mod = 1e9 + 7;
using base = modint<mod>;

void solve() {
    int n, m;
    cin >> n >> m;
    vector<base> a(n), b(m);
    copy_n(istream_iterator<base>(cin), n, begin(a));
    copy_n(istream_iterator<base>(cin), m, begin(b));
    fft::mul(a, b);
    ranges::copy(a, ostream_iterator<base>(cout, " "));
}

signed main() {
    //freopen("input.txt", "r", stdin);
    ios::sync_with_stdio(0);
    cin.tie(0);
    int t = 1;
    while(t--) {
        solve();
    }
}

Test cases

Env Name Status Elapsed Memory
g++ example_00 :heavy_check_mark: AC 5 ms 4 MB
g++ example_01 :heavy_check_mark: AC 4 ms 4 MB
g++ fft_killer_00 :heavy_check_mark: AC 158 ms 44 MB
g++ fft_killer_01 :heavy_check_mark: AC 161 ms 44 MB
g++ fft_killer_02 :heavy_check_mark: AC 148 ms 44 MB
g++ fft_killer_03 :heavy_check_mark: AC 157 ms 43 MB
g++ fft_killer_04 :heavy_check_mark: AC 151 ms 43 MB
g++ fft_killer_05 :heavy_check_mark: AC 153 ms 44 MB
g++ fft_killer_06 :heavy_check_mark: AC 159 ms 44 MB
g++ fft_killer_07 :heavy_check_mark: AC 164 ms 44 MB
g++ fft_killer_08 :heavy_check_mark: AC 146 ms 43 MB
g++ fft_killer_09 :heavy_check_mark: AC 154 ms 43 MB
g++ max_ans_zero_00 :heavy_check_mark: AC 170 ms 44 MB
g++ max_random_00 :heavy_check_mark: AC 156 ms 44 MB
g++ max_random_01 :heavy_check_mark: AC 173 ms 43 MB
g++ medium_00 :heavy_check_mark: AC 7 ms 4 MB
g++ medium_01 :heavy_check_mark: AC 6 ms 4 MB
g++ medium_02 :heavy_check_mark: AC 6 ms 4 MB
g++ medium_all_zero_00 :heavy_check_mark: AC 6 ms 4 MB
g++ random_00 :heavy_check_mark: AC 131 ms 42 MB
g++ random_01 :heavy_check_mark: AC 136 ms 42 MB
g++ random_02 :heavy_check_mark: AC 69 ms 24 MB
g++ signed_overflow_00 :heavy_check_mark: AC 6 ms 4 MB
g++ small_00 :heavy_check_mark: AC 5 ms 4 MB
g++ small_01 :heavy_check_mark: AC 5 ms 4 MB
g++ small_02 :heavy_check_mark: AC 5 ms 4 MB
g++ small_03 :heavy_check_mark: AC 5 ms 4 MB
g++ small_04 :heavy_check_mark: AC 5 ms 4 MB
g++ small_05 :heavy_check_mark: AC 5 ms 4 MB
g++ small_06 :heavy_check_mark: AC 5 ms 4 MB
g++ small_07 :heavy_check_mark: AC 5 ms 4 MB
g++ small_08 :heavy_check_mark: AC 5 ms 4 MB
g++ small_09 :heavy_check_mark: AC 5 ms 4 MB
g++ small_10 :heavy_check_mark: AC 5 ms 4 MB
g++ small_11 :heavy_check_mark: AC 5 ms 4 MB
g++ small_12 :heavy_check_mark: AC 5 ms 4 MB
g++ small_13 :heavy_check_mark: AC 4 ms 4 MB
g++ small_14 :heavy_check_mark: AC 5 ms 4 MB
g++ small_15 :heavy_check_mark: AC 5 ms 4 MB
g++ small_and_large_00 :heavy_check_mark: AC 95 ms 41 MB
g++ small_and_large_01 :heavy_check_mark: AC 100 ms 41 MB
g++ small_and_large_02 :heavy_check_mark: AC 99 ms 41 MB
g++ small_and_large_03 :heavy_check_mark: AC 93 ms 41 MB
g++ unsigned_overflow_00 :heavy_check_mark: AC 5 ms 4 MB
Back to top page