convolution
(src/convolution/convolution.hpp)

• View this file on GitHub
• Last update: 2024-11-09 02:16:49+09:00
• Include: #include "src/convolution/convolution.hpp"

convolution

vector<mint> convolution(vector<mint> a, vector<mint> b)

畳み込みを行います．
数列 $a_0, a_1, \cdots, a_{N - 1}$ と数列 $b_0, b_1, \cdots, b_{M - 1}$ から，長さ $N + M - 1$ の数列

\[c_i = \sum_{j = 0}^i a_j b_{i - j}\]

を計算します．

$a$ または $b$ が空配列の場合は空配列を返します．

制約

mint の法を $m$ として，

• $\mathrm{m}$ は素数
• $(\mathrm{m} - 1) \% 2^c = 0$ かつ $N + M - 1 \leq 2^c$ なる $c$ が存在する

計算量

$n = N + M$ として，

• $O(n \log n + \log m)$

butterfly

void butterfly(vector<mint>& a)

a のNumber-Theoretic Transformを計算します．
参照渡しであることに注意してください．

制約

• mint の法 $m$ が NTT-friendlyである

計算量

• $O(N \log N + \log m)$

butterfly_inv

void butterfly_inv(vector<mint>& a)

a のInverse Number-Theoretic Transformを計算します．
参照渡しであることに注意してください．

制約

• mint の法 $m$ が NTT-friendlyである．

計算量

• $O(N \log N + \log m)$

Depends on

• pow_mod (src/math/pow_mod.hpp)
• primitive_root (src/math/primitive_root.hpp)
• template (src/template/template.hpp)

Required by

• convolution_arbitrary (src/convolution/convolution_arbitrary.hpp)
• convolution_ll (src/convolution/convolution_ll.hpp)
• FormalPowerSeries (src/fps/formal_power_series.hpp)
• FormalPowerSeriesArbitrary (src/fps/formal_power_series_arbitrary.hpp)
• FormalPowerSeriesLL (src/fps/formal_power_series_ll.hpp)
• wildcard_pattern_matching (src/string/wildcard_pattern_matching.hpp)

Verified with

• verify/library_checker/convolution/convolution.test.cpp
• verify/library_checker/convolution/convolution_mod_1000000007.test.cpp
• verify/library_checker/other/find_linear_recurrence.test.cpp
• verify/library_checker/other/kth_term_of_linearly_recurrent_sequence.test.cpp
• verify/library_checker/polynomial/division_of_polynomial.test.cpp
• verify/library_checker/polynomial/exp_of_formal_power_series.test.cpp
• verify/library_checker/polynomial/inv_of_formal_power_series.test.cpp
• verify/library_checker/polynomial/log_of_formal_power_series.test.cpp
• verify/library_checker/polynomial/polynomial_taylor_shift.test.cpp
• verify/library_checker/polynomial/pow_of_formal_power_series.test.cpp
• verify/library_checker/polynomial/product_of_polynomial_sequence.test.cpp
• verify/library_checker/string/wildcard_pattern_matching.test.cpp
• verify/library_checker/tree/frequency_table_of_tree_distance.test.cpp
• verify/unit_test/convolution/convolution_ll.test.cpp
• verify/yukicoder/1559.test.cpp

Code

