#line 1 "000.cpp" #include using namespace std::literals::string_literals; using i64 = std::int_fast64_t; using std::cout; using std::cerr; using std::endl; using std::cin; #if defined(DONLINE_JUDGE) #define NDEBUG #elif defined(ONLINE_JUDGE) #define NDEBUG #endif template std::vector make_v(size_t a){return std::vector(a);} template auto make_v(size_t a, Ts... ts){ return std::vector(ts...))>(a, make_v(ts...)); } #line 1 "/home/ecasdqina/cpcpp/libs/library_cpp/math/formal_power_series.hpp" #line 8 "/home/ecasdqina/cpcpp/libs/library_cpp/math/formal_power_series.hpp" namespace cplib { template class formal_power_series: public T { public: using T::T; using value_type = typename T::value_type; using reference = typename T::reference; using const_reference = typename T::const_reference; using size_type = typename T::size_type; using i64 = std::int_fast64_t; public: formal_power_series(const T& p): T(p) {} public: formal_power_series inverse() const { assert((*this)[0] != value_type{}); formal_power_series ret(1, (*this)[0].inverse()); for(size_t i = 1; i < this->size(); i <<= 1) { auto tmp = T::convolution(ret, this->prefix(i << 1)); for(size_t j = 0; j < i; j++) { tmp[j] = value_type{}; if(j + i < tmp.size()) tmp[j + i] *= value_type(-1); } tmp = tmp * ret; for(size_t j = 0; j < i; j++) tmp[j] = ret[j]; ret = std::move(tmp).prefix(i << 1); } return ret.prefix(this->size()); } formal_power_series log() const { assert((*this)[0] == value_type(1)); return (formal_power_series(this->differential()) * this->inverse()).integral().prefix(this->size()); } formal_power_series exp() const { assert((*this)[0] == value_type{}); formal_power_series f(1, value_type(1)), g(1, value_type(1)); for(size_t i = 1; i < this->size(); i <<= 1) { g = (g * value_type(2) - f * g * g).prefix(i); formal_power_series q = this->differential().prefix(i - 1); formal_power_series w = (q + g * (f.differential() - f * q)).prefix((i << 1) - 1); f = (f + f * (*this - w.integral()).prefix(i << 1)).prefix(i << 1); } return f.prefix(this->size()); } formal_power_series power(i64 k) const { if(k < 0) return this->inverse().power(-k); for(size_type i = 0; i < this->size(); i++) { if((*this)[i] != value_type{}) { value_type inv = (*this)[i].inverse(); formal_power_series f(*this * inv); formal_power_series g(f >> i); g = formal_power_series(g.log() * value_type(k)).exp() * (*this)[i].power(k); if(i * k > this->size()) return formal_power_series(this->size()); return (g << (i * k)).prefix(this->size()); } } return *this; } // e^{e^x - 1}.prefix(n) constexpr formal_power_series bell_number(size_t n) { formal_power_series f(n); f[0] = 1; for(size_t i = 0; i < n; i++) f[i + 1] = f[i] / (i + 1); return (f - 1).exp(); } public: formal_power_series convolute(const formal_power_series& r) const { return ((*this) * r).prefix(this->size()); } }; } // @docs docs/formal_power_series.md #line 1 "/home/ecasdqina/cpcpp/libs/library_cpp/math/number_theoritic_transform.hpp" #line 1 "/home/ecasdqina/cpcpp/libs/library_cpp/math/modint.hpp" #line 5 "/home/ecasdqina/cpcpp/libs/library_cpp/math/modint.hpp" namespace cplib { template class modint { using u32 = std::uint_fast32_t; using u64 = std::uint_fast64_t; using i32 = std::int_fast32_t; using i64 = std::int_fast64_t; inline u64 apply(i64 x) { return (x < 0 ? x + Modulus : x); }; public: u64 a; static constexpr u64 mod = Modulus; constexpr modint(const i64& x = 0) noexcept: a(apply(x % (i64)Modulus)) {} constexpr modint operator+(const modint& rhs) const noexcept { return modint(*this) += rhs; } constexpr modint operator-(const modint& rhs) const noexcept { return modint(*this) -= rhs; } constexpr modint operator*(const modint& rhs) const noexcept { return modint(*this) *= rhs; } constexpr modint operator/(const modint& rhs) const noexcept { return modint(*this) /= rhs; } constexpr modint operator^(const u64& k) const noexcept { return modint(*this) ^= k; } constexpr modint operator^(const modint& k) const noexcept { return modint(*this) ^= k.value(); } constexpr modint operator-() const noexcept { return modint(Modulus - a); } constexpr modint operator++() noexcept { return (*this) = modint(*this) + 1; } constexpr modint operator--() noexcept { return (*this) = modint(*this) - 1; } const bool operator==(const modint& rhs) const noexcept { return a == rhs.a; }; const bool operator!=(const modint& rhs) const noexcept { return a != rhs.a; }; const bool operator<=(const modint& rhs) const noexcept { return a <= rhs.a; }; const bool operator>=(const modint& rhs) const noexcept { return a >= rhs.a; }; const bool operator<(const modint& rhs) const noexcept { return a < rhs.a; }; const bool operator>(const modint& rhs) const noexcept { return a > rhs.a; }; constexpr modint& operator+=(const modint& rhs) noexcept { a += rhs.a; if (a >= Modulus) a -= Modulus; return *this; } constexpr modint& operator-=(const modint& rhs) noexcept { if (a < rhs.a) a += Modulus; a -= rhs.a; return *this; } constexpr modint& operator*=(const modint& rhs) noexcept { a = a * rhs.a % Modulus; return *this; } constexpr modint& operator/=(modint rhs) noexcept { u64 exp = Modulus - 2; while (exp) { if (exp % 2) (*this) *= rhs; rhs *= rhs; exp /= 2; } return *this; } constexpr modint& operator^=(u64 k) noexcept { auto b = modint(1); while(k) { if(k & 1) b = b * (*this); (*this) *= (*this); k >>= 1; } return (*this) = b; } constexpr modint& operator=(const modint& rhs) noexcept { a = rhs.a; return (*this); } const modint inverse() const { return modint(1) / *this; } const modint power(i64 k) const { if(k < 0) return modint(*this).inverse() ^ (-k); return modint(*this) ^ k; } explicit operator bool() const { return a; } explicit operator u64() const { return a; } constexpr u64& value() noexcept { return a; } constexpr const u64& value() const noexcept { return a; } friend std::ostream& operator<<(std::ostream& os, const modint& p) { return os << p.a; } friend std::istream& operator>>(std::istream& is, modint& p) { u64 t; is >> t; p = modint(t); return is; } }; } #line 1 "/home/ecasdqina/cpcpp/libs/library_cpp/math/polynomial.hpp" #line 8 "/home/ecasdqina/cpcpp/libs/library_cpp/math/polynomial.hpp" namespace cplib { template class polynomial: public std::vector { public: using std::vector::vector; using value_type = typename std::vector::value_type; using reference = typename std::vector::reference; using const_reference = typename std::vector::const_reference; using size_type = typename std::vector::size_type; public: polynomial(const std::vector& r): std::vector(r) {} polynomial(size_t size, std::function f): std::vector(size) { for(size_t i = 0; i < size; i++) (*this)[i] = f(i); } polynomial operator+(const polynomial& r) const { return polynomial(*this) += r; } polynomial operator+(const_reference r) const { return polynomial(*this) += r; } friend polynomial operator+(const_reference l, polynomial r) { return r += l; } polynomial operator-(const polynomial& r) const { return polynomial(*this) -= r; } polynomial operator-(const_reference r) const { return polynomial(*this) -= r; } friend polynomial operator-(const_reference l, polynomial r) { return r -= l; } polynomial operator*(const_reference r) const { return polynomial(*this) *= r; } friend polynomial operator*(const_reference l, polynomial r) { return r *= l; } polynomial operator*(const polynomial& r) const { return polynomial(*this) *= r; } polynomial operator/(const_reference r) const { return polynomial(*this) /= r; } polynomial operator<<(size_type r) const { return polynomial(*this) <<= r; } polynomial operator>>(size_type r) const { return polynomial(*this) >>= r; } polynomial operator-() const { polynomial ret(this->size()); for(size_t i = 0; i < this->size(); i++) ret[i] = -(*this)[i]; return ret; } polynomial& operator+=(const polynomial& r) { if(r.size() > this->size()) this->resize(r.size()); for(size_t i = 0; i < r.size(); i++) (*this)[i] = (*this)[i] + r[i]; return *this; } polynomial& operator+=(const_reference r) { if(this->empty()) this->push_back(T{}); (*this)[0] += r; return *this; } polynomial& operator-=(const polynomial& r) { if(r.size() > this->size()) this->resize(r.size()); for(size_t i = 0; i < r.size(); i++) (*this)[i] = (*this)[i] - r[i]; return *this; } polynomial& operator-=(const_reference r) { if(this->empty()) this->push_back(T{}); (*this)[0] -= r; return *this; } polynomial& operator*=(const_reference r) { for(size_t i = 0; i < this->size(); i++) (*this)[i] = (*this)[i] * r; return *this; } polynomial& operator/=(const_reference r) { for(size_t i = 0; i < this->size(); i++) (*this)[i] = (*this)[i] / r; return *this; } polynomial& operator<<=(size_type r) { this->insert(begin(*this), r, T{}); return *this; } polynomial& operator>>=(size_type r) { if(r >= this->size()) this->clear(); else this->erase(begin(*this), begin(*this) + r); return *this; } polynomial& operator*=(const polynomial& r) { polynomial ret(this->degree() + r.degree() + 1); for(size_t i = 0; i < this->size(); i++) { if((*this)[i] == T{}) continue; for(size_t j = 0; j < r.size(); j++) ret[i + j] += (*this)[i] * r[j]; } return *this = std::move(ret); }; polynomial differential(size_type k) const { polynomial ret(*this); for(size_t i = 0; i < k; i++) ret = ret.differential(); return ret; } polynomial differential() const { if(degree() < 1) return polynomial(); polynomial ret(this->size() - 1); for(size_t i = 1; i < this->size(); i++) ret[i - 1] = (*this)[i] * T(i); return ret; } polynomial integral(size_type k) const { polynomial ret(*this); for(size_t i = 0; i < k; i++) ret = ret.integral(); return ret; } polynomial integral() const { polynomial ret(this->size() + 1); for(size_t i = 0; i < this->size(); i++) ret[i + 1] = (*this)[i] / T(i + 1); return ret; } polynomial prefix(size_type size) const { return polynomial(begin(*this), begin(*this) + std::min(this->size(), size)); } void shrink() { while(!this->empty() and this->back() == T{}) this->pop_back(); } public: polynomial sparse_convex(const polynomial& r) const { if(this->nonzeros() < r.nonzeros()) return (*this) * r; return r * (*this); } size_type degree() const { return std::max(1, this->size()) - 1; } size_type nonzeros() const { size_type ret = 0; for(size_t i = 0; i < this->size(); i++) { if((*this)[i] != T{}) ret++; } return ret; } T operator()(T x) const { return eval(x); } T eval(T x) const { T ret = (*this)[0], tmp = x; for(size_t i = 1; i < this->size(); i++) { ret = ret + (*this)[i] * tmp; tmp = tmp * x; } return ret; } static polynomial one() { return term(0); } static polynomial zero() { return polynomial({T{0}}); } static polynomial term(int k) { return polynomial({T{1}}) << k; } }; template T convex_all(std::vector polies, int size = -1) { if(polies.empty()) return T::zero(); std::deque qu; for(int i = 0; i < polies.size(); i++) qu.push_back(i); while(qu.size() > 1) { int a = qu.front(); qu.pop_front(); int b = qu.front(); qu.pop_front(); polies.push_back(polies[a] * polies[b]); if(size != -1) polies.back().resize(size); qu.push_back((int)polies.size() - 1); } return polies.back(); } } #line 7 "/home/ecasdqina/cpcpp/libs/library_cpp/math/number_theoritic_transform.hpp" namespace cplib { template> class number_theoritic_transform: public polynomial { public: using polynomial::polynomial; using value_type = typename polynomial::value_type; using reference = typename polynomial::reference; using const_reference = typename polynomial::const_reference; using size_type = typename polynomial::size_type; public: number_theoritic_transform(const polynomial& p): polynomial(p) {} protected: static void ntt(number_theoritic_transform& a) { size_t N = a.size(); static std::vector dw; if(dw.size() < N) { int n = dw.size(); dw.resize(N); for(size_t i = n; i < N; i++) dw[i] = -(T(primitive_root) ^ ((T::mod - 1) >> (i + 2))); } for(size_t m = N; m >>= 1;) { T w = 1; for(size_t s = 0, k = 0; s < N; s += 2 * m) { for(size_t i = s, j = s + m; i < s + m; i++, j++) { T x = a[i], y = a[j] * w; a[i] = x + y; a[j] = x - y; } w *= dw[__builtin_ctz(++k)]; } } } static void intt(number_theoritic_transform& a) { size_t N = a.size(); static std::vector idw; if(idw.size() < N) { size_t n = idw.size(); idw.resize(N); for(size_t i = n; i < N; i++) idw[i] = (-(T(primitive_root) ^ ((T::mod - 1) >> (i + 2)))).inverse(); } for(size_t m = 1; m < N; m *= 2) { T w = 1; for(size_t s = 0, k = 0; s < N; s += 2 * m) { for(size_t i = s, j = s + m; i < s + m; i++, j++) { T x = a[i], y = a[j]; a[i] = x + y; a[j] = (x - y) * w; } w *= idw[__builtin_ctz(++k)]; } } } static void transform(number_theoritic_transform& a, bool inverse = false) { size_type n = 0; while((1ul << n) < a.size()) n++; size_type N = 1 << n; a.resize(N); if(!inverse) { ntt(a); } else { intt(a); T inv = T(N).inverse(); for(size_t i = 0; i < a.size(); i++) a[i] *= inv; } } static number_theoritic_transform convolution(const number_theoritic_transform& ar, const number_theoritic_transform& br) { size_type size = ar.degree() + br.degree() + 1; number_theoritic_transform a(ar), b(br); a.resize(size); b.resize(size); transform(a, false); transform(b, false); for(size_t i = 0; i < a.size(); i++) a[i] *= b[i]; transform(a, true); a.resize(size); return a; } public: number_theoritic_transform convolute(const number_theoritic_transform& r) const { return (*this) * r; } number_theoritic_transform operator*(const number_theoritic_transform& r) const { return number_theoritic_transform(*this) *= r; } number_theoritic_transform& operator*=(const number_theoritic_transform& r) { return (*this) = convolution((*this), r); } number_theoritic_transform operator*(const_reference r) const { return number_theoritic_transform(*this) *= r; } number_theoritic_transform& operator*=(const_reference r) { for(size_t i = 0; i < this->size(); i++) (*this)[i] = (*this)[i] * r; return *this; } }; } // @docs docs/number_theoritic_transform.md #line 22 "000.cpp" using fps = cplib::formal_power_series>; using mint = fps::value_type; int main() { int x, y; scanf("%d%d", &x, &y); if(!(x < y)) std::swap(x, y); assert(x + y <= 100'000); std::vector fact(x + y + 1, 1), ifact(x + y + 1); for(int i = 1; i < fact.size(); i++) fact[i] = fact[i - 1] * i; for(int i = 0; i < fact.size(); i++) ifact[i] = fact[i].inverse(); mint ans = 0; for(int p = 0; p <= x; p++) { for(int q = 0; p + q <= x; q++) { const int A = x - (p + q); const int B = y - (p + q); const int C = p; ans += mint(2).power(C) * fact[A + B + C] * ifact[A] * ifact[B] * ifact[C]; } } printf("%lld\n", ans.value()); return 0; }