hly1204's library
This documentation is automatically generated by online-judge-tools/verification-helper
test/convolution/2d_convolution.0.test.cpp
- View this file on GitHub
- Last update: 2025-01-19 15:28:01+08:00
- Problem: https://judge.yosupo.jp/problem/aplusb
Depends on
Code
#define PROBLEM "https://judge.yosupo.jp/problem/aplusb"
#include "ks.hpp"
#include "modint.hpp"
#include "random.hpp"
#include "rng.hpp"
#include <cassert>
#include <iostream>
#include <random>
#include <vector>
bool verify0() {
std::uniform_int_distribution<> dis(2, 50);
xoshiro256starstar gen(std::random_device{}());
using mint = ModInt<998244353>;
const int maxLenA = dis(gen);
const int maxLenB = dis(gen);
std::vector<std::vector<mint>> A(maxLenA);
std::vector<std::vector<mint>> B(maxLenB);
for (int i = 0; i < maxLenA; ++i) A[i] = random_vector<mint>(dis(gen));
for (int i = 0; i < maxLenB; ++i) B[i] = random_vector<mint>(dis(gen));
return convolution_2d_ks(A, B) == convolution_2d_naive(A, B);
}
bool verify1() {
std::uniform_int_distribution<> dis(2, 50);
xoshiro256starstar gen(std::random_device{}());
using mint = ModInt<998244353>;
const int maxLenA = dis(gen);
const int maxLenB = dis(gen);
std::vector<std::vector<mint>> A(maxLenA);
std::vector<std::vector<mint>> B(maxLenB);
for (int i = 0; i < maxLenA; ++i) A[i] = random_vector<mint>(dis(gen));
for (int i = 0; i < maxLenB; ++i) B[i] = random_vector<mint>(dis(gen));
return convolution_2d_ks_reciprocal(A, B) == convolution_2d_naive(A, B);
}
bool verify2() {
std::uniform_int_distribution<> dis(2, 50);
xoshiro256starstar gen(std::random_device{}());
using mint = ModInt<998244353>;
const int maxLenA = dis(gen);
const int maxLenB = dis(gen);
std::vector<std::vector<mint>> A(maxLenA);
std::vector<std::vector<mint>> B(maxLenB);
for (int i = 0; i < maxLenA; ++i) A[i] = random_vector<mint>(dis(gen));
for (int i = 0; i < maxLenB; ++i) B[i] = random_vector<mint>(dis(gen));
return convolution_2d_ks_negated(A, B) == convolution_2d_naive(A, B);
}
int main() {
std::ios::sync_with_stdio(false);
std::cin.tie(nullptr);
{
const bool ok = verify0();
if (!ok) return 1;
}
{
const bool ok = verify1();
if (!ok) return 1;
}
{
const bool ok = verify2();
if (!ok) return 1;
}
long long a, b;
std::cin >> a >> b;
std::cout << a + b;
return 0;
}#line 1 "test/convolution/2d_convolution.0.test.cpp"
#define PROBLEM "https://judge.yosupo.jp/problem/aplusb"
#line 2 "ks.hpp"
#line 2 "fft.hpp"
#include <algorithm>
#include <cassert>
#include <iterator>
#include <memory>
#include <vector>
template<typename Tp> class FftInfo {
static Tp least_quadratic_nonresidue() {
for (int i = 2;; ++i)
if (Tp(i).pow((Tp::mod() - 1) / 2) == -1) return Tp(i);
}
const int ordlog2_;
const Tp zeta_;
const Tp invzeta_;
const Tp imag_;
const Tp invimag_;
mutable std::vector<Tp> root_;
mutable std::vector<Tp> invroot_;
FftInfo()
: ordlog2_(__builtin_ctzll(Tp::mod() - 1)),
zeta_(least_quadratic_nonresidue().pow((Tp::mod() - 1) >> ordlog2_)),
invzeta_(zeta_.inv()), imag_(zeta_.pow(1LL << (ordlog2_ - 2))), invimag_(-imag_),
root_{Tp(1), imag_}, invroot_{Tp(1), invimag_} {}
public:
static const FftInfo &get() {
static FftInfo info;
return info;
}
Tp imag() const { return imag_; }
Tp inv_imag() const { return invimag_; }
Tp zeta() const { return zeta_; }
Tp inv_zeta() const { return invzeta_; }
const std::vector<Tp> &root(int n) const {
// [0, n)
assert((n & (n - 1)) == 0);
if (const int s = root_.size(); s < n) {
root_.resize(n);
for (int i = __builtin_ctz(s); (1 << i) < n; ++i) {
const int j = 1 << i;
root_[j] = zeta_.pow(1LL << (ordlog2_ - i - 2));
for (int k = j + 1; k < j * 2; ++k) root_[k] = root_[k - j] * root_[j];
}
}
return root_;
}
const std::vector<Tp> &inv_root(int n) const {
// [0, n)
assert((n & (n - 1)) == 0);
if (const int s = invroot_.size(); s < n) {
invroot_.resize(n);
for (int i = __builtin_ctz(s); (1 << i) < n; ++i) {
const int j = 1 << i;
invroot_[j] = invzeta_.pow(1LL << (ordlog2_ - i - 2));
for (int k = j + 1; k < j * 2; ++k) invroot_[k] = invroot_[k - j] * invroot_[j];
}
}
return invroot_;
}
};
inline int fft_len(int n) {
--n;
n |= n >> 1, n |= n >> 2, n |= n >> 4, n |= n >> 8;
return (n | n >> 16) + 1;
}
namespace detail {
template<typename Iterator> inline void
butterfly_n(Iterator a, int n,
const std::vector<typename std::iterator_traits<Iterator>::value_type> &root) {
assert(n > 0);
assert((n & (n - 1)) == 0);
const int bn = __builtin_ctz(n);
if (bn & 1) {
for (int i = 0; i < n / 2; ++i) {
const auto a0 = a[i], a1 = a[i + n / 2];
a[i] = a0 + a1, a[i + n / 2] = a0 - a1;
}
}
for (int i = n >> (bn & 1); i >= 4; i /= 4) {
const int i4 = i / 4;
for (int k = 0; k < i4; ++k) {
const auto a0 = a[k + i4 * 0], a1 = a[k + i4 * 1];
const auto a2 = a[k + i4 * 2], a3 = a[k + i4 * 3];
const auto a02p = a0 + a2, a02m = a0 - a2;
const auto a13p = a1 + a3, a13m = (a1 - a3) * root[1];
a[k + i4 * 0] = a02p + a13p, a[k + i4 * 1] = a02p - a13p;
a[k + i4 * 2] = a02m + a13m, a[k + i4 * 3] = a02m - a13m;
}
for (int j = i, m = 2; j < n; j += i, m += 2) {
const auto r = root[m], r2 = r * r, r3 = r2 * r;
for (int k = j; k < j + i4; ++k) {
const auto a0 = a[k + i4 * 0], a1 = a[k + i4 * 1] * r;
const auto a2 = a[k + i4 * 2] * r2, a3 = a[k + i4 * 3] * r3;
const auto a02p = a0 + a2, a02m = a0 - a2;
const auto a13p = a1 + a3, a13m = (a1 - a3) * root[1];
a[k + i4 * 0] = a02p + a13p, a[k + i4 * 1] = a02p - a13p;
a[k + i4 * 2] = a02m + a13m, a[k + i4 * 3] = a02m - a13m;
}
}
}
}
template<typename Iterator> inline void
inv_butterfly_n(Iterator a, int n,
const std::vector<typename std::iterator_traits<Iterator>::value_type> &root) {
assert(n > 0);
assert((n & (n - 1)) == 0);
const int bn = __builtin_ctz(n);
for (int i = 4; i <= (n >> (bn & 1)); i *= 4) {
const int i4 = i / 4;
for (int k = 0; k < i4; ++k) {
const auto a0 = a[k + i4 * 0], a1 = a[k + i4 * 1];
const auto a2 = a[k + i4 * 2], a3 = a[k + i4 * 3];
const auto a01p = a0 + a1, a01m = a0 - a1;
const auto a23p = a2 + a3, a23m = (a2 - a3) * root[1];
a[k + i4 * 0] = a01p + a23p, a[k + i4 * 1] = a01m + a23m;
a[k + i4 * 2] = a01p - a23p, a[k + i4 * 3] = a01m - a23m;
}
for (int j = i, m = 2; j < n; j += i, m += 2) {
const auto r = root[m], r2 = r * r, r3 = r2 * r;
for (int k = j; k < j + i4; ++k) {
const auto a0 = a[k + i4 * 0], a1 = a[k + i4 * 1];
const auto a2 = a[k + i4 * 2], a3 = a[k + i4 * 3];
const auto a01p = a0 + a1, a01m = a0 - a1;
const auto a23p = a2 + a3, a23m = (a2 - a3) * root[1];
a[k + i4 * 0] = a01p + a23p, a[k + i4 * 1] = (a01m + a23m) * r;
a[k + i4 * 2] = (a01p - a23p) * r2, a[k + i4 * 3] = (a01m - a23m) * r3;
}
}
}
if (bn & 1) {
for (int i = 0; i < n / 2; ++i) {
const auto a0 = a[i], a1 = a[i + n / 2];
a[i] = a0 + a1, a[i + n / 2] = a0 - a1;
}
}
}
} // namespace detail
// FFT_n: A(x) |-> bit-reversed order of [A(1), A(zeta_n), ..., A(zeta_n^(n-1))]
template<typename Iterator> inline void fft_n(Iterator a, int n) {
using Tp = typename std::iterator_traits<Iterator>::value_type;
detail::butterfly_n(a, n, FftInfo<Tp>::get().root(n / 2));
}
template<typename Tp> inline void fft(std::vector<Tp> &a) { fft_n(a.begin(), a.size()); }
// IFFT_n: bit-reversed order of [A(1), A(zeta_n), ..., A(zeta_n^(n-1))] |-> A(x)
template<typename Iterator> inline void inv_fft_n(Iterator a, int n) {
using Tp = typename std::iterator_traits<Iterator>::value_type;
detail::inv_butterfly_n(a, n, FftInfo<Tp>::get().inv_root(n / 2));
const Tp iv = Tp::mod() - (Tp::mod() - 1) / n;
for (int i = 0; i < n; ++i) a[i] *= iv;
}
template<typename Tp> inline void inv_fft(std::vector<Tp> &a) { inv_fft_n(a.begin(), a.size()); }
// IFFT_n^T: A(x) |-> 1/n FFT_n((x^n A(x^(-1))) mod (x^n - 1))
template<typename Iterator> inline void transposed_inv_fft_n(Iterator a, int n) {
using Tp = typename std::iterator_traits<Iterator>::value_type;
const Tp iv = Tp::mod() - (Tp::mod() - 1) / n;
for (int i = 0; i < n; ++i) a[i] *= iv;
detail::butterfly_n(a, n, FftInfo<Tp>::get().inv_root(n / 2));
}
template<typename Tp> inline void transposed_inv_fft(std::vector<Tp> &a) {
transposed_inv_fft_n(a.begin(), a.size());
}
// FFT_n^T : FFT_n((x^n A(x^(-1))) mod (x^n - 1)) |-> n A(x)
template<typename Iterator> inline void transposed_fft_n(Iterator a, int n) {
using Tp = typename std::iterator_traits<Iterator>::value_type;
detail::inv_butterfly_n(a, n, FftInfo<Tp>::get().root(n / 2));
}
template<typename Tp> inline void transposed_fft(std::vector<Tp> &a) {
transposed_fft_n(a.begin(), a.size());
}
template<typename Tp> inline std::vector<Tp> convolution_fft(std::vector<Tp> a, std::vector<Tp> b) {
if (a.empty() || b.empty()) return {};
const int n = a.size();
const int m = b.size();
const int len = fft_len(n + m - 1);
a.resize(len);
b.resize(len);
fft(a);
fft(b);
for (int i = 0; i < len; ++i) a[i] *= b[i];
inv_fft(a);
a.resize(n + m - 1);
return a;
}
template<typename Tp> inline std::vector<Tp> square_fft(std::vector<Tp> a) {
if (a.empty()) return {};
const int n = a.size();
const int len = fft_len(n * 2 - 1);
a.resize(len);
fft(a);
for (int i = 0; i < len; ++i) a[i] *= a[i];
inv_fft(a);
a.resize(n * 2 - 1);
return a;
}
template<typename Tp>
inline std::vector<Tp> convolution_naive(const std::vector<Tp> &a, const std::vector<Tp> &b) {
if (a.empty() || b.empty()) return {};
const int n = a.size();
const int m = b.size();
std::vector<Tp> res(n + m - 1);
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 Tp>
inline std::vector<Tp> convolution(const std::vector<Tp> &a, const std::vector<Tp> &b) {
if (std::min(a.size(), b.size()) < 60) return convolution_naive(a, b);
if (std::addressof(a) == std::addressof(b)) return square_fft(a);
return convolution_fft(a, b);
}
#line 7 "ks.hpp"
// returns max[0 <= j < a.size()]{a[j].size()}
template<typename Tp> inline int max_len_x_ks(const std::vector<std::vector<Tp>> &a) {
int len = -1;
for (int i = 0; i < (int)a.size(); ++i) len = std::max<int>(len, a[i].size());
return len;
}
// returns a(x, x^N) where a(x,y) in R[x][y]
template<typename Tp>
inline std::vector<Tp> pack_2d_ks(const std::vector<std::vector<Tp>> &a, int N) {
assert(N > 0);
// y |-> x^N
std::vector<Tp> b(N * a.size());
for (int i = 0; i < (int)a.size(); ++i)
for (int j = 0; j < (int)a[i].size(); ++j) b[i * N + j] += a[i][j]; // if N < a[i].size()
return b;
}
template<typename Tp>
inline std::vector<std::vector<Tp>> unpack_2d_ks(const std::vector<Tp> &a, int N) {
assert(N > 0);
// x^N |-> y
std::vector<std::vector<Tp>> b((a.size() + (N - 1)) / N, std::vector<Tp>(N));
for (int i = 0; i < (int)((a.size() + N - 1) / N); ++i)
for (int j = 0; j < N && i * N + j < (int)a.size(); ++j) b[i][j] = a[i * N + j];
return b;
}
template<typename Tp> inline std::vector<std::vector<Tp>>
convolution_2d_naive(const std::vector<std::vector<Tp>> &a, const std::vector<std::vector<Tp>> &b) {
if (a.empty() || b.empty()) return {};
const int lenA = max_len_x_ks(a);
const int lenB = max_len_x_ks(b);
if (lenA == 0 || lenB == 0) return std::vector<std::vector<Tp>>(a.size() + b.size() - 1);
const int N = lenA + lenB - 1;
std::vector<std::vector<Tp>> ab(a.size() + b.size() - 1, std::vector<Tp>(N));
for (int i = 0; i < (int)a.size(); ++i)
for (int j = 0; j < (int)b.size(); ++j) {
const auto aibj = convolution(a[i], b[j]);
for (int k = 0; k < (int)aibj.size(); ++k) ab[i + j][k] += aibj[k];
}
return ab;
}
// standard Kronecker substitution
template<typename Tp> inline std::vector<std::vector<Tp>>
convolution_2d_ks(const std::vector<std::vector<Tp>> &a, const std::vector<std::vector<Tp>> &b) {
if (a.empty() || b.empty()) return {};
const int lenA = max_len_x_ks(a);
const int lenB = max_len_x_ks(b);
if (lenA == 0 || lenB == 0) return std::vector<std::vector<Tp>>(a.size() + b.size() - 1);
const int N = lenA + lenB - 1;
auto ab = convolution(pack_2d_ks(a, N), pack_2d_ks(b, N));
ab.resize((a.size() + b.size() - 1) * N);
return unpack_2d_ks(ab, N);
}
// see:
// [1]: David Harvey. Faster polynomial multiplication via multipoint Kronecker substitution.
// https://doi.org/10.1016/j.jsc.2009.05.004
template<typename Tp> inline std::vector<std::vector<Tp>>
convolution_2d_ks_reciprocal(const std::vector<std::vector<Tp>> &a,
const std::vector<std::vector<Tp>> &b) {
if (a.empty() || b.empty()) return {};
const int lenA = max_len_x_ks(a);
const int lenB = max_len_x_ks(b);
if (lenA == 0 || lenB == 0) return std::vector<std::vector<Tp>>(a.size() + b.size() - 1);
const int N = std::max(lenA, lenB);
// original version: compute a(x, x^(-N)) b(x, x^(-N))
// modified version (this version): compute x^(2N-2) a(x^(-1), x^N) b(x^(-1), x^N)
// ab0 = a(x, x^N) b(x, x^N)
auto ab0 = convolution(pack_2d_ks(a, N), pack_2d_ks(b, N));
// returns x^(N-1) a(x^(-1), y)
auto make_reciprocal = [](const std::vector<std::vector<Tp>> &a, int N) {
std::vector<std::vector<Tp>> b(a.size());
for (int i = 0; i < (int)a.size(); ++i) {
b[i] = a[i];
b[i].resize(N);
std::reverse(b[i].begin(), b[i].end());
}
return b;
};
// ab1 = x^(2N-2) a(x^(-1), x^N) b(x^(-1), x^N)
auto ab1 =
convolution(pack_2d_ks(make_reciprocal(a, N), N), pack_2d_ks(make_reciprocal(b, N), N));
std::vector<std::vector<Tp>> ab(a.size() + b.size() - 1, std::vector<Tp>(N * 2 - 1));
/*
example:
let pack: a(x,y) |-> a(x,x^3)
rec: a(x,y) |-> x^2 a(x^(-1),y)
a = 1 + 2*x + 4*x^2 + (3 + 4*x)*y
b = 2 + 4*x + 5*x^2 + (5 + 8*x)*y
pack(a) = 1 + 2*x + 4*x^2 + 3*x^3 + 4*x^4
pack(b) = 2 + 4*x + 5*x^2 + 5*x^3 + 8*x^4
pack(rec(a)) = 4 + 2*x + x^2 + (0 + 4*x + 3*x^2)*y
pack(rec(b)) = 5 + 4*x + 2*x^2 + (0 + 8*x + 5*x^2)*y
ab = 2 + 8*x + 21*x^2 + 26*x^3 + 20*x^4 +
(11 + 38*x + 67*x^2 + 52*x^3)*y +
+ (15 + 44*x + 32*x^2)*y^2
pack(a)pack(b) = 2 + 8*x + 21*x^2 + 37*x^3 + 58*x^4 + 67*x^5 + 67*x^6 + 44*x^7 + 32*x^8
pack(rec(a))pack(rec(b))
= 20 + 26*x + 21*x^2 + 8*x^3 + 2*x^4 +
+ ... (overlap)
*/
// restore ab[0]
for (int i = 0; i < N; ++i) ab[0][i] = ab0[i];
// ab1[0] = [x^0](x^(2N - 2) a(x^(-1), x^N) b(x^(-1), x^N))
for (int i = 0; i < N; ++i) ab[0][(N - 1) * 2 - i] = ab1[i];
// restore ab[1..] by subtracting the overlaped coefficients
for (int i = 1; i < (int)(a.size() + b.size() - 1); ++i) {
// TODO: remove redundant assignment/subtraction
for (int j = 0; j < N * 2 - 1; ++j) {
ab0[(i - 1) * N + j] -= ab[i - 1][j];
ab1[(i - 1) * N + j] -= ab[i - 1][(N - 1) * 2 - j];
}
for (int j = 0; j < N; ++j) ab[i][j] = ab0[i * N + j];
for (int j = 0; j < N; ++j) ab[i][(N - 1) * 2 - j] = ab1[i * N + j];
}
for (int i = 0; i < (int)(a.size() + b.size() - 1); ++i) ab[i].resize(lenA + lenB - 1);
return ab;
}
// see:
// [1]: David Harvey. Faster polynomial multiplication via multipoint Kronecker substitution.
// https://doi.org/10.1016/j.jsc.2009.05.004
template<typename Tp> inline std::vector<std::vector<Tp>>
convolution_2d_ks_negated(const std::vector<std::vector<Tp>> &a,
const std::vector<std::vector<Tp>> &b) {
if (a.empty() || b.empty()) return {};
const int lenA = max_len_x_ks(a);
const int lenB = max_len_x_ks(b);
if (lenA == 0 || lenB == 0) return std::vector<std::vector<Tp>>(a.size() + b.size() - 1);
const int N = std::max(lenA, lenB);
// ab0 = a(x, x^N) b(x, x^N)
const auto ab0 = convolution(pack_2d_ks(a, N), pack_2d_ks(b, N));
// returns a(x, -y)
auto make_negated = [](const std::vector<std::vector<Tp>> &a) {
auto b = a;
for (int i = 1; i < (int)b.size(); i += 2)
for (int j = 0; j < (int)b[i].size(); ++j) b[i][j] = -b[i][j];
return b;
};
// ab1 = a(x, -x^N) b(x, -x^N)
const auto ab1 = convolution(pack_2d_ks(make_negated(a), N), pack_2d_ks(make_negated(b), N));
/*
example:
let pack: a(x,y) |-> a(x,x^3)
a = 1 + 2*x + 4*x^2 + (3 + 4*x)*y
b = 2 + 4*x + 5*x^2 + (5 + 8*x)*y
pack(a) = 1 + 2*x + 4*x^2 + 3*x^3 + 4*x^4
pack(b) = 2 + 4*x + 5*x^2 + 5*x^3 + 8*x^4
pack(a(x,-y)) = 1 + 2*x + 4*x^2 + -(3*x^3 + 4*x^4)
pack(b(x,-y)) = 2 + 4*x + 5*x^2 + -(5*x^3 + 8*x^4)
ab = 2 + 8*x + 21*x^2 + 26*x^3 + 20*x^4 +
(11 + 38*x + 67*x^2 + 52*x^3)*y +
+ (15 + 44*x + 32*x^2)*y^2
pack(a)pack(b) = 2 + 8*x + 21*x^2 + 37*x^3 + 58*x^4 + 67*x^5 + 67*x^6 + 44*x^7 + 32*x^8
(26+11) (20+38) ...
pack(a(x,-y))pack(b(x,-y))
= 2 + 8*x + 21*x^2 + 15*x^3 + -18*x^4 + ...
(26-11) (20-38) ...
*/
std::vector<std::vector<Tp>> ab(a.size() + b.size() - 1, std::vector<Tp>(lenA + lenB - 1));
for (int i = 0; i < (int)(a.size() + b.size() - 1); ++i) {
if (i & 1) {
for (int j = 0; j < lenA + lenB - 1; ++j)
ab[i][j] = (ab0[i * N + j] - ab1[i * N + j]).div_by_2();
} else {
for (int j = 0; j < lenA + lenB - 1; ++j)
ab[i][j] = (ab0[i * N + j] + ab1[i * N + j]).div_by_2();
}
}
return ab;
}
#line 2 "modint.hpp"
#include <iostream>
#include <type_traits>
// clang-format off
template<unsigned Mod> class ModInt {
static_assert((Mod >> 31) == 0, "`Mod` must less than 2^(31)");
template<typename Int>
static std::enable_if_t<std::is_integral_v<Int>, unsigned> safe_mod(Int v) { using D = std::common_type_t<Int, unsigned>; return (v %= (int)Mod) < 0 ? (D)(v + (int)Mod) : (D)v; }
struct PrivateConstructor {} static inline private_constructor{};
ModInt(PrivateConstructor, unsigned v) : v_(v) {}
unsigned v_;
public:
static unsigned mod() { return Mod; }
static ModInt from_raw(unsigned v) { return ModInt(private_constructor, v); }
static ModInt zero() { return from_raw(0); }
static ModInt one() { return from_raw(1); }
bool is_zero() const { return v_ == 0; }
bool is_one() const { return v_ == 1; }
ModInt() : v_() {}
template<typename Int, typename std::enable_if_t<std::is_signed_v<Int>, int> = 0> ModInt(Int v) : v_(safe_mod(v)) {}
template<typename Int, typename std::enable_if_t<std::is_unsigned_v<Int>, int> = 0> ModInt(Int v) : v_(v % Mod) {}
unsigned val() const { return v_; }
ModInt operator-() const { return from_raw(v_ == 0 ? v_ : Mod - v_); }
ModInt pow(long long e) const { if (e < 0) return inv().pow(-e); for (ModInt x(*this), res(from_raw(1));; x *= x) { if (e & 1) res *= x; if ((e >>= 1) == 0) return res; }}
ModInt inv() const { int x1 = 1, x3 = 0, a = val(), b = Mod; while (b) { const int q = a / b, x1_old = x1, a_old = a; x1 = x3, x3 = x1_old - x3 * q, a = b, b = a_old - b * q; } return from_raw(x1 < 0 ? x1 + (int)Mod : x1); }
template<bool Odd = (Mod & 1)> std::enable_if_t<Odd, ModInt> div_by_2() const { if (v_ & 1) return from_raw((v_ + Mod) >> 1); return from_raw(v_ >> 1); }
ModInt &operator+=(const ModInt &a) { if ((v_ += a.v_) >= Mod) v_ -= Mod; return *this; }
ModInt &operator-=(const ModInt &a) { if ((v_ += Mod - a.v_) >= Mod) v_ -= Mod; return *this; }
ModInt &operator*=(const ModInt &a) { v_ = (unsigned long long)v_ * a.v_ % Mod; return *this; }
ModInt &operator/=(const ModInt &a) { return *this *= a.inv(); }
ModInt &operator++() { return *this += one(); }
ModInt operator++(int) { ModInt o(*this); *this += one(); return o; }
ModInt &operator--() { return *this -= one(); }
ModInt operator--(int) { ModInt o(*this); *this -= one(); return o; }
friend ModInt operator+(const ModInt &a, const ModInt &b) { return ModInt(a) += b; }
friend ModInt operator-(const ModInt &a, const ModInt &b) { return ModInt(a) -= b; }
friend ModInt operator*(const ModInt &a, const ModInt &b) { return ModInt(a) *= b; }
friend ModInt operator/(const ModInt &a, const ModInt &b) { return ModInt(a) /= b; }
friend bool operator==(const ModInt &a, const ModInt &b) { return a.v_ == b.v_; }
friend bool operator!=(const ModInt &a, const ModInt &b) { return a.v_ != b.v_; }
friend std::istream &operator>>(std::istream &a, ModInt &b) { int v; a >> v; b.v_ = safe_mod(v); return a; }
friend std::ostream &operator<<(std::ostream &a, const ModInt &b) { return a << b.val(); }
};
// clang-format on
#line 2 "random.hpp"
#line 2 "rng.hpp"
#include <cstdint>
#include <limits>
// see: https://prng.di.unimi.it/xoshiro256starstar.c
// original license CC0 1.0
class xoshiro256starstar {
using u64 = std::uint64_t;
static inline u64 rotl(const u64 x, int k) { return (x << k) | (x >> (64 - k)); }
u64 s_[4];
u64 next() {
const u64 res = rotl(s_[1] * 5, 7) * 9;
const u64 t = s_[1] << 17;
s_[2] ^= s_[0];
s_[3] ^= s_[1];
s_[1] ^= s_[2];
s_[0] ^= s_[3];
s_[2] ^= t;
s_[3] = rotl(s_[3], 45);
return res;
}
public:
// see: https://prng.di.unimi.it/splitmix64.c
// original license CC0 1.0
explicit xoshiro256starstar(u64 seed) {
for (int i = 0; i < 4; ++i) {
u64 z = (seed += 0x9e3779b97f4a7c15);
z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9;
z = (z ^ (z >> 27)) * 0x94d049bb133111eb;
s_[i] = z ^ (z >> 31);
}
}
// see: https://en.cppreference.com/w/cpp/named_req/UniformRandomBitGenerator
using result_type = u64;
static constexpr u64 min() { return std::numeric_limits<u64>::min(); }
static constexpr u64 max() { return std::numeric_limits<u64>::max(); }
u64 operator()() { return next(); }
};
#line 4 "random.hpp"
#include <random>
#line 6 "random.hpp"
template<typename Tp> inline std::vector<Tp> random_vector(int n) {
std::vector<Tp> res(n);
xoshiro256starstar rng(std::random_device{}());
std::uniform_int_distribution<decltype(Tp::mod())> dis(0, Tp::mod() - 1);
for (int i = 0; i < n; ++i) res[i] = dis(rng);
return res;
}
template<typename Tp> inline std::vector<Tp> random_vector_without_zero(int n) {
std::vector<Tp> res(n);
xoshiro256starstar rng(std::random_device{}());
std::uniform_int_distribution<decltype(Tp::mod())> dis(1, Tp::mod() - 1);
for (int i = 0; i < n; ++i) res[i] = dis(rng);
return res;
}
#line 11 "test/convolution/2d_convolution.0.test.cpp"
bool verify0() {
std::uniform_int_distribution<> dis(2, 50);
xoshiro256starstar gen(std::random_device{}());
using mint = ModInt<998244353>;
const int maxLenA = dis(gen);
const int maxLenB = dis(gen);
std::vector<std::vector<mint>> A(maxLenA);
std::vector<std::vector<mint>> B(maxLenB);
for (int i = 0; i < maxLenA; ++i) A[i] = random_vector<mint>(dis(gen));
for (int i = 0; i < maxLenB; ++i) B[i] = random_vector<mint>(dis(gen));
return convolution_2d_ks(A, B) == convolution_2d_naive(A, B);
}
bool verify1() {
std::uniform_int_distribution<> dis(2, 50);
xoshiro256starstar gen(std::random_device{}());
using mint = ModInt<998244353>;
const int maxLenA = dis(gen);
const int maxLenB = dis(gen);
std::vector<std::vector<mint>> A(maxLenA);
std::vector<std::vector<mint>> B(maxLenB);
for (int i = 0; i < maxLenA; ++i) A[i] = random_vector<mint>(dis(gen));
for (int i = 0; i < maxLenB; ++i) B[i] = random_vector<mint>(dis(gen));
return convolution_2d_ks_reciprocal(A, B) == convolution_2d_naive(A, B);
}
bool verify2() {
std::uniform_int_distribution<> dis(2, 50);
xoshiro256starstar gen(std::random_device{}());
using mint = ModInt<998244353>;
const int maxLenA = dis(gen);
const int maxLenB = dis(gen);
std::vector<std::vector<mint>> A(maxLenA);
std::vector<std::vector<mint>> B(maxLenB);
for (int i = 0; i < maxLenA; ++i) A[i] = random_vector<mint>(dis(gen));
for (int i = 0; i < maxLenB; ++i) B[i] = random_vector<mint>(dis(gen));
return convolution_2d_ks_negated(A, B) == convolution_2d_naive(A, B);
}
int main() {
std::ios::sync_with_stdio(false);
std::cin.tie(nullptr);
{
const bool ok = verify0();
if (!ok) return 1;
}
{
const bool ok = verify1();
if (!ok) return 1;
}
{
const bool ok = verify2();
if (!ok) return 1;
}
long long a, b;
std::cin >> a >> b;
std::cout << a + b;
return 0;
}