#pragma once
#include "../template/template.hpp"
#include "../math/primitive_root.hpp"
constexpr int countr_zero(const unsigned int n) {
    int res = 0;
    while(!(n & (1 << res))) ++res;
    return res;
}
template <typename mint, int g = primitive_root(mint::mod())>
struct FFT_Info {
    static constexpr int rank2 = countr_zero(mint::mod() - 1);
    array<mint, rank2 + 1> root;
    array<mint, rank2 + 1> iroot;
    array<mint, max(0, rank2 - 2 + 1)> rate2;
    array<mint, max(0, rank2 - 2 + 1)> irate2;
    array<mint, max(0, rank2 - 3 + 1)> rate3;
    array<mint, max(0, rank2 - 3 + 1)> irate3;
    FFT_Info() {
        root[rank2] = mint(g).pow((mint::mod() - 1) >> rank2);
        iroot[rank2] = root[rank2].inv();
        for(int i = rank2 - 1; i >= 0; --i) {
            root[i] = root[i + 1] * root[i + 1];
            iroot[i] = iroot[i + 1] * iroot[i + 1];
        }
        {
            mint prod = 1, iprod = 1;
            for(int i = 0; i <= rank2 - 2; ++i) {
                rate2[i] = root[i + 2] * prod;
                irate2[i] = iroot[i + 2] * iprod;
                prod *= iroot[i + 2];
                iprod *= root[i + 2];
            }
        }
        {
            mint prod = 1, iprod = 1;
            for(int i = 0; i <= rank2 - 3; ++i) {
                rate3[i] = root[i + 3] * prod;
                irate3[i] = iroot[i + 3] * iprod;
                prod *= iroot[i + 3];
                iprod *= root[i + 3];
            }
        }
    }
};
template <typename mint>
void butterfly(vector<mint>& a) {
    const int n = (int)a.size();
    const int h = __builtin_ctz((unsigned int)n);
    static const FFT_Info<mint> info;
    int len = 0;
    while(len < h) {
        if(h - len == 1) {
            const int p = 1 << (h - len - 1);
            mint rot = 1;
            for(int s = 0; s < (1 << len); ++s) {
                const int offset = s << (h - len);
                for(int i = 0; i < p; ++i) {
                    const auto l = a[i + offset];
                    const auto r = a[i + offset + p] * rot;
                    a[i + offset] = l + r;
                    a[i + offset + p] = l - r;
                }
                if(s + 1 != (1 << len)) rot *= info.rate2[__builtin_ctz(~(unsigned int)(s))];
            }
            ++len;
        } else {
            const int p = 1 << (h - len - 2);
            mint rot = 1, imag = info.root[2];
            for(int s = 0; s < (1 << len); ++s) {
                const mint rot2 = rot * rot;
                const mint rot3 = rot2 * rot;
                const int offset = s << (h - len);
                for(int i = 0; i < p; ++i) {
                    const auto mod2 = 1ULL * mint::mod() * mint::mod();
                    const auto a0 = 1ULL * a[i + offset].val();
                    const auto a1 = 1ULL * a[i + offset + p].val() * rot.val();
                    const auto a2 = 1ULL * a[i + offset + 2 * p].val() * rot2.val();
                    const auto a3 = 1ULL * a[i + offset + 3 * p].val() * rot3.val();
                    const auto a1na3imag = 1ULL * mint(a1 + mod2 - a3).val() * imag.val();
                    const auto na2 = mod2 - a2;
                    a[i + offset] = a0 + a2 + a1 + a3;
                    a[i + offset + 1 * p] = a0 + a2 + (2 * mod2 - (a1 + a3));
                    a[i + offset + 2 * p] = a0 + na2 + a1na3imag;
                    a[i + offset + 3 * p] = a0 + na2 + (mod2 - a1na3imag);
                }
                if(s + 1 != (1 << len)) rot *= info.rate3[__builtin_ctz(~(unsigned int)(s))];
            }
            len += 2;
        }
    }
}
template <typename mint>
void butterfly_inv(vector<mint>& a) {
    const int n = (int)a.size();
    const int h = __builtin_ctz((unsigned int)n);
    static const FFT_Info<mint> info;
    int len = h;
    while(len) {
        if(len == 1) {
            const int p = 1 << (h - len);
            mint irot = 1;
            for(int s = 0; s < (1 << (len - 1)); ++s) {
                const int offset = s << (h - len + 1);
                for(int i = 0; i < p; ++i) {
                    const auto l = a[i + offset];
                    const auto r = a[i + offset + p];
                    a[i + offset] = l + r;
                    a[i + offset + p] = (unsigned long long)(mint::mod() + l.val() - r.val()) * irot.val();
                }
                if(s + 1 != (1 << (len - 1))) irot *= info.irate2[__builtin_ctz(~(unsigned int)(s))];
            }
            --len;
        } else {
            const int p = 1 << (h - len);
            mint irot = 1, iimag = info.iroot[2];
            for(int s = 0; s < (1 << (len - 2)); ++s) {
                const mint irot2 = irot * irot;
                const mint irot3 = irot2 * irot;
                const int offset = s << (h - len + 2);
                for(int i = 0; i < p; ++i) {
                    const auto a0 = 1ULL * a[i + offset + 0 * p].val();
                    const auto a1 = 1ULL * a[i + offset + 1 * p].val();
                    const auto a2 = 1ULL * a[i + offset + 2 * p].val();
                    const auto a3 = 1ULL * a[i + offset + 3 * p].val();
                    const auto a2na3iimag = 1ULL * mint((mint::mod() + a2 - a3) * iimag.val()).val();
                    a[i + offset] = a0 + a1 + a2 + a3;
                    a[i + offset + 1 * p] = (a0 + (mint::mod() - a1) + a2na3iimag) * irot.val();
                    a[i + offset + 2 * p] = (a0 + a1 + (mint::mod() - a2) + (mint::mod() - a3)) * irot2.val();
                    a[i + offset + 3 * p] = (a0 + (mint::mod() - a1) + (mint::mod() - a2na3iimag)) * irot3.val();
                }
                if(s + 1 != (1 << (len - 2))) irot *= info.irate3[__builtin_ctz(~(unsigned int)(s))];
            }
            len -= 2;
        }
    }
}
template <typename mint>
vector<mint> convolution_naive(const vector<mint>& a, const vector<mint>& b) {
    const int n = (int)a.size(), m = (int)b.size();
    vector<mint> res(n + m - 1);
    if(n < m) {
        for(int j = 0; j < m; ++j) {
            for(int i = 0; i < n; ++i) {
                res[i + j] += a[i] * b[j];
            }
        }
    } else {
        for(int i = 0; i < n; ++i) {
            for(int j = 0; j < m; ++j) {
                res[i + j] += a[i] * b[j];
            }
        }
    }
    return res;
}
template <typename mint>
vector<mint> convolution(vector<mint> a, vector<mint> b) {
    const int n = (int)a.size(), m = (int)b.size();
    if(n == 0 or m == 0) return {};
    int z = 1;
    while(z < n + m - 1) z *= 2;
    assert((mint::mod() - 1) % z == 0);
    if(min(n, m) <= 60) return convolution_naive(a, b);
    a.resize(z);
    b.resize(z);
    butterfly(a);
    butterfly(b);
    for(int i = 0; i < z; ++i) a[i] *= b[i];
    butterfly_inv(a);
    a.resize(n + m - 1);
    const mint iz = mint(z).inv();
    for(int i = 0; i < n + m - 1; ++i) a[i] *= iz;
    return a;
}

