Convolution (Large)
(library/convolution/convolution_large.hpp)

• View this file on GitHub
• Last update: 2023-05-11 13:19:06+09:00
• Include: #include "library/convolution/convolution_large.hpp"

Convolution (Large)

$\mathrm{mod}\ p$ における $1$ の $2 ^ k$ 乗根が存在しない場合に、列を分割して畳み込む。

Depends on

• Naive Convolution (library/convolution/convolution_naive.hpp)

Verified with

• test/src/convolution/convolution_large/convolution_large.test.cpp

Code

#ifndef SUISEN_CONVOLTION_LARGE
#define SUISEN_CONVOLTION_LARGE

#include <atcoder/convolution>

#include "library/convolution/convolution_naive.hpp"

namespace suisen {
    template <typename mint, atcoder::internal::is_static_modint_t<mint>* = nullptr>
    std::vector<mint> convolution_large(const std::vector<mint>& a, const std::vector<mint>& b) {
        static constexpr int max_size = (mint::mod() - 1) & -(mint::mod() - 1);
        static constexpr int half_size = max_size >> 1;
        static constexpr int inv_max_size = atcoder::internal::inv_gcd(max_size, mint::mod()).second;

        const int n = int(a.size()), m = int(b.size());
        if (n + m - 1 <= max_size) return atcoder::convolution(a, b);
        if (n == 0 or m == 0) return {};
        if (std::min(n, m) <= 60) return internal::convolution_naive(a, b);

        const int dn = (n + half_size - 1) / half_size;
        const int dm = (m + half_size - 1) / half_size;

        std::vector<std::vector<mint>> as(dn), bs(dm);
        for (int i = 0; i < dn; ++i) {
            const int offset = half_size * i;
            as[i] = std::vector<mint>(a.begin() + offset, a.begin() + std::min(n, offset + half_size));
            as[i].resize(max_size);
            atcoder::internal::butterfly(as[i]);
        }
        for (int j = 0; j < dm; ++j) {
            const int offset = half_size * j;
            bs[j] = std::vector<mint>(b.begin() + offset, b.begin() + std::min(m, offset + half_size));
            bs[j].resize(max_size);
            atcoder::internal::butterfly(bs[j]);
        }
        std::vector<std::vector<mint>> cs(dn + dm - 1, std::vector<mint>(max_size));
        for (int i = 0; i < dn; ++i) for (int j = 0; j < dm; ++j) {
            for (int k = 0; k < max_size; ++k) cs[i + j][k] += as[i][k] * bs[j][k];
        }
        std::vector<mint> res(n + m - 1);
        for (int i = 0; i < dn + dm - 1; ++i) {
            atcoder::internal::butterfly_inv(cs[i]);
            const int offset = half_size * i;
            const int jmax = std::min(n + m - 1 - offset, max_size);
            for (int j = 0; j < jmax; ++j) {
                res[offset + j] += cs[i][j] * mint::raw(inv_max_size);
            }
        }
        return res;
    }
} // namespace suisen


#endif // SUISEN_CONVOLTION_LARGE

#line 1 "library/convolution/convolution_large.hpp"



#include <atcoder/convolution>

#line 1 "library/convolution/convolution_naive.hpp"



#include <vector>

namespace suisen::internal {
    template <typename T, typename R = T>
    std::vector<R> convolution_naive(const std::vector<T>& a, const std::vector<T>& b) {
        const int n = a.size(), m = b.size();
        std::vector<R> c(n + m - 1);
        if (n < m) {
            for (int j = 0; j < m; j++) for (int i = 0; i < n; i++) c[i + j] += R(a[i]) * b[j];
        } else {
            for (int i = 0; i < n; i++) for (int j = 0; j < m; j++) c[i + j] += R(a[i]) * b[j];
        }
        return c;
    }
} // namespace suisen



#line 7 "library/convolution/convolution_large.hpp"

namespace suisen {
    template <typename mint, atcoder::internal::is_static_modint_t<mint>* = nullptr>
    std::vector<mint> convolution_large(const std::vector<mint>& a, const std::vector<mint>& b) {
        static constexpr int max_size = (mint::mod() - 1) & -(mint::mod() - 1);
        static constexpr int half_size = max_size >> 1;
        static constexpr int inv_max_size = atcoder::internal::inv_gcd(max_size, mint::mod()).second;

        const int n = int(a.size()), m = int(b.size());
        if (n + m - 1 <= max_size) return atcoder::convolution(a, b);
        if (n == 0 or m == 0) return {};
        if (std::min(n, m) <= 60) return internal::convolution_naive(a, b);

        const int dn = (n + half_size - 1) / half_size;
        const int dm = (m + half_size - 1) / half_size;

        std::vector<std::vector<mint>> as(dn), bs(dm);
        for (int i = 0; i < dn; ++i) {
            const int offset = half_size * i;
            as[i] = std::vector<mint>(a.begin() + offset, a.begin() + std::min(n, offset + half_size));
            as[i].resize(max_size);
            atcoder::internal::butterfly(as[i]);
        }
        for (int j = 0; j < dm; ++j) {
            const int offset = half_size * j;
            bs[j] = std::vector<mint>(b.begin() + offset, b.begin() + std::min(m, offset + half_size));
            bs[j].resize(max_size);
            atcoder::internal::butterfly(bs[j]);
        }
        std::vector<std::vector<mint>> cs(dn + dm - 1, std::vector<mint>(max_size));
        for (int i = 0; i < dn; ++i) for (int j = 0; j < dm; ++j) {
            for (int k = 0; k < max_size; ++k) cs[i + j][k] += as[i][k] * bs[j][k];
        }
        std::vector<mint> res(n + m - 1);
        for (int i = 0; i < dn + dm - 1; ++i) {
            atcoder::internal::butterfly_inv(cs[i]);
            const int offset = half_size * i;
            const int jmax = std::min(n + m - 1 - offset, max_size);
            for (int j = 0; j < jmax; ++j) {
                res[offset + j] += cs[i][j] * mint::raw(inv_max_size);
            }
        }
        return res;
    }
} // namespace suisen

Back to top page