#include // created [2020/02/05] 11:35:57 #pragma GCC diagnostic ignored "-Wsign-compare" #pragma GCC diagnostic ignored "-Wsign-conversion" using i32 = int32_t; using i64 = int64_t; using u32 = uint32_t; using u64 = uint64_t; using uint = unsigned int; using usize = std::size_t; using ll = long long; using ull = unsigned long long; using ld = long double; template using arr = T (&)[n]; template using c_arr = const T (&)[n]; template using max_heap = std::priority_queue; template using min_heap = std::priority_queue, std::greater>; template constexpr T popcount(const T u) { return u ? static_cast(__builtin_popcountll(static_cast(u))) : static_cast(0); } template constexpr T log2p1(const T u) { return u ? static_cast(64 - __builtin_clzll(static_cast(u))) : static_cast(0); } template constexpr T msbp1(const T u) { return log2p1(u); } template constexpr T lsbp1(const T u) { return __builtin_ffsll(u); } template constexpr T clog(const T u) { return u ? log2p1(u - 1) : static_cast(u); } template constexpr bool ispow2(const T u) { return u and (static_cast(u) & static_cast(u - 1)) == 0; } template constexpr T ceil2(const T u) { return static_cast(1) << clog(u); } template constexpr T floor2(const T u) { return u == 0 ? static_cast(0) : static_cast(1) << (log2p1(u) - 1); } template constexpr bool btest(const T mask, const usize ind) { return static_cast((static_cast(mask) >> ind) & static_cast(1)); } template void bset(T& mask, const usize ind) { mask |= (static_cast(1) << ind); } template void breset(T& mask, const usize ind) { mask &= ~(static_cast(1) << ind); } template void bflip(T& mask, const usize ind) { mask ^= (static_cast(1) << ind); } template void bset(T& mask, const usize ind, const bool b) { (b ? bset(mask, ind) : breset(mask, ind)); } template constexpr T bcut(const T mask, const usize ind) { return ind == 0 ? static_cast(0) : static_cast((static_cast(mask) << (64 - ind)) >> (64 - ind)); } template bool chmin(T& a, const T& b) { return (a > b ? a = b, true : false); } template bool chmax(T& a, const T& b) { return (a < b ? a = b, true : false); } constexpr unsigned int mod = 1000000007; template constexpr T inf_v = std::numeric_limits::max() / 4; template constexpr Real pi_v = Real{3.141592653589793238462643383279502884}; auto mfp = [](auto&& f) { return [=](auto&&... args) { return f(f, std::forward(args)...); }; }; template T in() { T v; return std::cin >> v, v; } template T in_v(typename std::enable_if<(i == n), c_arr>::type) { return in(); } template auto in_v(typename std::enable_if<(i < n), c_arr>::type& szs) { const usize s = (usize)szs[i]; std::vector(szs))> ans(s); for (usize j = 0; j < s; j++) { ans[j] = in_v(szs); } return ans; } template auto in_v(c_arr szs) { return in_v(szs); } template auto in_t() { return std::tuple...>{in()...}; } struct io_init { io_init() { std::cin.tie(nullptr), std::ios::sync_with_stdio(false); std::cout << std::fixed << std::setprecision(20); } void clear() { std::cin.tie(), std::ios::sync_with_stdio(true); } } io_setting; int out() { return 0; } template int out(const T& v) { return std::cout << v, 0; } template int out(const std::vector& v) { for (usize i = 0; i < v.size(); i++) { if (i > 0) { std::cout << ' '; } out(v[i]); } return 0; } template int out(const std::pair& v) { return out(v.first), std::cout << ' ', out(v.second), 0; } template int out(const T& v, const Args... args) { return out(v), std::cout << ' ', out(args...), 0; } template int outln(const Args... args) { return out(args...), std::cout << '\n', 0; } template int outel(const Args... args) { return out(args...), std::cout << std::endl, 0; } # define SHOW(...) static_cast(0) constexpr ull TEN(const usize n) { return n == 0 ? 1ULL : TEN(n - 1) * 10ULL; } template auto make_v(typename std::enable_if<(i == n), c_arr>::type, const T& v = T{}) { return v; } template auto make_v(typename std::enable_if<(i < n), c_arr>::type szs, const T& v = T{}) { const usize s = (usize)szs[i]; return std::vector(szs, v))>(s, make_v(szs, v)); } template auto make_v(c_arr szs, const T& t = T{}) { return make_v(szs, t); } template struct complex { using value_type = Real; complex() : real{Real{0}}, imag{Real{0}} {} complex(const complex&) = default; complex(const Real& theta) : real(std::cos(theta)), imag(std::sin(theta)) {} complex(const Real& r, const Real& i) : real{r}, imag{i} {} ~complex() = default; friend complex operator+(const complex& c) { return c; } friend complex operator-(const complex& c) { return complex{-c.real, -c.imag}; } friend complex operator+(const complex& c1, const complex& c2) { return complex{c1.real + c2.real, c1.imag + c2.imag}; } friend complex operator-(const complex& c1, const complex& c2) { return complex{c1.real - c2.real, c1.imag - c2.imag}; } friend complex operator*(const complex& c1, const complex& c2) { return complex{c1.real * c2.real - c1.imag * c2.imag, c1.real * c2.imag + c1.imag * c2.real}; } friend complex operator*(const complex& c, const Real& r) { return complex{c.real * r, c.imag * r}; } friend complex operator/(complex& c1, complex& c2) { c1* c2.conj() / c2.norm(); } friend bool operator==(const complex& c1, const complex& c2) { return c1.real == c2.real and c1.imag == c2.imag; } friend bool operator!=(const complex& c1, const complex& c2) { return not(c1 == c2); } friend complex& operator+=(complex& c1, const complex& c2) { return c1.real += c2.real, c1.imag += c2.imag, c1; } friend complex& operator-=(complex& c1, const complex& c2) { return c1.real += c2.real, c1.imag += c2.imag, c1; } friend complex& operator*=(complex& c1, const complex& c2) { return c1 = c1 * c2; } friend complex& operator*=(complex& c, const Real& r) { return c = c * r; } friend complex& operator/=(complex& c1, const complex& c2) { return c1 = c1 / c2; } complex conj() const { return complex{real, -imag}; } Real norm() const { return real * real + imag * imag; } Real abs() const { return std::sqrt(norm()); } Real arg() const { return std::atan2(imag, real); } friend std::ostream& operator<<(std::ostream& os, const complex& c) { return os << c.real << "+" << c.imag << "i"; } Real real, imag; }; template T gcd(const T& a, const T& b) { return a < 0 ? gcd(-a, b) : b < 0 ? gcd(a, -b) : (a > b ? gcd(b, a) : a == 0 ? b : gcd(b % a, a)); } template T lcm(const T& a, const T& b) { return a / gcd(a, b) * b; } template constexpr std::pair extgcd(const T a, const T b) { if (b == 0) { return std::pair{1, 0}; } const auto g = gcd(a, b), da = std::abs(b) / g; const auto p = extgcd(b, a % b); const auto x = (da + p.second % da) % da, y = (g - a * x) / b; return {x, y}; } template constexpr T inverse(const T a, const T mod) { return extgcd(a, mod).first; } template class modint_base { public: template static std::enable_if_t mod() { return mod_ref(); } template static constexpr std::enable_if_t mod() { return mod_value; } template static void set_mod(const std::enable_if_t mod) { mod_ref() = mod, inv_ref() = {1, 1}; } modint_base() : v{0} {} modint_base(const ll val) : v{norm(static_cast(val % static_cast(mod()) + static_cast(mod())))} {} modint_base(const modint_base& n) : v{n()} {} explicit operator bool() const { return v != 0; } bool operator!() const { return not static_cast(*this); } modint_base& operator=(const modint_base& m) { return v = m(), (*this); } modint_base& operator=(const ll val) { return v = norm(uint(val % static_cast(mod()) + static_cast(mod()))), (*this); } friend modint_base operator+(const modint_base& m) { return m; } friend modint_base operator-(const modint_base& m) { return make(norm(mod() - m.v)); } friend modint_base operator+(const modint_base& m1, const modint_base& m2) { return make(norm(m1.v + m2.v)); } friend modint_base operator-(const modint_base& m1, const modint_base& m2) { return make(norm(m1.v + mod() - m2.v)); } friend modint_base operator*(const modint_base& m1, const modint_base& m2) { return make(static_cast(static_cast(m1.v) * static_cast(m2.v) % static_cast(mod()))); } friend modint_base operator/(const modint_base& m1, const modint_base& m2) { return m1 * inv(m2.v); } friend modint_base operator+(const modint_base& m, const ll val) { return modint_base{static_cast(m.v) + val}; } friend modint_base operator-(const modint_base& m, const ll val) { return modint_base{static_cast(m.v) - val}; } friend modint_base operator*(const modint_base& m, const ll val) { return modint_base{static_cast(m.v) * (val % static_cast(mod()))}; } friend modint_base operator/(const modint_base& m, const ll val) { return modint_base{static_cast(m.v) * inv(val)}; } friend modint_base operator+(const ll val, const modint_base& m) { return modint_base{static_cast(m.v) + val}; } friend modint_base operator-(const ll val, const modint_base& m) { return modint_base{-static_cast(m.v) + val}; } friend modint_base operator*(const ll val, const modint_base& m) { return modint_base{static_cast(m.v) * (val % static_cast(mod()))}; } friend modint_base operator/(const ll val, const modint_base& m) { return modint_base{val * inv(static_cast(m.v))}; } friend modint_base& operator+=(modint_base& m1, const modint_base& m2) { return m1 = m1 + m2; } friend modint_base& operator-=(modint_base& m1, const modint_base& m2) { return m1 = m1 - m2; } friend modint_base& operator*=(modint_base& m1, const modint_base& m2) { return m1 = m1 * m2; } friend modint_base& operator/=(modint_base& m1, const modint_base& m2) { return m1 = m1 / m2; } friend modint_base& operator+=(modint_base& m, const ll val) { return m = m + val; } friend modint_base& operator-=(modint_base& m, const ll val) { return m = m - val; } friend modint_base& operator*=(modint_base& m, const ll val) { return m = m * val; } friend modint_base& operator/=(modint_base& m, const ll val) { return m = m / val; } friend modint_base operator^(const modint_base& m, const ll n) { return power(m.v, n); } friend modint_base& operator^=(modint_base& m, const ll n) { return m = m ^ n; } friend bool operator==(const modint_base& m1, const modint_base& m2) { return m1.v == m2.v; } friend bool operator!=(const modint_base& m1, const modint_base& m2) { return not(m1 == m2); } friend bool operator==(const modint_base& m, const ll val) { return m.v == norm(static_cast(static_cast(mod()) + val % static_cast(mod()))); } friend bool operator!=(const modint_base& m, const ll val) { return not(m == val); } friend bool operator==(const ll val, const modint_base& m) { return m.v == norm(static_cast(static_cast(mod()) + val % static_cast(mod()))); } friend bool operator!=(const ll val, const modint_base& m) { return not(m == val); } friend std::istream& operator>>(std::istream& is, modint_base& m) { ll v; return is >> v, m = v, is; } friend std::ostream& operator<<(std::ostream& os, const modint_base& m) { return os << m(); } uint operator()() const { return v; } static modint_base small_inv(const usize n) { auto& in = inv_ref(); if (n < in.size()) { return in[n]; } for (usize i = in.size(); i <= n; i++) { in.push_back(-in[modint_base::mod() % i] * (modint_base::mod() / i)); } return in.back(); } std::pair quad() const { const auto ans = quad_r(v, mod()); ll x = std::get<0>(ans), y = std::get<1>(ans); if (y < 0) { x = -x, y = -y; } return {x, y}; } private: static std::tuple quad_r(const ll r, const ll p) // r = x/y (mod p), (x,y,z) s.t. x=yr+pz { if (std::abs(r) <= 1000) { return {r, 1, 0}; } ll nr = p % r, q = p / r; if (nr * 2LL >= r) { nr -= r, q++; } if (nr * 2LL <= -r) { nr += r, q--; } const auto sub = quad_r(nr, r); const ll x = std::get<0>(sub), z = std::get<1>(sub), y = std::get<2>(sub); return {x, y - q * z, z}; } template static std::enable_if_t mod_ref() { static UInt mod = 0; return mod; } static uint norm(const uint x) { return x < mod() ? x : x - mod(); } static modint_base make(const uint x) { modint_base m; return m.v = x, m; } static modint_base power(modint_base x, ull n) { modint_base ans = 1; for (; n; n >>= 1, x *= x) { if (n & 1) { ans *= x; } } return ans; } static modint_base inv(const ll v) { return v <= 2000000 ? small_inv(static_cast(v)) : modint_base{inverse(v, static_cast(mod()))}; } static std::vector& inv_ref() { static std::vector in{1, 1}; return in; } uint v; }; template using modint = modint_base; template using dynamic_modint = modint_base; template class fft { private: static constexpr usize depth = 30; static constexpr Real pi = pi_v; static void transform(std::vector>& a, const usize lg, const bool rev) { static std::vector> root[depth]; const usize sz = a.size(); assert((1UL << lg) == sz); if (root[lg].empty()) { root[lg].reserve(sz), root[lg].resize(sz); for (usize i = 0; i < sz; i++) { root[lg][i] = complex(pi * Real(2 * i) / Real(sz)); } } std::vector> tmp(sz); for (usize w = (sz >> 1); w > 0; w >>= 1) { for (usize y = 0; y < (sz >> 1); y += w) { const complex r = rev ? root[lg][y].conj() : root[lg][y]; for (usize x = 0; x < w; x++) { const auto u = a[y << 1 | x], v = a[y << 1 | x | w] * r; tmp[y | x] = u + v, tmp[y | x | (sz >> 1)] = u - v; } } std::swap(tmp, a); } } public: using value_type = Real; fft() = delete; template static std::vector simple_convolute(const std::vector& a, const std::vector& b) { const usize need = a.size() + b.size() - 1, lg = clog(need), sz = 1UL << lg; std::vector> x(sz), y(sz); for (usize i = 0; i < a.size(); i++) { x[i] = {(Real)a[i], (Real)0}; } for (usize i = 0; i < b.size(); i++) { y[i] = {(Real)b[i], (Real)0}; } transform(x, lg, false), transform(y, lg, false); for (usize i = 0; i < sz; i++) { x[i] *= y[i]; } transform(x, lg, true); std::vector ans(need); for (usize i = 0; i < need; i++) { ans[i] = (T)std::round(x[i].real / (Real)sz); } return ans; } template static std::vector convolute(const std::vector& a, const std::vector& b) { constexpr usize bitnum = (depth + division - 1) / division; const usize need = a.size() + b.size() - 1, lg = clog(need), sz = 1UL << lg; std::vector> x[division], y[division], tmp(sz); for (usize i = 0; i < division; i++) { x[i].reserve(sz), x[i].resize(sz), y[i].reserve(sz), y[i].resize(sz); std::fill(tmp.begin() + std::min(a.size(), b.size()), tmp.end(), complex{}); for (usize j = 0; j < a.size(); j++) { tmp[j].real = value_type((a[j] >> (bitnum * i)) & ((1 << bitnum) - 1)); } for (usize j = 0; j < b.size(); j++) { tmp[j].imag = value_type((b[j] >> (bitnum * i)) & ((1 << bitnum) - 1)); } transform(tmp, lg, false); for (usize j = 0; j < sz; j++) { tmp[j] *= value_type(0.5); } for (usize j = 0; j < sz; j++) { const usize k = j == 0 ? 0UL : sz - j; x[i][j] = complex{tmp[j].real + tmp[k].real, tmp[j].imag - tmp[k].imag}, y[i][j] = complex{tmp[j].imag + tmp[k].imag, -tmp[j].real + tmp[k].real}; } } std::vector> z[division]; for (usize i = 0; i < division; i++) { z[i].reserve(sz), z[i].resize(sz); } for (usize a = 0; a < division; a++) { for (usize b = 0; b < division; b++) { for (usize i = 0; i < sz; i++) { if (a + b < division) { z[a + b][i] += x[a][i] * y[b][i]; } else { z[a + b - division][i] += x[a][i] * y[b][i] * complex(0, 1); } } } } for (usize i = 0; i < division; i++) { transform(z[i], lg, true); } std::vector ans(need); T base = 1; for (usize k = 0; k < 2 * division - 1; k++, base *= (1LL << bitnum)) { for (usize i = 0; i < need; i++) { if (k < division) { ans[i] += base * T(std::round(z[k][i].real / value_type(sz))); } else { ans[i] += base * T(std::round(z[k - division][i].imag / value_type(sz))); } } } return ans; } template static std::vector> convolute(const std::vector>& a, const std::vector>& b) { using mint = modint_base; constexpr usize bitnum = (depth + division - 1) / division; const usize need = a.size() + b.size() - 1, lg = clog(need), sz = 1UL << lg; std::vector> x[division], y[division], tmp(sz); for (usize i = 0; i < division; i++) { x[i].reserve(sz), x[i].resize(sz), y[i].reserve(sz), y[i].resize(sz); std::fill(tmp.begin() + std::min(a.size(), b.size()), tmp.end(), complex{}); for (usize j = 0; j < a.size(); j++) { tmp[j].real = value_type((a[j]() >> (bitnum * i)) & ((1 << bitnum) - 1)); } for (usize j = 0; j < b.size(); j++) { tmp[j].imag = value_type((b[j]() >> (bitnum * i)) & ((1 << bitnum) - 1)); } transform(tmp, lg, false); for (usize j = 0; j < sz; j++) { tmp[j] *= value_type(0.5); } for (usize j = 0; j < sz; j++) { const usize k = j == 0 ? 0UL : sz - j; x[i][j] = complex{tmp[j].real + tmp[k].real, tmp[j].imag - tmp[k].imag}, y[i][j] = complex{tmp[j].imag + tmp[k].imag, -tmp[j].real + tmp[k].real}; } } std::vector> z[division]; for (usize i = 0; i < division; i++) { z[i].reserve(sz), z[i].resize(sz); } for (usize a = 0; a < division; a++) { for (usize b = 0; b < division; b++) { for (usize i = 0; i < sz; i++) { if (a + b < division) { z[a + b][i] += x[a][i] * y[b][i]; } else { z[a + b - division][i] += x[a][i] * y[b][i] * complex(0, 1); } } } } for (usize i = 0; i < division; i++) { transform(z[i], lg, true); } std::vector ans(need); mint base = 1; for (usize k = 0; k < 2 * division - 1; k++, base *= (1LL << bitnum)) { for (usize i = 0; i < need; i++) { if (k < division) { ans[i] += int((base * ll(std::round(z[k][i].real / value_type(sz))))()); } else { ans[i] += int((base * ll(std::round(z[k - division][i].imag / value_type(sz))))()); } } } return ans; } }; template class bigint { public: using bucket_type = ll; static constexpr bucket_type bucket_value = TEN(bucket_digits); bigint() {} bigint(const bucket_type x) : sign(x >= 0), d(x == 0 ? 0 : 1, x < 0 ? -x : x) { normalize(); } bigint(const bool sign, const std::vector& d) : sign(sign), d(d) { normalize(); } bigint(const bigint& n) : sign(n.sign), d(n.d) {} bigint(const std::string& str) { assert(not str.empty()); sign = str[0] != '-'; for (usize i = 0; i < str.size(); i += bucket_digits) { const usize r = str.size() - i, l = (sign ? (r >= bucket_digits ? r - bucket_digits : 0) : (r > bucket_digits ? r - bucket_digits : 1)); if (r > l) { d.push_back(static_cast(std::stoll(str.substr(l, r - l)))); } } } friend bigint operator+(const bigint& m) { return m; } friend bigint operator-(const bigint& m) { return bigint{not m.sign, m.d}; } explicit operator bool() const { return not d.empty(); } bool operator!() const { return not static_cast(*this); } friend bigint operator+(const bigint& m, const bigint& n) { if (m.sign == n.sign) { usize sz = std::max(m.d.size(), n.d.size()); std::vector v(sz, 0); for (usize i = 0; i < m.d.size(); i++) { v[i] += m.d[i]; } for (usize i = 0; i < n.d.size(); i++) { v[i] += n.d[i]; } return bigint{m.sign, v}; } else { usize sz = std::max(m.d.size(), n.d.size()); std::vector v(sz, 0); if (abs_comp(m, n) == -1) { for (usize i = 0; i < m.d.size(); i++) { v[i] -= m.d[i]; } for (usize i = 0; i < n.d.size(); i++) { v[i] += n.d[i]; } return bigint{n.sign, v}; } else { for (usize i = 0; i < m.d.size(); i++) { v[i] += m.d[i]; } for (usize i = 0; i < n.d.size(); i++) { v[i] -= n.d[i]; } return bigint{m.sign, v}; } } } friend bigint operator-(const bigint& m, const bigint& n) { return m + (-n); } friend bigint operator*(const bigint& m, const bigint& n) { return bigint{not(m.sign ^ n.sign), fft<>::convolute(m.d, n.d)}; } friend bigint operator/(const bigint& m, const bigint& n) { return div(m, n).first; } friend bigint operator%(const bigint& m, const bigint& n) { return div(m, n).second; } friend bigint operator^(const bigint& m, const bigint& n) { return n == 0 ? bigint{1} : n % 2 == 1 ? (m ^ (n - 1)) * m : (m * m) ^ (n / 2); } friend bool operator<(const bigint& m, const bigint& n) { return comp(m, n) == -1; } friend bool operator>(const bigint& m, const bigint& n) { return comp(m, n) == 1; } friend bool operator==(const bigint& m, const bigint& n) { return comp(m, n) == 0; } friend bool operator!=(const bigint& m, const bigint& n) { return not(m == n); } friend bool operator<=(const bigint& m, const bigint& n) { return not(m > n); } friend bool operator>=(const bigint& m, const bigint& n) { return not(m < n); } friend bigint operator<<(const bigint& m, const usize n) // *(B^n) { std::vector ans(m.d.size() + n, 0); for (usize i = 0; i < m.d.size(); i++) { ans[i + n] = m.d[i]; } return bigint{m.sign, ans}; } friend bigint operator>>(const bigint& m, const usize n) // /(B^n) { if (m.d.size() <= n) { return 0; } std::vector ans(m.d.size() - n, 0); for (usize i = 0; i < m.d.size() - n; i++) { ans[i] = m.d[i + n]; } return bigint{m.sign, ans}; } bigint& operator=(const bigint& n) { return sign = n.sign, d = n.d, *this; } bigint& operator=(const bucket_type n) { return *this = bigint{n}; } friend bigint& operator+=(bigint& m, const bigint& n) { return m = m + n; } friend bigint& operator-=(bigint& m, const bigint& n) { return m = m - n; } friend bigint& operator*=(bigint& m, const bigint& n) { return m = m * n; } friend bigint& operator/=(bigint& m, const bigint& n) { return m = m / n; } friend bigint& operator%=(bigint& m, const bigint& n) { return m = m % n; } friend bigint& operator^=(bigint m, const bigint& n) { return m = m ^ n; } friend bigint& operator<<=(bigint& m, const usize n) { return m = m << n; } friend bigint& operator>>=(bigint& m, const usize n) { return m = m >> n; } friend bigint operator+(const bigint& m, const bucket_type n) { return m + bigint{n}; } friend bigint operator-(const bigint& m, const bucket_type n) { return m - bigint{n}; } friend bigint operator*(const bigint& m, const bucket_type n) { return small_multiply(m, n); } friend bigint operator/(const bigint& m, const bucket_type n) { return m / bigint{n}; } friend bigint operator%(const bigint& m, const bucket_type n) { return m % bigint{n}; } friend bool operator<(const bigint& m, const bucket_type n) { return comp(m, bigint{n}) == -1; } friend bool operator>(const bigint& m, const bucket_type n) { return comp(m, bigint{n}) == 1; } friend bool operator==(const bigint& m, const bucket_type n) { return comp(m, bigint{n}) == 0; } friend bool operator!=(const bigint& m, const bucket_type n) { return not(m == n); } friend bool operator<=(const bigint& m, const bucket_type n) { return not(m > n); } friend bool operator>=(const bigint& m, const bucket_type n) { return not(m < n); } friend bigint& operator+=(const bigint& m, const bucket_type n) { return m += bigint{n}; } friend bigint& operator-=(const bigint& m, const bucket_type n) { return m -= bigint{n}; } friend bigint& operator*=(const bigint& m, const bucket_type n) { return m *= bigint{n}; } friend bigint& operator/=(const bigint& m, const bucket_type n) { return m /= bigint{n}; } friend bigint& operator%=(const bigint& m, const bucket_type n) { return m %= bigint{n}; } friend bigint operator+(const bucket_type m, const bigint& n) { return bigint{m} + n; } friend bigint operator-(const bucket_type m, const bigint& n) { return bigint{m} - n; } friend bigint operator*(const bucket_type m, const bigint& n) { return bigint{m} * n; } friend bigint operator/(const bucket_type m, const bigint& n) { return bigint{m} / n; } friend bigint operator%(const bucket_type m, const bigint& n) { return bigint{m} % n; } friend bool operator<(const bucket_type m, const bigint& n) { return comp(bigint{m}, n) == -1; } friend bool operator>(const bucket_type m, const bigint& n) { return comp(bigint{m}, n) == 1; } friend bool operator==(const bucket_type m, const bigint& n) { return comp(bigint{m}, n) == 0; } friend bool operator!=(const bucket_type m, const bigint& n) { return not(m == n); } friend bool operator<=(const bucket_type m, const bigint& n) { return not(m > n); } friend bool operator>=(const bucket_type m, const bigint& n) { return not(m < n); } bigint& operator++() { return *this += bigint{1}; } bigint operator++(int) { return std::exchange(*this, *this + bigint{1}); } bigint& operator--() { return *this -= bigint{1}; } bigint operator--(int) { return std::exchange(*this, *this - bigint{1}); } friend std::pair div(bigint m, bigint n) { assert(not n.d.empty()); if (abs_comp(m, n) == -1) { return {0, m}; } const bool ms = m.sign, ns = n.sign; m.sign = n.sign = true; const usize sho_size = m.d.size() - n.d.size(); std::vector sho(sho_size + 1); for (usize i = 0; i <= sho_size; i++) { const usize ind = sho_size - i; const auto L = m >> ind; bucket_type inf = -1, sup = bucket_value; while (sup - inf > 1) { const bucket_type mid = (sup + inf) >> 1; (abs_comp(L, small_multiply(n, mid)) == -1 ? sup : inf) = mid; } sho[ind] = inf, m -= (small_multiply(n, inf) << ind); } return m.sign = ms, std::pair{bigint{not(ms ^ ns), sho}, m}; } friend std::istream& operator>>(std::istream& is, bigint& n) { std::string str; is >> str; return n = bigint{str}, is; } friend std::ostream& operator<<(std::ostream& os, const bigint& n) { if (n.d.empty()) { return os << '0'; } if (n.sign == false) { os << '-'; } const usize sz = n.d.size(); for (usize i = 0; i < sz; i++) { if (i == 0) { os << n.d[sz - i - 1]; } else { usize nz = 0; for (bucket_type base = 1; n.d[sz - i - 1] >= base; nz++, base *= bucket_type(10)) {} for (usize i = 0; i + nz < bucket_digits; i++) { os << '0'; } if (n.d[sz - i - 1] != 0LL) { os << n.d[sz - i - 1]; } } } return os; } friend usize log10p1(const bigint& m) { if (m.d.empty()) { return 0; } assert(m.sign); return (m.d.size() - 1) * bucket_digits + log10p1(m.d.back()); } friend bigint abs(const bigint& m) { auto ans = m; return ans.sign = true, ans; } long long to_ll() const { long long ans = 0, base = 1; for (usize i = 0; i < d.size(); i++, base *= bucket_value) { ans += base * d[i]; } return sign ? ans : -ans; } std::string to_string() const { if (d.empty()) { return "0"; } std::string ans; if (sign == false) { ans.push_back('-'); } const usize N = d.size(); for (usize i = 0; i < N; i++) { if (i == 0) { ans += std::to_string(d[N - i - 1]); } else { usize nz = 0; for (bucket_type U = 1; d[N - i - 1] >= U; nz++, U *= bucket_type(10)) {} for (usize i = 0; i + nz < bucket_digits; i++) { ans.push_back('0'); } if (d[N - i - 1] != 0LL) { ans += std::to_string(d[N - i - 1]); } } } return ans; } private: friend bigint small_multiply(const bigint& m, bucket_type n) // [TODO] karatsuba { const bool s = n < 0 ? not m.sign : m.sign; n = n < 0 ? -n : n; assert(n < bucket_value); std::vector ans(m.d.size()); for (usize i = 0; i < m.d.size(); i++) { ans[i] = m.d[i] * n; } return bigint{s, ans}; } friend int comp(const bigint& m, const bigint& n) { if (m.sign != n.sign) { return m.sign ? 1 : -1; } else { return m.sign ? abs_comp(m, n) : -abs_comp(m, n); } } static usize log10p1(const bucket_type n) { usize ans = 0; for (bucket_type base = 1; base <= n; base *= static_cast(10), ans++) {} return ans; } friend int abs_comp(const bigint& m, const bigint& n) // mn:1 { if (m.d.size() != n.d.size()) { return m.d.size() < n.d.size() ? -1 : 1; } else { for (usize i = 0; i < m.d.size(); i++) { const auto m_back = m.d[m.d.size() - i - 1], n_back = n.d[n.d.size() - i - 1]; if (m_back != n_back) { return m_back < n_back ? -1 : 1; } } } return 0; } void normalize() { for (usize i = 0; i < d.size(); i++) { if (d[i] >= bucket_value) { if (i + 1 == d.size()) { d.push_back(d[i] / bucket_value); } else { d[i + 1] += d[i] / bucket_value; } d[i] %= bucket_value; } else if (d[i] < 0) { assert(i + 1 != d.size()); while (d[i] < 0) { d[i + 1]--, d[i] += bucket_value; } } } for (; not d.empty() and d.back() == 0; d.pop_back()) {} } bool sign = true; std::vector d; }; int main() { const auto N = in>(); auto pow = mfp([&](auto&& self, const bigint<>& x, const bigint<>& n) -> bigint<> { if (n == 0) { return 1; } if (n % 2 == 1) { return self(self, x, n - 1) * x % N; } else { return self(self, x * x % N, n / 2); } }); if (N % 2 == 0) { return outel('!', N / 2, 2); } for (bigint a = 2; a * a <= N; a++) { if (N % a == 0) { return outel('!', N / a, a); } outel('?', a); const auto p = in>(); if (p % 2 == 1) { continue; } const auto b = pow(a, p / 2) + 1; if (b != N and b != 1 and N % b == 0) { return outel('!', N / b, b); } } return 0; }