#line 2 "src/template/template.hpp"
#include <bits/stdc++.h>
using namespace std;
using ll = long long;
using P = pair<long long, long long>;
#define rep(i, a, b) for(long long i = (a); i < (b); ++i)
#define rrep(i, a, b) for(long long i = (a); i >= (b); --i)
constexpr long long inf = 4e18;
struct SetupIO {
    SetupIO() {
        ios::sync_with_stdio(0);
        cin.tie(0);
        cout << fixed << setprecision(30);
    }
} setup_io;
#line 3 "src/math/pow_mod.hpp"
constexpr long long pow_mod(long long x, long long n, const long long mod) {
    assert(n >= 0 and mod >= 1);
    x %= mod;
    if(x < 0) x += mod;
    long long res = 1;
    while(n > 0) {
        if(n & 1) res = res * x % mod;
        x = x * x % mod;
        n >>= 1;
    }
    return res;
}
#line 4 "src/math/primitive_root.hpp"
constexpr int primitive_root(const int m) {
    if(m == 2) return 1;
    if(m == 167772161) return 3;
    if(m == 469762049) return 3;
    if(m == 754974721) return 11;
    if(m == 998244353) return 3;
    int divs[20] = {};
    divs[0] = 2;
    int cnt = 1;
    int x = (m - 1) / 2;
    while(x % 2 == 0) x /= 2;
    for(int i = 3; (long long)(i)*i <= x; i += 2) {
        if(x % i == 0) {
            divs[cnt++] = i;
            while(x % i == 0) {
                x /= i;
            }
        }
    }
    if(x > 1) {
        divs[cnt++] = x;
    }
    for(int g = 2;; ++g) {
        bool ok = true;
        for(int i = 0; i < cnt; ++i) {
            if(pow_mod(g, (m - 1) / divs[i], m) == 1) {
                ok = false;
                break;
            }
        }
        if(ok) return g;
    }
}
#line 4 "src/convolution/convolution.hpp"
constexpr int countr_zero(const unsigned int n) {
    int res = 0;
    while(!(n & (1 << res))) ++res;
    return res;
}
template <typename mint, int g = primitive_root(mint::mod())>
struct FFT_Info {
    static constexpr int rank2 = countr_zero(mint::mod() - 1);
    array<mint, rank2 + 1> root;
    array<mint, rank2 + 1> iroot;
    array<mint, max(0, rank2 - 2 + 1)> rate2;
    array<mint, max(0, rank2 - 2 + 1)> irate2;
    array<mint, max(0, rank2 - 3 + 1)> rate3;
    array<mint, max(0, rank2 - 3 + 1)> irate3;
    FFT_Info() {
        root[rank2] = mint(g).pow((mint::mod() - 1) >> rank2);
        iroot[rank2] = root[rank2].inv();
        for(int i = rank2 - 1; i >= 0; --i) {
            root[i] = root[i + 1] * root[i + 1];
            iroot[i] = iroot[i + 1] * iroot[i + 1];
        }
        {
            mint prod = 1, iprod = 1;
            for(int i = 0; i <= rank2 - 2; ++i) {
                rate2[i] = root[i + 2] * prod;
                irate2[i] = iroot[i + 2] * iprod;
                prod *= iroot[i + 2];
                iprod *= root[i + 2];
            }
        }
        {
            mint prod = 1, iprod = 1;
            for(int i = 0; i <= rank2 - 3; ++i) {
                rate3[i] = root[i + 3] * prod;
                irate3[i] = iroot[i + 3] * iprod;
                prod *= iroot[i + 3];
                iprod *= root[i + 3];
            }
        }
    }
};
template <typename mint>
void butterfly(vector<mint>& a) {
    const int n = (int)a.size();
    const int h = __builtin_ctz((unsigned int)n);
    static const FFT_Info<mint> info;
    int len = 0;
    while(len < h) {
        if(h - len == 1) {
            const int p = 1 << (h - len - 1);
            mint rot = 1;
            for(int s = 0; s < (1 << len); ++s) {
                const int offset = s << (h - len);
                for(int i = 0; i < p; ++i) {
                    const auto l = a[i + offset];
                    const auto r = a[i + offset + p] * rot;
                    a[i + offset] = l + r;
                    a[i + offset + p] = l - r;
                }
                if(s + 1 != (1 << len)) rot *= info.rate2[__builtin_ctz(~(unsigned int)(s))];
            }
            ++len;
        } else {
            const int p = 1 << (h - len - 2);
            mint rot = 1, imag = info.root[2];
            for(int s = 0; s < (1 << len); ++s) {
                const mint rot2 = rot * rot;
                const mint rot3 = rot2 * rot;
                const int offset = s << (h - len);
                for(int i = 0; i < p; ++i) {
                    const auto mod2 = 1ULL * mint::mod() * mint::mod();
                    const auto a0 = 1ULL * a[i + offset].val();
                    const auto a1 = 1ULL * a[i + offset + p].val() * rot.val();
                    const auto a2 = 1ULL * a[i + offset + 2 * p].val() * rot2.val();
                    const auto a3 = 1ULL * a[i + offset + 3 * p].val() * rot3.val();
                    const auto a1na3imag = 1ULL * mint(a1 + mod2 - a3).val() * imag.val();
                    const auto na2 = mod2 - a2;
                    a[i + offset] = a0 + a2 + a1 + a3;
                    a[i + offset + 1 * p] = a0 + a2 + (2 * mod2 - (a1 + a3));
                    a[i + offset + 2 * p] = a0 + na2 + a1na3imag;
                    a[i + offset + 3 * p] = a0 + na2 + (mod2 - a1na3imag);
                }
                if(s + 1 != (1 << len)) rot *= info.rate3[__builtin_ctz(~(unsigned int)(s))];
            }
            len += 2;
        }
    }
}
template <typename mint>
void butterfly_inv(vector<mint>& a) {
    const int n = (int)a.size();
    const int h = __builtin_ctz((unsigned int)n);
    static const FFT_Info<mint> info;
    int len = h;
    while(len) {
        if(len == 1) {
            const int p = 1 << (h - len);
            mint irot = 1;
            for(int s = 0; s < (1 << (len - 1)); ++s) {
                const int offset = s << (h - len + 1);
                for(int i = 0; i < p; ++i) {
                    const auto l = a[i + offset];
                    const auto r = a[i + offset + p];
                    a[i + offset] = l + r;
                    a[i + offset + p] = (unsigned long long)(mint::mod() + l.val() - r.val()) * irot.val();
                }
                if(s + 1 != (1 << (len - 1))) irot *= info.irate2[__builtin_ctz(~(unsigned int)(s))];
            }
            --len;
        } else {
            const int p = 1 << (h - len);
            mint irot = 1, iimag = info.iroot[2];
            for(int s = 0; s < (1 << (len - 2)); ++s) {
                const mint irot2 = irot * irot;
                const mint irot3 = irot2 * irot;
                const int offset = s << (h - len + 2);
                for(int i = 0; i < p; ++i) {
                    const auto a0 = 1ULL * a[i + offset + 0 * p].val();
                    const auto a1 = 1ULL * a[i + offset + 1 * p].val();
                    const auto a2 = 1ULL * a[i + offset + 2 * p].val();
                    const auto a3 = 1ULL * a[i + offset + 3 * p].val();
                    const auto a2na3iimag = 1ULL * mint((mint::mod() + a2 - a3) * iimag.val()).val();
                    a[i + offset] = a0 + a1 + a2 + a3;
                    a[i + offset + 1 * p] = (a0 + (mint::mod() - a1) + a2na3iimag) * irot.val();
                    a[i + offset + 2 * p] = (a0 + a1 + (mint::mod() - a2) + (mint::mod() - a3)) * irot2.val();
                    a[i + offset + 3 * p] = (a0 + (mint::mod() - a1) + (mint::mod() - a2na3iimag)) * irot3.val();
                }
                if(s + 1 != (1 << (len - 2))) irot *= info.irate3[__builtin_ctz(~(unsigned int)(s))];
            }
            len -= 2;
        }
    }
}
template <typename mint>
vector<mint> convolution_naive(const vector<mint>& a, const vector<mint>& b) {
    const int n = (int)a.size(), m = (int)b.size();
    vector<mint> res(n + m - 1);
    if(n < m) {
        for(int j = 0; j < m; ++j) {
            for(int i = 0; i < n; ++i) {
                res[i + j] += a[i] * b[j];
            }
        }
    } else {
        for(int i = 0; i < n; ++i) {
            for(int j = 0; j < m; ++j) {
                res[i + j] += a[i] * b[j];
            }
        }
    }
    return res;
}
template <typename mint>
vector<mint> convolution(vector<mint> a, vector<mint> b) {
    const int n = (int)a.size(), m = (int)b.size();
    if(n == 0 or m == 0) return {};
    int z = 1;
    while(z < n + m - 1) z *= 2;
    assert((mint::mod() - 1) % z == 0);
    if(min(n, m) <= 60) return convolution_naive(a, b);
    a.resize(z);
    b.resize(z);
    butterfly(a);
    butterfly(b);
    for(int i = 0; i < z; ++i) a[i] *= b[i];
    butterfly_inv(a);
    a.resize(n + m - 1);
    const mint iz = mint(z).inv();
    for(int i = 0; i < n + m - 1; ++i) a[i] *= iz;
    return a;
}

Back to top page