CP-Library
This documentation is automatically generated by online-judge-tools/verification-helper
math/convolution/Convolution.hpp
- View this file on GitHub
- Last update: 2026-08-24 14:55:40+07:00
- Include:
#include "math/convolution/Convolution.hpp" - Link: https://judge.yosupo.jp/problem/convolution_mod_2_64
Depends on
Verified with
- test/1_library_checker/math/convolution/Convolution107.test.cpp
- test/1_library_checker/math/convolution/Convolution2^64.test.cpp
- test/1_library_checker/math/convolution/Convolution998.test.cpp
- test/1_library_checker/math/convolution/ConvolutionLarge.test.cpp
Code
#pragma once
#include <vector>
#include <cstdint>
#include <type_traits>
#include <algorithm>
#include <cassert>
#include "math/ModInt.hpp"
namespace Convolution {
constexpr uint64_t primitive_root_for_convolution(uint64_t p) {
if (p == 2) return 1;
if (p == 998244353) return 3;
if (p == 469762049) return 3;
if (p == 1811939329) return 11;
if (p == 2013265921) return 11;
for (uint64_t g = 2; g < p; ++g) {
uint64_t base = g, exp = (p - 1) >> 1, res = 1;
while (exp > 0) {
if (exp & 1) res = (res * base) % p;
base = (base * base) % p;
exp >>= 1;
}
if (res != 1) return g;
}
return -1;
}
namespace internal {
constexpr int msb(uint32_t x) { return x == 0 ? -1 : 31 - __builtin_clz(x); }
constexpr int ceil_log2(uint32_t x) { return x <= 1 ? 0 : 32 - __builtin_clz(x - 1); }
template <class T>
class NthRoot {
private:
static constexpr unsigned int lg = msb((T::get_mod() - 1) & (1 - T::get_mod()));
T root[lg + 1];
T inv_root[lg + 1];
T rate[lg + 1];
T inv_rate[lg + 1];
public:
constexpr NthRoot() : root{}, inv_root{}, rate{}, inv_rate{} {
root[lg] = T(primitive_root_for_convolution(T::get_mod())).pow((T::get_mod() - 1) >> lg);
inv_root[lg] = root[lg].inv();
for (int i = lg - 1; i >= 0; i--) {
root[i] = root[i + 1] * root[i + 1];
inv_root[i] = inv_root[i + 1] * inv_root[i + 1];
}
T r = 1;
for (int i = 2; i <= (int)lg; i++) {
rate[i - 2] = r * root[i];
r = r * inv_root[i];
}
r = 1;
for (int i = 2; i <= (int)lg; i++) {
inv_rate[i - 2] = r * inv_root[i];
r = r * root[i];
}
}
static constexpr unsigned int get_lg() { return lg; }
constexpr T get(int n) const { return root[n]; }
constexpr T inv(int n) const { return inv_root[n]; }
constexpr T get_rate(int n) const { return rate[n]; }
constexpr T get_inv_rate(int n) const { return inv_rate[n]; }
};
template <class T>
void number_theoretic_transform(std::vector<T>& a) {
static constexpr NthRoot<T> nth_root;
int n = (int)a.size();
for (int i = n >> 1; i > 0; i >>= 1) {
T z = T::raw(1);
for (int j = 0; j < n; j += (i << 1)) {
for (int k = 0; k < i; k++) {
const T x = a[j + k];
const T y = a[j + i + k] * z;
a[j + k] = x + y;
a[j + i + k] = x - y;
}
z *= nth_root.get_rate(__builtin_popcount(j & ~(j + (i << 1))));
}
}
}
template <class T>
void inverse_number_theoretic_transform(std::vector<T>& a) {
static constexpr NthRoot<T> nth_root;
int n = (int)a.size();
for (int i = 1; i < n; i <<= 1) {
T z = T::raw(1);
for (int j = 0; j < n; j += (i << 1)) {
for (int k = 0; k < i; k++) {
const T x = a[j + k];
const T y = a[j + i + k];
a[j + k] = x + y;
a[j + i + k] = (x - y) * z;
}
z *= nth_root.get_inv_rate(__builtin_popcount(j & ~(j + (i << 1))));
}
}
T inv_n = T(1) / T(n);
for (auto& x : a) x *= inv_n;
}
template <class T>
std::vector<T> convolution_naive(const std::vector<T>& a, const std::vector<T>& b) {
int n = (int)a.size(), m = (int)b.size();
std::vector<T> c(n + m - 1);
for (int i = 0; i < n; i++) {
for (int j = 0; j < m; j++) c[i + j] += a[i] * b[j];
}
return c;
}
template <class T>
std::vector<T> convolution_pow2(std::vector<T> a) {
int n = (int)a.size() * 2 - 1;
if (n <= 1) return convolution_naive(a, a);
int lg = msb(n - 1) + 1;
if (n - (1 << (lg - 1)) <= 5) {
--lg;
int m = (int)a.size() - (1 << (lg - 1));
std::vector<T> a1(a.begin(), a.begin() + m);
std::vector<T> a2(a.begin() + m, a.end());
std::vector<T> c(n);
std::vector<T> c1 = convolution_naive(a1, a1);
std::vector<T> c2 = convolution_naive(a1, a2);
std::vector<T> c3 = convolution_pow2(a2);
for (int i = 0; i < (int)c1.size(); i++) c[i] += c1[i];
for (int i = 0; i < (int)c2.size(); i++) c[i + m] += c2[i] * 2;
for (int i = 0; i < (int)c3.size(); i++) c[i + m * 2] += c3[i];
return c;
}
int m = 1 << lg;
a.resize(m);
number_theoretic_transform(a);
for (int i = 0; i < m; i++) a[i] *= a[i];
inverse_number_theoretic_transform(a);
a.resize(n);
return a;
}
template <class T>
std::vector<T> convolution(std::vector<T> a, std::vector<T> b) {
int n = (int)a.size() + (int)b.size() - 1;
if (n <= 1) return convolution_naive(a, b);
int lg = ceil_log2(n);
int m = 1 << lg;
if (n - (1 << (lg - 1)) <= 5) {
--lg;
if (a.size() < b.size()) std::swap(a, b);
int half_m = n - (1 << lg);
std::vector<T> a1(a.begin(), a.begin() + half_m);
std::vector<T> a2(a.begin() + half_m, a.end());
std::vector<T> c(n);
std::vector<T> c1 = convolution_naive(a1, b);
std::vector<T> c2 = convolution(a2, b);
for (int i = 0; i < (int)c1.size(); i++) c[i] += c1[i];
for (int i = 0; i < (int)c2.size(); i++) c[i + half_m] += c2[i];
return c;
}
a.resize(m);
b.resize(m);
number_theoretic_transform(a);
number_theoretic_transform(b);
for (int i = 0; i < m; i++) a[i] *= b[i];
inverse_number_theoretic_transform(a);
a.resize(n);
return a;
}
} // namespace internal
template <uint32_t p>
struct is_ntt_friendly : std::false_type {};
template <> struct is_ntt_friendly<998244353> : std::true_type {};
template <> struct is_ntt_friendly<469762049> : std::true_type {};
template <> struct is_ntt_friendly<1811939329> : std::true_type {};
template <> struct is_ntt_friendly<2013265921> : std::true_type {};
template <uint32_t p>
std::vector<StaticModInt<p>> convolution_for_any_mod(const std::vector<StaticModInt<p>>& a, const std::vector<StaticModInt<p>>& b);
template <uint32_t p>
std::vector<StaticModInt<p>> convolution_large(const std::vector<StaticModInt<p>>& a, const std::vector<StaticModInt<p>>& b);
template <uint32_t p>
std::vector<StaticModInt<p>> convolution(const std::vector<StaticModInt<p>>& a, const std::vector<StaticModInt<p>>& b) {
int n = (int)a.size(), m = (int)b.size();
if (n == 0 || m == 0) return {};
if (n <= 60 || m <= 60) return internal::convolution_naive(a, b);
if (n + m - 1 <= (int)((1 - p) & (p - 1))) {
if (n == m && a == b) return internal::convolution_pow2(a);
return internal::convolution(a, b);
}
if constexpr (is_ntt_friendly<p>::value) {
return convolution_large<p>(a, b);
}
return convolution_for_any_mod(a, b);
}
template <uint32_t p>
std::vector<StaticModInt<p>> convolution_for_any_mod(const std::vector<StaticModInt<p>>& a, const std::vector<StaticModInt<p>>& b) {
int n = (int)a.size(), m = (int)b.size();
assert(n + m - 1 <= (1 << 26));
static constexpr uint32_t MOD1 = 469762049;
static constexpr uint32_t MOD2 = 1811939329;
static constexpr uint32_t MOD3 = 2013265921;
using mint1 = StaticModInt<MOD1>;
using mint2 = StaticModInt<MOD2>;
using mint3 = StaticModInt<MOD3>;
std::vector<mint1> a1(n), b1(m);
std::vector<mint2> a2(n), b2(m);
std::vector<mint3> a3(n), b3(m);
for (int i = 0; i < n; i++) {
uint32_t v = a[i].get();
a1[i] = v; a2[i] = v; a3[i] = v;
}
for (int i = 0; i < m; i++) {
uint32_t v = b[i].get();
b1[i] = v; b2[i] = v; b3[i] = v;
}
auto c1 = internal::convolution(a1, b1);
auto c2 = internal::convolution(a2, b2);
auto c3 = internal::convolution(a3, b3);
static constexpr uint32_t INV1_2 = mint2(MOD1).inv().get();
static constexpr uint32_t INV1_3 = mint3(MOD1).inv().get();
static constexpr uint32_t INV2_3 = mint3(MOD2).inv().get();
std::vector<StaticModInt<p>> res(n + m - 1);
for (int i = 0; i < n + m - 1; i++) {
int64_t t1 = c1[i].get();
int64_t t2 = (c2[i].get() - t1 + MOD2) * INV1_2 % MOD2;
if (t2 < 0) t2 += MOD2;
int64_t t3 = ((c3[i].get() - t1 + MOD3) * INV1_3 % MOD3 - t2 + MOD3) * INV2_3 % MOD3;
if (t3 < 0) t3 += MOD3;
res[i] = StaticModInt<p>(t1 + (t2 + t3 * MOD2) % p * (MOD1 % p));
}
return res;
}
template <uint32_t p>
std::vector<StaticModInt<p>> convolution_large(const std::vector<StaticModInt<p>>& a, const std::vector<StaticModInt<p>>& b) {
int n = (int)a.size(), m = (int)b.size();
int limit = (int)((1 - p) & (p - 1));
int L = limit >> 1;
int k = (n + L - 1) / L;
int l = (m + L - 1) / L;
std::vector<std::vector<StaticModInt<p>>> A_hat(k, std::vector<StaticModInt<p>>(limit));
std::vector<std::vector<StaticModInt<p>>> B_hat(l, std::vector<StaticModInt<p>>(limit));
for (int i = 0; i < k; ++i) {
int start = i * L;
int end = std::min(n, start + L);
for (int j = start; j < end; ++j) {
A_hat[i][j - start] = a[j];
}
internal::number_theoretic_transform(A_hat[i]);
}
for (int i = 0; i < l; ++i) {
int start = i * L;
int end = std::min(m, start + L);
for (int j = start; j < end; ++j) {
B_hat[i][j - start] = b[j];
}
internal::number_theoretic_transform(B_hat[i]);
}
std::vector<StaticModInt<p>> res(n + m - 1);
for (int c = 0; c < k + l - 1; ++c) {
std::vector<StaticModInt<p>> C_chunk(limit);
for (int i = 0; i < k; ++i) {
int j = c - i;
if (0 <= j && j < l) {
for (int x = 0; x < limit; ++x) {
C_chunk[x] += A_hat[i][x] * B_hat[j][x];
}
}
}
internal::inverse_number_theoretic_transform(C_chunk);
int start = c * L;
int end = std::min(n + m - 1, start + limit);
for (int j = start; j < end; ++j) {
res[j] += C_chunk[j - start];
}
}
return res;
}
template <uint32_t p>
std::vector<int64_t> convolution(const std::vector<int64_t>& a, const std::vector<int64_t>& b) {
int n = (int)a.size(), m = (int)b.size();
std::vector<StaticModInt<p>> a2(n), b2(m);
for (int i = 0; i < n; i++) a2[i] = a[i];
for (int i = 0; i < m; i++) b2[i] = b[i];
auto c2 = convolution<p>(a2, b2);
std::vector<int64_t> c(c2.size());
for (int i = 0; i < (int)c2.size(); i++) c[i] = c2[i].get();
return c;
}
inline std::vector<int64_t> convolution_ll(const std::vector<int64_t>& a, const std::vector<int64_t>& b) {
int n = (int)a.size(), m = (int)b.size();
if (n == 0 || m == 0) return {};
static constexpr uint32_t MOD1 = 469762049;
static constexpr uint32_t MOD2 = 1811939329;
static constexpr uint32_t MOD3 = 2013265921;
using mint1 = StaticModInt<MOD1>;
using mint2 = StaticModInt<MOD2>;
using mint3 = StaticModInt<MOD3>;
std::vector<mint1> a1(n), b1(m);
std::vector<mint2> a2(n), b2(m);
std::vector<mint3> a3(n), b3(m);
for (int i = 0; i < n; i++) {
int64_t v = a[i] % MOD1; if (v < 0) v += MOD1; a1[i] = v;
v = a[i] % MOD2; if (v < 0) v += MOD2; a2[i] = v;
v = a[i] % MOD3; if (v < 0) v += MOD3; a3[i] = v;
}
for (int i = 0; i < m; i++) {
int64_t v = b[i] % MOD1; if (v < 0) v += MOD1; b1[i] = v;
v = b[i] % MOD2; if (v < 0) v += MOD2; b2[i] = v;
v = b[i] % MOD3; if (v < 0) v += MOD3; b3[i] = v;
}
auto c1 = internal::convolution(a1, b1);
auto c2 = internal::convolution(a2, b2);
auto c3 = internal::convolution(a3, b3);
static constexpr uint32_t INV1_2 = mint2(MOD1).inv().get();
static constexpr uint32_t INV1_3 = mint3(MOD1).inv().get();
static constexpr uint32_t INV2_3 = mint3(MOD2).inv().get();
std::vector<int64_t> res(n + m - 1);
for (int i = 0; i < n + m - 1; i++) {
int64_t t1 = c1[i].get();
int64_t t2 = c2[i].get() - t1;
t2 %= MOD2;
if (t2 < 0) t2 += MOD2;
t2 = (t2 * INV1_2) % MOD2;
int64_t t3 = c3[i].get() - t1;
t3 %= MOD3;
if (t3 < 0) t3 += MOD3;
t3 = (t3 * INV1_3) % MOD3;
t3 = (t3 - t2) % MOD3;
if (t3 < 0) t3 += MOD3;
t3 = (t3 * INV2_3) % MOD3;
__int128_t result = t1;
result += (__int128_t)t2 * MOD1;
result += (__int128_t)t3 * MOD1 * MOD2;
res[i] = (int64_t)result;
}
return res;
}
// https://judge.yosupo.jp/problem/convolution_mod_2_64
inline std::vector<uint64_t> convolution_64(const std::vector<uint64_t>& a, const std::vector<uint64_t>& b) {
int n = (int)a.size(), m = (int)b.size();
if (n == 0 || m == 0) return {};
std::vector<int64_t> a_lo(n), a_hi(n), b_lo(m), b_hi(m);
for (int i = 0; i < n; i++) {
a_lo[i] = a[i] & 0xFFFFFFFF;
a_hi[i] = a[i] >> 32;
}
for (int i = 0; i < m; i++) {
b_lo[i] = b[i] & 0xFFFFFFFF;
b_hi[i] = b[i] >> 32;
}
auto c_lolo = convolution_ll(a_lo, b_lo);
auto c_hilo = convolution_ll(a_hi, b_lo);
auto c_lohi = convolution_ll(a_lo, b_hi);
std::vector<uint64_t> res(n + m - 1);
for (int i = 0; i < n + m - 1; i++) {
uint64_t lolo = (uint64_t)c_lolo[i];
uint64_t hilo = (uint64_t)c_hilo[i];
uint64_t lohi = (uint64_t)c_lohi[i];
res[i] = lolo + (hilo << 32) + (lohi << 32);
}
return res;
}
template <class T>
void ntt_doubling_(std::vector<T>& a, std::vector<T> b) {
static constexpr internal::NthRoot<T> nth_root;
int n = (int)a.size();
const T z = nth_root.get(internal::msb(n) + 1);
T r = 1;
for (int i = 0; i < n; i++) {
b[i] *= r;
r *= z;
}
internal::number_theoretic_transform(b);
a.reserve(2 * n);
a.insert(a.end(), b.begin(), b.end());
}
template <class T>
void ntt_doubling_(std::vector<T>& a) {
auto b = a;
internal::inverse_number_theoretic_transform(b);
ntt_doubling_(a, std::move(b));
}
template <int ID>
std::vector<DynamicModInt<ID>> convolution(const std::vector<DynamicModInt<ID>>& a, const std::vector<DynamicModInt<ID>>& b) {
int n = (int)a.size(), m = (int)b.size();
if (n == 0 || m == 0) return {};
if (n <= 60 || m <= 60) return internal::convolution_naive(a, b);
uint32_t mod = DynamicModInt<ID>::get_mod();
if (n + m - 1 <= (int)((1 - mod) & (mod - 1))) {
if (n == m && a == b) return internal::convolution_pow2(a);
return internal::convolution(a, b);
}
static constexpr uint32_t MOD1 = 469762049;
static constexpr uint32_t MOD2 = 1811939329;
static constexpr uint32_t MOD3 = 2013265921;
using mint1 = StaticModInt<MOD1>;
using mint2 = StaticModInt<MOD2>;
using mint3 = StaticModInt<MOD3>;
std::vector<mint1> a1(n), b1(m);
std::vector<mint2> a2(n), b2(m);
std::vector<mint3> a3(n), b3(m);
for (int i = 0; i < n; i++) {
uint32_t v = a[i].get();
a1[i] = v; a2[i] = v; a3[i] = v;
}
for (int i = 0; i < m; i++) {
uint32_t v = b[i].get();
b1[i] = v; b2[i] = v; b3[i] = v;
}
auto c1 = internal::convolution(a1, b1);
auto c2 = internal::convolution(a2, b2);
auto c3 = internal::convolution(a3, b3);
static constexpr uint32_t INV1_2 = mint2(MOD1).inv().get();
static constexpr uint32_t INV1_3 = mint3(MOD1).inv().get();
static constexpr uint32_t INV2_3 = mint3(MOD2).inv().get();
std::vector<DynamicModInt<ID>> res(n + m - 1);
for (int i = 0; i < n + m - 1; i++) {
int64_t t1 = c1[i].get();
int64_t t2 = (c2[i].get() - t1 + MOD2) * INV1_2 % MOD2;
if (t2 < 0) t2 += MOD2;
int64_t t3 = ((c3[i].get() - t1 + MOD3) * INV1_3 % MOD3 - t2 + MOD3) * INV2_3 % MOD3;
if (t3 < 0) t3 += MOD3;
res[i] = DynamicModInt<ID>(t1 + (t2 + t3 * MOD2) % mod * MOD1);
}
return res;
}
} // namespace Convolution#line 2 "math/convolution/Convolution.hpp"
#include <vector>
#include <cstdint>
#include <type_traits>
#include <algorithm>
#include <cassert>
#line 3 "math/ModInt.hpp"
#include <iostream>
#line 5 "math/ModInt.hpp"
#include <limits>
#line 3 "math/Barrett.hpp"
#include <utility>
// reference: https://maspypy.github.io/library/mod/barrett.hpp
struct Barrett {
uint32_t mod;
uint64_t inv_mod;
explicit Barrett(uint32_t m = 1) : mod(m), inv_mod((uint64_t)(-1) / m + 1) {}
uint32_t umod() const { return mod; }
uint32_t modulo(uint64_t z) const {
if (mod == 1) return 0;
uint64_t x = (uint64_t)(((__uint128_t)(z)*inv_mod) >> 64);
uint64_t y = x * mod;
return (uint32_t)(z - y + (z < y ? mod : 0));
}
uint64_t floor(uint64_t z) const {
if (mod == 1) return z;
uint64_t x = (uint64_t)(((__uint128_t)(z)*inv_mod) >> 64);
uint64_t y = x * mod;
return (z < y ? x - 1 : x);
}
std::pair<uint64_t, uint32_t> divmod(uint64_t z) const {
if (mod == 1) return {z, 0};
uint64_t x = (uint64_t)(((__uint128_t)(z)*inv_mod) >> 64);
uint64_t y = x * mod;
if (z < y) return {x - 1, z - y + mod};
return {x, z - y};
}
uint32_t mul(uint32_t a, uint32_t b) const { return modulo((uint64_t)(a)*b); }
};
#line 8 "math/ModInt.hpp"
template <uint32_t Mod>
class StaticModInt {
static_assert(Mod > 0, "Modulus must be positive");
static_assert(Mod <= std::numeric_limits<uint32_t>::max() / 2,
"Modulus * 2 must be <= uint32_t::max()");
private:
uint32_t val;
static constexpr uint32_t inv1000000007[] = {
0, 1, 500000004, 333333336, 250000002, 400000003,
166666668, 142857144, 125000001, 111111112, 700000005
};
static constexpr uint32_t inv998244353[] = {
0, 1, 499122177, 332748118, 748683265, 598946612,
166374059, 855638017, 873463809, 443664157, 299473306
};
static constexpr uint32_t mod_inv(int64_t a) {
int64_t b = Mod, x = 1, u = 0, t = 0, tmp = 0;
while (b) {
t = a / b;
tmp = a - t * b;
a = b;
b = tmp;
tmp = x - t * u;
x = u;
u = tmp;
}
if (x < 0) x += Mod;
return (uint32_t)x;
}
static constexpr uint32_t mod_reduce(int64_t v) {
int64_t x = v % (int64_t)Mod;
if (x < 0) x += Mod;
return (uint32_t)x;
}
public:
constexpr StaticModInt() : val(0) {}
template <class U, std::enable_if_t<std::is_integral_v<U> && std::is_signed_v<U>, int> = 0>
constexpr StaticModInt(U v) : val(mod_reduce((int64_t)v)) {}
template <class U, std::enable_if_t<std::is_integral_v<U> && std::is_unsigned_v<U>, int> = 0>
constexpr StaticModInt(U v) : val((uint32_t)(v % Mod)) {}
constexpr uint32_t get() const { return val; }
static constexpr uint32_t get_mod() { return Mod; }
static constexpr StaticModInt raw(uint32_t v) {
StaticModInt res;
res.val = v;
return res;
}
constexpr StaticModInt inv() const {
if constexpr (Mod == 1000000007) {
if (val <= 10) return raw(inv1000000007[val]);
} else if constexpr (Mod == 998244353) {
if (val <= 10) return raw(inv998244353[val]);
}
return raw(mod_inv(val));
}
constexpr StaticModInt& operator++() {
if (++val == Mod) val = 0;
return *this;
}
constexpr StaticModInt operator++(int) {
StaticModInt res = *this;
++*this;
return res;
}
constexpr StaticModInt& operator--() {
if (val == 0) val = Mod;
--val;
return *this;
}
constexpr StaticModInt operator--(int) {
StaticModInt res = *this;
--*this;
return res;
}
constexpr StaticModInt& operator+=(const StaticModInt& other) {
val += other.val;
if (val >= Mod) val -= Mod;
return *this;
}
constexpr StaticModInt& operator-=(const StaticModInt& other) {
if (val < other.val) val += Mod;
val -= other.val;
return *this;
}
constexpr StaticModInt& operator*=(const StaticModInt& other) {
val = (uint32_t)((uint64_t)val * other.val % Mod);
return *this;
}
constexpr StaticModInt& operator/=(const StaticModInt& other) {
*this *= other.inv();
return *this;
}
friend constexpr StaticModInt operator+(const StaticModInt& lhs, const StaticModInt& rhs) { return StaticModInt(lhs) += rhs; }
friend constexpr StaticModInt operator-(const StaticModInt& lhs, const StaticModInt& rhs) { return StaticModInt(lhs) -= rhs; }
friend constexpr StaticModInt operator*(const StaticModInt& lhs, const StaticModInt& rhs) { return StaticModInt(lhs) *= rhs; }
friend constexpr StaticModInt operator/(const StaticModInt& lhs, const StaticModInt& rhs) { return StaticModInt(lhs) /= rhs; }
constexpr StaticModInt operator+() const { return *this; }
constexpr StaticModInt operator-() const { return StaticModInt() - *this; }
friend constexpr bool operator==(const StaticModInt& lhs, const StaticModInt& rhs) { return lhs.val == rhs.val; }
friend constexpr bool operator!=(const StaticModInt& lhs, const StaticModInt& rhs) { return lhs.val != rhs.val; }
constexpr StaticModInt pow(int64_t a) const {
StaticModInt v = *this, res = 1;
while (a) {
if (a & 1) res *= v;
a >>= 1;
v *= v;
}
return res;
}
friend std::ostream& operator<<(std::ostream& os, const StaticModInt& m) { return os << m.val; }
friend std::istream& operator>>(std::istream& is, StaticModInt& m) {
int64_t v; is >> v;
m = StaticModInt(v);
return is;
}
};
using modint107 = StaticModInt<1000000007>;
using modint998 = StaticModInt<998244353>;
template <int ID>
class DynamicModInt {
private:
uint32_t val;
static Barrett bt;
static uint32_t mod_inv(int64_t a) {
int64_t b = bt.mod, x = 1, u = 0, t = 0, tmp = 0;
while (b) {
t = a / b;
tmp = a - t * b;
a = b;
b = tmp;
tmp = x - t * u;
x = u;
u = tmp;
}
if (x < 0) x += bt.mod;
return (uint32_t)x;
}
public:
DynamicModInt() : val(0) {}
template <class U, std::enable_if_t<std::is_integral_v<U> && std::is_signed_v<U>, int> = 0>
DynamicModInt(U v) {
int64_t x = (int64_t)(v % (int64_t)bt.mod);
if (x < 0) x += bt.mod;
val = (uint32_t)x;
}
template <class U, std::enable_if_t<std::is_integral_v<U> && std::is_unsigned_v<U>, int> = 0>
DynamicModInt(U v) : val((uint32_t)(v % bt.mod)) {}
uint32_t get() const { return val; }
static uint32_t get_mod() { return bt.mod; }
static void set_mod(uint32_t m) {
assert(m > 0 && m <= std::numeric_limits<uint32_t>::max() / 2);
bt = Barrett(m);
}
static DynamicModInt raw(uint32_t v) {
DynamicModInt res;
res.val = v;
return res;
}
DynamicModInt inv() const { return raw(mod_inv(val)); }
DynamicModInt& operator++() { if (++val == bt.mod) val = 0; return *this; }
DynamicModInt operator++(int) { DynamicModInt res = *this; ++*this; return res; }
DynamicModInt& operator--() { if (val == 0) val = bt.mod; --val; return *this; }
DynamicModInt operator--(int) { DynamicModInt res = *this; --*this; return res; }
DynamicModInt& operator+=(const DynamicModInt& other) { val += other.val; if (val >= bt.mod) val -= bt.mod; return *this; }
DynamicModInt& operator-=(const DynamicModInt& other) { if (val < other.val) val += bt.mod; val -= other.val; return *this; }
DynamicModInt& operator*=(const DynamicModInt& other) { val = bt.mul(val, other.val); return *this; }
DynamicModInt& operator/=(const DynamicModInt& other) { *this *= other.inv(); return *this; }
friend DynamicModInt operator+(const DynamicModInt& lhs, const DynamicModInt& rhs) { return DynamicModInt(lhs) += rhs; }
friend DynamicModInt operator-(const DynamicModInt& lhs, const DynamicModInt& rhs) { return DynamicModInt(lhs) -= rhs; }
friend DynamicModInt operator*(const DynamicModInt& lhs, const DynamicModInt& rhs) { return DynamicModInt(lhs) *= rhs; }
friend DynamicModInt operator/(const DynamicModInt& lhs, const DynamicModInt& rhs) { return DynamicModInt(lhs) /= rhs; }
DynamicModInt operator+() const { return *this; }
DynamicModInt operator-() const { return DynamicModInt() - *this; }
friend bool operator==(const DynamicModInt& lhs, const DynamicModInt& rhs) { return lhs.val == rhs.val; }
friend bool operator!=(const DynamicModInt& lhs, const DynamicModInt& rhs) { return lhs.val != rhs.val; }
DynamicModInt pow(int64_t a) const {
DynamicModInt v = *this, res = 1;
while (a) {
if (a & 1) res *= v;
a >>= 1;
v *= v;
}
return res;
}
friend std::ostream& operator<<(std::ostream& os, const DynamicModInt& m) { return os << m.val; }
friend std::istream& operator>>(std::istream& is, DynamicModInt& m) {
int64_t v; is >> v;
m = DynamicModInt(v);
return is;
}
};
template <int ID>
Barrett DynamicModInt<ID>::bt = Barrett(998244353);
using dmint = DynamicModInt<0>;
#line 8 "math/convolution/Convolution.hpp"
namespace Convolution {
constexpr uint64_t primitive_root_for_convolution(uint64_t p) {
if (p == 2) return 1;
if (p == 998244353) return 3;
if (p == 469762049) return 3;
if (p == 1811939329) return 11;
if (p == 2013265921) return 11;
for (uint64_t g = 2; g < p; ++g) {
uint64_t base = g, exp = (p - 1) >> 1, res = 1;
while (exp > 0) {
if (exp & 1) res = (res * base) % p;
base = (base * base) % p;
exp >>= 1;
}
if (res != 1) return g;
}
return -1;
}
namespace internal {
constexpr int msb(uint32_t x) { return x == 0 ? -1 : 31 - __builtin_clz(x); }
constexpr int ceil_log2(uint32_t x) { return x <= 1 ? 0 : 32 - __builtin_clz(x - 1); }
template <class T>
class NthRoot {
private:
static constexpr unsigned int lg = msb((T::get_mod() - 1) & (1 - T::get_mod()));
T root[lg + 1];
T inv_root[lg + 1];
T rate[lg + 1];
T inv_rate[lg + 1];
public:
constexpr NthRoot() : root{}, inv_root{}, rate{}, inv_rate{} {
root[lg] = T(primitive_root_for_convolution(T::get_mod())).pow((T::get_mod() - 1) >> lg);
inv_root[lg] = root[lg].inv();
for (int i = lg - 1; i >= 0; i--) {
root[i] = root[i + 1] * root[i + 1];
inv_root[i] = inv_root[i + 1] * inv_root[i + 1];
}
T r = 1;
for (int i = 2; i <= (int)lg; i++) {
rate[i - 2] = r * root[i];
r = r * inv_root[i];
}
r = 1;
for (int i = 2; i <= (int)lg; i++) {
inv_rate[i - 2] = r * inv_root[i];
r = r * root[i];
}
}
static constexpr unsigned int get_lg() { return lg; }
constexpr T get(int n) const { return root[n]; }
constexpr T inv(int n) const { return inv_root[n]; }
constexpr T get_rate(int n) const { return rate[n]; }
constexpr T get_inv_rate(int n) const { return inv_rate[n]; }
};
template <class T>
void number_theoretic_transform(std::vector<T>& a) {
static constexpr NthRoot<T> nth_root;
int n = (int)a.size();
for (int i = n >> 1; i > 0; i >>= 1) {
T z = T::raw(1);
for (int j = 0; j < n; j += (i << 1)) {
for (int k = 0; k < i; k++) {
const T x = a[j + k];
const T y = a[j + i + k] * z;
a[j + k] = x + y;
a[j + i + k] = x - y;
}
z *= nth_root.get_rate(__builtin_popcount(j & ~(j + (i << 1))));
}
}
}
template <class T>
void inverse_number_theoretic_transform(std::vector<T>& a) {
static constexpr NthRoot<T> nth_root;
int n = (int)a.size();
for (int i = 1; i < n; i <<= 1) {
T z = T::raw(1);
for (int j = 0; j < n; j += (i << 1)) {
for (int k = 0; k < i; k++) {
const T x = a[j + k];
const T y = a[j + i + k];
a[j + k] = x + y;
a[j + i + k] = (x - y) * z;
}
z *= nth_root.get_inv_rate(__builtin_popcount(j & ~(j + (i << 1))));
}
}
T inv_n = T(1) / T(n);
for (auto& x : a) x *= inv_n;
}
template <class T>
std::vector<T> convolution_naive(const std::vector<T>& a, const std::vector<T>& b) {
int n = (int)a.size(), m = (int)b.size();
std::vector<T> c(n + m - 1);
for (int i = 0; i < n; i++) {
for (int j = 0; j < m; j++) c[i + j] += a[i] * b[j];
}
return c;
}
template <class T>
std::vector<T> convolution_pow2(std::vector<T> a) {
int n = (int)a.size() * 2 - 1;
if (n <= 1) return convolution_naive(a, a);
int lg = msb(n - 1) + 1;
if (n - (1 << (lg - 1)) <= 5) {
--lg;
int m = (int)a.size() - (1 << (lg - 1));
std::vector<T> a1(a.begin(), a.begin() + m);
std::vector<T> a2(a.begin() + m, a.end());
std::vector<T> c(n);
std::vector<T> c1 = convolution_naive(a1, a1);
std::vector<T> c2 = convolution_naive(a1, a2);
std::vector<T> c3 = convolution_pow2(a2);
for (int i = 0; i < (int)c1.size(); i++) c[i] += c1[i];
for (int i = 0; i < (int)c2.size(); i++) c[i + m] += c2[i] * 2;
for (int i = 0; i < (int)c3.size(); i++) c[i + m * 2] += c3[i];
return c;
}
int m = 1 << lg;
a.resize(m);
number_theoretic_transform(a);
for (int i = 0; i < m; i++) a[i] *= a[i];
inverse_number_theoretic_transform(a);
a.resize(n);
return a;
}
template <class T>
std::vector<T> convolution(std::vector<T> a, std::vector<T> b) {
int n = (int)a.size() + (int)b.size() - 1;
if (n <= 1) return convolution_naive(a, b);
int lg = ceil_log2(n);
int m = 1 << lg;
if (n - (1 << (lg - 1)) <= 5) {
--lg;
if (a.size() < b.size()) std::swap(a, b);
int half_m = n - (1 << lg);
std::vector<T> a1(a.begin(), a.begin() + half_m);
std::vector<T> a2(a.begin() + half_m, a.end());
std::vector<T> c(n);
std::vector<T> c1 = convolution_naive(a1, b);
std::vector<T> c2 = convolution(a2, b);
for (int i = 0; i < (int)c1.size(); i++) c[i] += c1[i];
for (int i = 0; i < (int)c2.size(); i++) c[i + half_m] += c2[i];
return c;
}
a.resize(m);
b.resize(m);
number_theoretic_transform(a);
number_theoretic_transform(b);
for (int i = 0; i < m; i++) a[i] *= b[i];
inverse_number_theoretic_transform(a);
a.resize(n);
return a;
}
} // namespace internal
template <uint32_t p>
struct is_ntt_friendly : std::false_type {};
template <> struct is_ntt_friendly<998244353> : std::true_type {};
template <> struct is_ntt_friendly<469762049> : std::true_type {};
template <> struct is_ntt_friendly<1811939329> : std::true_type {};
template <> struct is_ntt_friendly<2013265921> : std::true_type {};
template <uint32_t p>
std::vector<StaticModInt<p>> convolution_for_any_mod(const std::vector<StaticModInt<p>>& a, const std::vector<StaticModInt<p>>& b);
template <uint32_t p>
std::vector<StaticModInt<p>> convolution_large(const std::vector<StaticModInt<p>>& a, const std::vector<StaticModInt<p>>& b);
template <uint32_t p>
std::vector<StaticModInt<p>> convolution(const std::vector<StaticModInt<p>>& a, const std::vector<StaticModInt<p>>& b) {
int n = (int)a.size(), m = (int)b.size();
if (n == 0 || m == 0) return {};
if (n <= 60 || m <= 60) return internal::convolution_naive(a, b);
if (n + m - 1 <= (int)((1 - p) & (p - 1))) {
if (n == m && a == b) return internal::convolution_pow2(a);
return internal::convolution(a, b);
}
if constexpr (is_ntt_friendly<p>::value) {
return convolution_large<p>(a, b);
}
return convolution_for_any_mod(a, b);
}
template <uint32_t p>
std::vector<StaticModInt<p>> convolution_for_any_mod(const std::vector<StaticModInt<p>>& a, const std::vector<StaticModInt<p>>& b) {
int n = (int)a.size(), m = (int)b.size();
assert(n + m - 1 <= (1 << 26));
static constexpr uint32_t MOD1 = 469762049;
static constexpr uint32_t MOD2 = 1811939329;
static constexpr uint32_t MOD3 = 2013265921;
using mint1 = StaticModInt<MOD1>;
using mint2 = StaticModInt<MOD2>;
using mint3 = StaticModInt<MOD3>;
std::vector<mint1> a1(n), b1(m);
std::vector<mint2> a2(n), b2(m);
std::vector<mint3> a3(n), b3(m);
for (int i = 0; i < n; i++) {
uint32_t v = a[i].get();
a1[i] = v; a2[i] = v; a3[i] = v;
}
for (int i = 0; i < m; i++) {
uint32_t v = b[i].get();
b1[i] = v; b2[i] = v; b3[i] = v;
}
auto c1 = internal::convolution(a1, b1);
auto c2 = internal::convolution(a2, b2);
auto c3 = internal::convolution(a3, b3);
static constexpr uint32_t INV1_2 = mint2(MOD1).inv().get();
static constexpr uint32_t INV1_3 = mint3(MOD1).inv().get();
static constexpr uint32_t INV2_3 = mint3(MOD2).inv().get();
std::vector<StaticModInt<p>> res(n + m - 1);
for (int i = 0; i < n + m - 1; i++) {
int64_t t1 = c1[i].get();
int64_t t2 = (c2[i].get() - t1 + MOD2) * INV1_2 % MOD2;
if (t2 < 0) t2 += MOD2;
int64_t t3 = ((c3[i].get() - t1 + MOD3) * INV1_3 % MOD3 - t2 + MOD3) * INV2_3 % MOD3;
if (t3 < 0) t3 += MOD3;
res[i] = StaticModInt<p>(t1 + (t2 + t3 * MOD2) % p * (MOD1 % p));
}
return res;
}
template <uint32_t p>
std::vector<StaticModInt<p>> convolution_large(const std::vector<StaticModInt<p>>& a, const std::vector<StaticModInt<p>>& b) {
int n = (int)a.size(), m = (int)b.size();
int limit = (int)((1 - p) & (p - 1));
int L = limit >> 1;
int k = (n + L - 1) / L;
int l = (m + L - 1) / L;
std::vector<std::vector<StaticModInt<p>>> A_hat(k, std::vector<StaticModInt<p>>(limit));
std::vector<std::vector<StaticModInt<p>>> B_hat(l, std::vector<StaticModInt<p>>(limit));
for (int i = 0; i < k; ++i) {
int start = i * L;
int end = std::min(n, start + L);
for (int j = start; j < end; ++j) {
A_hat[i][j - start] = a[j];
}
internal::number_theoretic_transform(A_hat[i]);
}
for (int i = 0; i < l; ++i) {
int start = i * L;
int end = std::min(m, start + L);
for (int j = start; j < end; ++j) {
B_hat[i][j - start] = b[j];
}
internal::number_theoretic_transform(B_hat[i]);
}
std::vector<StaticModInt<p>> res(n + m - 1);
for (int c = 0; c < k + l - 1; ++c) {
std::vector<StaticModInt<p>> C_chunk(limit);
for (int i = 0; i < k; ++i) {
int j = c - i;
if (0 <= j && j < l) {
for (int x = 0; x < limit; ++x) {
C_chunk[x] += A_hat[i][x] * B_hat[j][x];
}
}
}
internal::inverse_number_theoretic_transform(C_chunk);
int start = c * L;
int end = std::min(n + m - 1, start + limit);
for (int j = start; j < end; ++j) {
res[j] += C_chunk[j - start];
}
}
return res;
}
template <uint32_t p>
std::vector<int64_t> convolution(const std::vector<int64_t>& a, const std::vector<int64_t>& b) {
int n = (int)a.size(), m = (int)b.size();
std::vector<StaticModInt<p>> a2(n), b2(m);
for (int i = 0; i < n; i++) a2[i] = a[i];
for (int i = 0; i < m; i++) b2[i] = b[i];
auto c2 = convolution<p>(a2, b2);
std::vector<int64_t> c(c2.size());
for (int i = 0; i < (int)c2.size(); i++) c[i] = c2[i].get();
return c;
}
inline std::vector<int64_t> convolution_ll(const std::vector<int64_t>& a, const std::vector<int64_t>& b) {
int n = (int)a.size(), m = (int)b.size();
if (n == 0 || m == 0) return {};
static constexpr uint32_t MOD1 = 469762049;
static constexpr uint32_t MOD2 = 1811939329;
static constexpr uint32_t MOD3 = 2013265921;
using mint1 = StaticModInt<MOD1>;
using mint2 = StaticModInt<MOD2>;
using mint3 = StaticModInt<MOD3>;
std::vector<mint1> a1(n), b1(m);
std::vector<mint2> a2(n), b2(m);
std::vector<mint3> a3(n), b3(m);
for (int i = 0; i < n; i++) {
int64_t v = a[i] % MOD1; if (v < 0) v += MOD1; a1[i] = v;
v = a[i] % MOD2; if (v < 0) v += MOD2; a2[i] = v;
v = a[i] % MOD3; if (v < 0) v += MOD3; a3[i] = v;
}
for (int i = 0; i < m; i++) {
int64_t v = b[i] % MOD1; if (v < 0) v += MOD1; b1[i] = v;
v = b[i] % MOD2; if (v < 0) v += MOD2; b2[i] = v;
v = b[i] % MOD3; if (v < 0) v += MOD3; b3[i] = v;
}
auto c1 = internal::convolution(a1, b1);
auto c2 = internal::convolution(a2, b2);
auto c3 = internal::convolution(a3, b3);
static constexpr uint32_t INV1_2 = mint2(MOD1).inv().get();
static constexpr uint32_t INV1_3 = mint3(MOD1).inv().get();
static constexpr uint32_t INV2_3 = mint3(MOD2).inv().get();
std::vector<int64_t> res(n + m - 1);
for (int i = 0; i < n + m - 1; i++) {
int64_t t1 = c1[i].get();
int64_t t2 = c2[i].get() - t1;
t2 %= MOD2;
if (t2 < 0) t2 += MOD2;
t2 = (t2 * INV1_2) % MOD2;
int64_t t3 = c3[i].get() - t1;
t3 %= MOD3;
if (t3 < 0) t3 += MOD3;
t3 = (t3 * INV1_3) % MOD3;
t3 = (t3 - t2) % MOD3;
if (t3 < 0) t3 += MOD3;
t3 = (t3 * INV2_3) % MOD3;
__int128_t result = t1;
result += (__int128_t)t2 * MOD1;
result += (__int128_t)t3 * MOD1 * MOD2;
res[i] = (int64_t)result;
}
return res;
}
// https://judge.yosupo.jp/problem/convolution_mod_2_64
inline std::vector<uint64_t> convolution_64(const std::vector<uint64_t>& a, const std::vector<uint64_t>& b) {
int n = (int)a.size(), m = (int)b.size();
if (n == 0 || m == 0) return {};
std::vector<int64_t> a_lo(n), a_hi(n), b_lo(m), b_hi(m);
for (int i = 0; i < n; i++) {
a_lo[i] = a[i] & 0xFFFFFFFF;
a_hi[i] = a[i] >> 32;
}
for (int i = 0; i < m; i++) {
b_lo[i] = b[i] & 0xFFFFFFFF;
b_hi[i] = b[i] >> 32;
}
auto c_lolo = convolution_ll(a_lo, b_lo);
auto c_hilo = convolution_ll(a_hi, b_lo);
auto c_lohi = convolution_ll(a_lo, b_hi);
std::vector<uint64_t> res(n + m - 1);
for (int i = 0; i < n + m - 1; i++) {
uint64_t lolo = (uint64_t)c_lolo[i];
uint64_t hilo = (uint64_t)c_hilo[i];
uint64_t lohi = (uint64_t)c_lohi[i];
res[i] = lolo + (hilo << 32) + (lohi << 32);
}
return res;
}
template <class T>
void ntt_doubling_(std::vector<T>& a, std::vector<T> b) {
static constexpr internal::NthRoot<T> nth_root;
int n = (int)a.size();
const T z = nth_root.get(internal::msb(n) + 1);
T r = 1;
for (int i = 0; i < n; i++) {
b[i] *= r;
r *= z;
}
internal::number_theoretic_transform(b);
a.reserve(2 * n);
a.insert(a.end(), b.begin(), b.end());
}
template <class T>
void ntt_doubling_(std::vector<T>& a) {
auto b = a;
internal::inverse_number_theoretic_transform(b);
ntt_doubling_(a, std::move(b));
}
template <int ID>
std::vector<DynamicModInt<ID>> convolution(const std::vector<DynamicModInt<ID>>& a, const std::vector<DynamicModInt<ID>>& b) {
int n = (int)a.size(), m = (int)b.size();
if (n == 0 || m == 0) return {};
if (n <= 60 || m <= 60) return internal::convolution_naive(a, b);
uint32_t mod = DynamicModInt<ID>::get_mod();
if (n + m - 1 <= (int)((1 - mod) & (mod - 1))) {
if (n == m && a == b) return internal::convolution_pow2(a);
return internal::convolution(a, b);
}
static constexpr uint32_t MOD1 = 469762049;
static constexpr uint32_t MOD2 = 1811939329;
static constexpr uint32_t MOD3 = 2013265921;
using mint1 = StaticModInt<MOD1>;
using mint2 = StaticModInt<MOD2>;
using mint3 = StaticModInt<MOD3>;
std::vector<mint1> a1(n), b1(m);
std::vector<mint2> a2(n), b2(m);
std::vector<mint3> a3(n), b3(m);
for (int i = 0; i < n; i++) {
uint32_t v = a[i].get();
a1[i] = v; a2[i] = v; a3[i] = v;
}
for (int i = 0; i < m; i++) {
uint32_t v = b[i].get();
b1[i] = v; b2[i] = v; b3[i] = v;
}
auto c1 = internal::convolution(a1, b1);
auto c2 = internal::convolution(a2, b2);
auto c3 = internal::convolution(a3, b3);
static constexpr uint32_t INV1_2 = mint2(MOD1).inv().get();
static constexpr uint32_t INV1_3 = mint3(MOD1).inv().get();
static constexpr uint32_t INV2_3 = mint3(MOD2).inv().get();
std::vector<DynamicModInt<ID>> res(n + m - 1);
for (int i = 0; i < n + m - 1; i++) {
int64_t t1 = c1[i].get();
int64_t t2 = (c2[i].get() - t1 + MOD2) * INV1_2 % MOD2;
if (t2 < 0) t2 += MOD2;
int64_t t3 = ((c3[i].get() - t1 + MOD3) * INV1_3 % MOD3 - t2 + MOD3) * INV2_3 % MOD3;
if (t3 < 0) t3 += MOD3;
res[i] = DynamicModInt<ID>(t1 + (t2 + t3 * MOD2) % mod * MOD1);
}
return res;
}
} // namespace Convolution