using i64 = long long;
using Int = int;
using Long = i64;
const Long P = 998244353;
const int LEN2 = 20;
using COMPLEX_TYPE = double;
struct Complex {
COMPLEX_TYPE x, y;
Complex(COMPLEX_TYPE x = 0.0, COMPLEX_TYPE y = 0.0) : x(x), y(y) {}
friend Complex operator+(const Complex &A, const Complex &B) {
return Complex(A.x + B.x, A.y + B.y);
}
friend Complex operator-(const Complex &A, const Complex &B) {
return Complex(A.x - B.x, A.y - B.y);
}
friend Complex operator*(const Complex &A, const Complex &B) {
return Complex(A.x * B.x - A.y * B.y, A.x * B.y + A.y * B.x);
}
friend Complex operator*(const Complex &A, const double &B) {
return Complex(A.x * B, A.y * B);
}
friend Complex operator/(const Complex &A, const double &B) {
return Complex(A.x / B, A.y / B);
}
Complex conj() const {
return Complex(x, -y);
}
}; // 复数类
const Long NTT_G = 3; // P 的原根
const Long NTT_I = 86583718; // 根号 P - 1
Long qpow(Long a, Long b) {
Long res = 1;
while (b) {
if (b & 1) res = res * a % P;
a = a * a % P;
b >>= 1;
}
return res;
}
Long qpow(Long a, Long b, const Long P) {
Long res = 1;
while (b) {
if (b & 1) res = res * a % P;
a = a * a % P;
b >>= 1;
}
return res;
}
Int add(Int a, const Int b) {
a += b;
if (a >= P) {
a -= P;
}
return a;
}
Int add(Int a, const Int b, const Int P) {
a += b;
if (a >= P) {
a -= P;
}
return a;
}
Int sub(Int a, const Int b) {
a -= b;
if (a < 0) {
a += P;
}
return a;
}
Int sub(Int a, const Int b, const Int P) {
a -= b;
if (a < 0) {
a += P;
}
return a;
}
template<typename T>
struct Poly : vector<T> {
Poly() = default;
Poly(int n) : vector<T>(n) {}
Poly(int n, T x) : vector<T>(n, x) {}
Poly(vector<T>::iterator s, vector<T>::iterator t) : vector<T>(s, t) {}
Poly(vector<T>::const_iterator s, vector<T>::const_iterator t) : vector<T>(s, t) {}
Poly(initializer_list<T> lst) : vector<T>(lst) {}
};
namespace Polynomial {
int norm(int x) {
return 1 << (__lg(x - 1) + 1);
}
template<typename T>
void norm(Poly<T> &a) {
assert(!a.empty());
a.resize(norm(a.size()));
}
template<typename T>
T eval(const Poly<T> &a, T x) {
x %= P;
int n = a.size();
T res = a[n - 1];
for (int i = n - 2 ; i >= 0 ; i--) {
res = (res * x + a[i]) % P;
}
return res;
}
template<typename T>
Poly<T> operator~(Poly<T> a) {
reverse(begin(a), end(a));
return a;
}
Poly<Int> operator*(Poly<Int> a, Int b) {
for (auto &x : a) {
x = Long(x) * b % P;
}
return a;
}
Poly<Int> operator-(Poly<Int> a, Poly<Int> b) {
int n = max(a.size(), b.size());
a.resize(n);
for (int i = 0 ; i < n ; i++) {
a[i] = sub(a[i], b[i]);
}
return a;
}
Poly<Int> operator-(Poly<Int> a) {
int n = a.size();
for (int i = 0 ; i < n ; i++) {
a[i] = sub(0, a[i]);
}
return a;
}
Poly<Int> operator+(Poly<Int> a, Poly<Int> b) {
int n = max(a.size(), b.size());
a.resize(n);
for (int i = 0 ; i < n ; i++) {
a[i] = add(a[i], b[i]);
}
return a;
}
template<typename T>
istream& operator>>(istream &in, Poly<T> &x) {
for (auto &x : x) {
in >> x;
}
return in;
}
template<typename T>
ostream& operator<<(ostream &out, const Poly<T> &x) {
for (auto &x : x) {
out << x << " ";
}
return out;
}
}
using namespace Polynomial;
namespace FFT {
int binary_upper_bound_in_bits(int w) {
if (w <= 0) return 1;
return 32 - __builtin_clz(w);
}
const COMPLEX_TYPE PI = acos((COMPLEX_TYPE) -1.0);
const COMPLEX_TYPE PI2 = PI / 2;
Complex I(0, 1); // 单位虚数
vector<Complex> r;
void DFT(Complex *a, int n) {
assert((n & (n - 1)) == 0);
if (r.empty()) {
r.resize(1 << LEN2);
for (int i = 0 ; i < LEN2 ; i++) {
int L = 1 << i;
r[L] = Complex{cos(PI2 / L), sin(PI2 / L)};
for (int j = L + 1 ; j < (L << 1) ; j++) {
r[j] = r[j - L] * r[L];
}
}
}
for (int i = n ; i >= 2 ; i /= 2) {
int L = i / 2;
for (int j = 0 ; j < L ; j++) {
Complex z = a[j + L];
a[j + L] = a[j] - z;
a[j] = a[j] + z;
}
for (int j = i, res = 1 ; j < n ; j += i, res++) {
for (int k = 0 ; k < L ; k++) {
Complex z = a[j + k + L] * r[res];
a[j + k + L] = a[j + k] - z;
a[j + k] = a[j + k] + z;
}
}
}
}
void IDFT(Complex *a, int n) {
assert((n & (n - 1)) == 0);
for (int i = 2 ; i <= n ; i *= 2) {
int L = i / 2;
for (int j = 0 ; j < L ; j++) {
Complex z = a[j + L];
a[j + L] = a[j] - z;
a[j] = a[j] + z;
}
for (int j = i, res = 1 ; j < n ; j += i, res++) {
for (int k = 0 ; k < L ; k++) {
Complex z = a[j + k + L];
a[j + k + L] = (a[j + k] - z) * r[res];
a[j + k] = a[j + k] + z;
}
}
}
for (int i = 0 ; i < n ; i++) {
a[i] = a[i] / n;
}
reverse(a + 1, a + n);
}
}
namespace Polynomial {
void DFT(Poly<Complex> &a) {
FFT::DFT(a.data(), a.size());
}
void IDFT(Poly<Complex> &a) {
FFT::IDFT(a.data(), a.size());
}
void dot(Poly<Complex> &a, Poly<Complex> &b) {
for (int i = 0 ; i < (int) a.size() ; i++) {
a[i] = a[i] * b[i];
}
}
Poly<Complex> operator*(Poly<Complex> a, Poly<Complex> b) {
int n = a.size() + b.size() - 1;
if (a.size() <= 16 || b.size() <= 16) {
Poly<Complex> c(n);
for (int i = 0 ; i < (int) a.size() ; i++) {
for (int j = 0; j < (int) b.size() ; j++) {
c[i + j] = c[i + j] + a[i] * b[j];
}
}
return c;
}
int L = norm(n);
a.resize(L); b.resize(L);
DFT(a); DFT(b); dot(a, b); IDFT(a);
return a.resize(n), a;
}
}
namespace NTT {
const Long g = NTT_G; // P 的原根
const Long I = NTT_I; // 根号 P - 1
Poly<Long> W;
Poly<Long> __inv;
Long inv(Long x) {
if (__inv.empty()) {
int len = 1 << LEN2;
__inv.resize(len);
__inv[0] = __inv[1] = 1;
for (int i = 2 ; i < len ; i++) {
__inv[i] = Long(P - P / i) * __inv[P % i] % P;
}
}
if (x < (int) __inv.size()) {
return __inv[x];
}
return qpow(x, P - 2);
}
Long inv(Long x, Long P) {
return qpow(x, P - 2, P);
}
void DIF(Int *a, int n) {
if (W.empty()) {
int L = 1 << LEN2;
W.resize(L);
Int wn = qpow(g, P / L);
W[L >> 1] = 1;
for (int i = L / 2 + 1 ; i < L ; i++) {
W[i] = Long(W[i - 1]) * wn % P;
}
for (int i = L / 2 - 1 ; i >= 1 ; i--) {
W[i] = W[i << 1];
}
}
for (int k = n >> 1 ; k ; k >>= 1) {
for (int i = 0 ; i < n ; i += k << 1) {
for (int j = 0 ; j < k ; j++) {
Int x = a[i + j], y = a[i + j + k];
a[i + j] = add(a[i + j], y);
a[i + j + k] = Long(sub(x, y)) * W[j + k] % P;
}
}
}
}
void IDIT(Int *a, int n) {
for (int k = 1 ; k < n ; k <<= 1) {
for (int i = 0 ; i < n ; i += k << 1) {
for (int j = 0 ; j < k ; j++) {
Int x = a[i + j], y = Long(a[i + j + k]) * W[j + k] % P;
a[i + j] = add(x, y);
a[i + j + k] = sub(x, y);
}
}
}
Int inv = P - (P - 1) / n;
for (int i = 0 ; i < n ; i++) {
a[i] = Long(a[i]) * inv % P;
}
reverse(a + 1, a + n);
}
}
namespace Polynomial {
void DIF(Poly<Int> &a) {
NTT::DIF(a.data(), a.size());
}
void IDIT(Poly<Int> &a) {
NTT::IDIT(a.data(), a.size());
}
void dot(Poly<Int> &a, Poly<Int> &b) {
for (int i = 0 ; i < (int) a.size() ; i++) {
a[i] = Long(a[i]) * b[i] % P;
}
}
Poly<Int> operator*(Poly<Int> a, Poly<Int> b) {
int n = a.size() + b.size() - 1;
if (a.size() <= 16 || b.size() <= 16) {
Poly<Int> c(n);
for (int i = 0 ; i < (int) a.size() ; i++) {
for (int j = 0; j < (int) b.size() ; j++) {
c[i + j] = add(c[i + j], Long(a[i]) * b[j] % P);
}
}
return c;
}
int L = norm(n);
a.resize(L); b.resize(L);
DIF(a); DIF(b); dot(a, b); IDIT(a);
return a.resize(n), a;
}
Poly<Int> mtt(const Poly<Int> &a, const Poly<Int> &b, const Long P) {
int n = a.size() + b.size() - 1;
if (a.size() <= 16 || b.size() <= 16) {
Poly<Int> c(n);
for (int i = 0 ; i < (int) a.size() ; i++) {
for (int j = 0; j < (int) b.size() ; j++) {
c[i + j] = add(c[i + j], Long(a[i]) * b[j] % P, P);
}
}
return c;
}
int L = norm(n);
Poly<Complex> x(L), y(L), z(L);
for (int i = 0 ; i < (int) a.size() ; i++) {
x[i] = Complex(a[i] & 65535, a[i] >> 16);
}
for (int i = 0 ; i < (int) b.size() ; i++) {
y[i] = Complex(b[i] & 65535, b[i] >> 16);
}
DFT(x); DFT(y);
z = x;
x[0] = x[0] * (y[0].x * 2);
z[0] = z[0] * (y[0].y * 2);
int w = FFT::binary_upper_bound_in_bits(L - 1);
for (int k = 0 ; k < w ; k++) {
int l = 1 << k, mask = l - 1;
for (int i = l ; i < (l << 1) ; i++) {
Complex p = y[i], q = y[i ^ mask].conj(), tx = (p + q), ty = (q - p) * FFT::I;
x[i] = x[i] * tx;
z[i] = z[i] * ty;
}
}
IDFT(x); IDFT(z);
Poly<Int> res(n);
auto get = [&](double x) {
Long v = x + .5;
return (__int128) ((v < 0) ? (v - 1) / 2 : (v + 1) / 2) % P;
};
for (int i = 0 ; i < n ; i++) {
res[i] = add((get(z[i].y) << (16 << 1)) % P, (get(z[i].x + x[i].y) << 16) % P, P);
res[i] = add(res[i], get(x[i].x), P);
}
return res;
}
Poly<Int> inv2k(Poly<Int> a) {
int n = a.size(), m = n >> 1;;
if (n == 1) {
return {Int(qpow(a[0], P - 2))};
}
auto b = inv2k(Poly<Int>(begin(a), begin(a) + m));
auto c = b;
b.resize(n);
DIF(a); DIF(b); dot(a, b); IDIT(a);
for (int i = 0 ; i < n ; i++) {
a[i] = (i < m ? 0 : sub(0, a[i]));
}
DIF(a); dot(a, b); IDIT(a);
return move(begin(c), end(c), begin(a)), a;
}
Poly<Int> inv2k(Poly<Int> a, const Long P) {
int n = a.size(), m = n >> 1;;
if (n == 1) {
return {Int(qpow(a[0], P - 2, P))};
}
auto b = inv2k(Poly<Int>(begin(a), begin(a) + m), P);
auto c = b;
b.resize(n);
a = mtt(a, b, P); a.resize(n);
for (int i = 0 ; i < n ; i++) {
a[i] = (i < m ? 0 : sub(0, a[i], P));
}
a = mtt(a, b, P); a.resize(n);
return move(begin(c), end(c), begin(a)), a;
}
Poly<Int> Inv(Poly<Int> a) {
int n = a.size();
norm(a);
a = inv2k(a);
return a.resize(n), a;
}
Poly<Int> Inv(Poly<Int> a, const Long P) {
int n = a.size();
norm(a);
a = inv2k(a, P);
return a.resize(n), a;
}
Poly<Int> deriv(Poly<Int> a) {
for (int i = 1 ; i < (int) a.size() ; i++) {
a[i - 1] = Long(a[i]) * i % P;
}
return a.pop_back(), a;
}
Poly<Int> deriv(Poly<Int> a, const Long P) {
for (int i = 1 ; i < (int) a.size() ; i++) {
a[i - 1] = Long(a[i]) * i % P;
}
return a.pop_back(), a;
}
Poly<Int> integ(Poly<Int> a) {
a.push_back(0);
for (int i = (int) a.size() - 1 ; i > 0 ; i--) {
a[i] = Long(a[i - 1]) * NTT::inv(i) % P;
}
return a[0] = 0, a;
}
Poly<Int> integ(Poly<Int> a, const Long P) {
a.push_back(0);
for (int i = (int) a.size() - 1 ; i > 0 ; i--) {
a[i] = Long(a[i - 1]) * NTT::inv(i, P) % P;
}
return a[0] = 0, a;
}
Poly<Int> Ln(Poly<Int> a) {
assert(a[0] == 1);
int n = a.size();
a = deriv(a) * Inv(a);
return a.resize(n - 1), integ(a);
}
Poly<Int> Ln(Poly<Int> a, const Long P) {
assert(a[0] == 1);
int n = a.size();
a = mtt(deriv(a, P), Inv(a, P), P);
return a.resize(n - 1), integ(a, P);
}
Poly<Int> Exp(Poly<Int> a) {
assert(a[0] == 0);
int n = a.size(), k = norm(n);
Poly<Int> b(1, 1), c, d;
a.resize(k);
for (int L = 2 ; L <= k ; L <<= 1) {
d = b; b.resize(L);
c = Ln(b); c.resize(L);
for (int i = 0 ; i < L ; i++) {
c[i] = sub(a[i], c[i]);
}
c[0] = add(c[0], 1);
DIF(b); DIF(c); dot(b, c); IDIT(b);
move(begin(d), end(d), begin(b));
}
return b.resize(n), b;
}
Poly<Int> Exp(Poly<Int> a, const Long P) {
assert(a[0] == 0);
int n = a.size(), k = norm(n);
Poly<Int> b(1, 1), c, d;
a.resize(k);
for (int L = 2 ; L <= k ; L <<= 1) {
d = b; b.resize(L);
c = Ln(b, P); c.resize(L);
for (int i = 0 ; i < L ; i++) {
c[i] = sub(a[i], c[i], P);
}
c[0] = add(c[0], 1, P);
b = mtt(b, c, P); b.resize(L);
move(begin(d), end(d), begin(b));
}
return b.resize(n), b;
}
void Pow(Poly<Int> &a, Long m0, Long m1, Long d) {
int n = a.size();
a.erase(begin(a), begin(a) + d);
int a0 = a[0]; assert(a0);
a = a * NTT::inv(a0);
a = Exp(Ln(a) * m0) * qpow(a0, m1);
a.resize(n); d *= m0;
for (int i = n - 1 ; i >= 0 ; i--) {
a[i] = (i >= d ? a[i - d] : 0);
}
}
void Pow(Poly<Int> &a, Long m0, Long m1, Long d, const Long P) {
int n = a.size();
a.erase(begin(a), begin(a) + d);
int a0 = a[0]; assert(a0);
for (auto &x : a) {
x = Long(x) * NTT::inv(a0, P) % P;
}
a = Ln(a, P);
for (auto &x : a) {
x = Long(x) * m0 % P;
}
a = Exp(a, P); a.resize(n);
for (auto &x : a) {
x = Long(x) * qpow(a0, m1, P) % P;
}
d *= m0;
for (int i = n - 1 ; i >= 0 ; i--) {
a[i] = (i >= d ? a[i - d] : 0);
}
}
Poly<Int> Pow(Poly<Int> a, Long m) {
int n = a.size(), d = 0;
while (d < n && a[d] == 0) {
d++;
}
if (d == n) return a;
if (Long(d) * m >= n) {
return Poly<Int>(n);
}
Pow(a, m, m, d);
return a;
}
Poly<Int> Pow(Poly<Int> a, Long m, const Long P) {
int n = a.size(), d = 0;
while (d < n && a[d] == 0) {
d++;
}
if (d == n) return a;
if (Long(d) * m >= n) {
return Poly<Int>(n);
}
Pow(a, m, m, d, P);
return a;
}
Poly<Int> Pow(Poly<Int> a, string str) {
int n = a.size(), d = 0;
while (d < n && a[d] == 0) {
d++;
}
if (d == n) return a;
if (d) {
if (str.size() > 8 || Long(d) * stoll(str) >= n) {
return Poly<Int>(n);
}
}
Long m0 = 0, m1 = 0;
const Long phi = P - 1;
for (int i = 0 ; i < (int) str.size() ; i++) {
m0 = (m0 * 10 + str[i] - '0') % P;
m1 = (m1 * 10 + str[i] - '0') % phi;
}
Pow(a, m0, m1, d);
return a;
}
Poly<Int> Pow(Poly<Int> a, string str, const Long P) {
int n = a.size(), d = 0;
while (d < n && a[d] == 0) {
d++;
}
if (d == n) return a;
if (d) {
if (str.size() > 8 || Long(d) * stoll(str) >= n) {
return Poly<Int>(n);
}
}
Long m0 = 0, m1 = 0;
const Long phi = P - 1;
for (int i = 0 ; i < (int) str.size() ; i++) {
m0 = (m0 * 10 + str[i] - '0') % P;
m1 = (m1 * 10 + str[i] - '0') % phi;
}
Pow(a, m0, m1, d, P);
return a;
}
pair<Poly<Int>, Poly<Int>> Div(Poly<Int> f, Poly<Int> g) {
int n = f.size() - 1, m = g.size() - 1;
if (n < m) {
return {Poly<Int>{1}, f};
}
Poly<Int> f_rev = ~f, g_rev = ~g;
g_rev.resize(n - m + 1);
Poly<Int> q = f_rev * Inv(g_rev);
q.resize(n - m + 1); q = ~q;
g = g * q;
f.resize(m); g.resize(m);
Poly<Int> r = f - g;
return {q, r};
}
pair<Poly<Int>, Poly<Int>> Div(Poly<Int> f, Poly<Int> g, const Long P) {
int n = f.size() - 1, m = g.size() - 1;
if (n < m) {
return {Poly<Int>{1}, f};
}
Poly<Int> f_rev = ~f, g_rev = ~g;
g_rev.resize(n - m + 1);
Poly<Int> q = mtt(f_rev, Inv(g_rev), P);
q.resize(n - m + 1); q = ~q;
g = g * q;
f.resize(m); g.resize(m);
Poly<Int> r = f;
for (int i = 0 ; i < m ; i++) {
r[i] = sub(r[i], g[i], P);
}
return {q, r};
}
Poly<Int> Sqrt(Poly<Int> a) {
Int w = 0; // x^2 - n
auto mul = [&](array<Long, 2> a, array<Long, 2> b) -> array<Long, 2> {
return {(a[0] * b[0] + a[1] * b[1] % P * w) % P, (a[0] * b[1] + a[1] * b[0]) % P};
};
auto qqpow = [&](array<Long, 2> a, Long b) {
array<Long, 2> res{1, 0};
while (b) {
if (b & 1) {
res = mul(res, a);
}
a = mul(a, a);
b >>= 1;
}
return res[0];
};
auto cipolla = [&](Int a) -> Int {
// x^2\equiv a\pmod p
if (qpow(a, (P - 1) / 2) == P - 1) {
// 无解
return -1;
}
Int x = 0;
while (true) {
x = rand() % P;
w = sub(Long(x) * x % P, a);
if (qpow(w, (P - 1) / 2) == P - 1) {
break;
}
}
Int Z = qqpow({Long(x), Long(1)}, (P + 1) / 2);
return min<Int>(Z, sub(0, Z));
};
int n = a.size();
Poly<Int> res{cipolla(a[0])};
for (int len = 1 ; len < (n << 1) ; len <<= 1) {
Poly<Int> b = a; b.resize(len);
res.resize(len);
Poly<Int> c = b * Inv(res); c.resize(len);
res = res + c;
res = res * NTT::inv(2);
}
return res.resize(n), res;
}
Poly<Int> Sqrt(Poly<Int> a, const Long P) {
Int w = 0; // x^2 - n
auto mul = [&](array<Long, 2> a, array<Long, 2> b) -> array<Long, 2> {
return {(a[0] * b[0] + a[1] * b[1] % P * w) % P, (a[0] * b[1] + a[1] * b[0]) % P};
};
auto qqpow = [&](array<Long, 2> a, Long b) {
array<Long, 2> res{1, 0};
while (b) {
if (b & 1) {
res = mul(res, a);
}
a = mul(a, a);
b >>= 1;
}
return res[0];
};
auto cipolla = [&](Int a) -> Int {
// x^2\equiv a\pmod p
if (qpow(a, (P - 1) / 2) == P - 1) {
// 无解
return -1;
}
Int x = 0;
while (true) {
x = rand() % P;
w = sub(Long(x) * x % P, a, P);
if (qpow(w, (P - 1) / 2, P) == P - 1) {
break;
}
}
Int Z = qqpow({Long(x), Long(1)}, (P + 1) / 2);
return min<Int>(Z, sub(0, Z, P));
};
int n = a.size();
Poly<Int> res{cipolla(a[0])};
for (int len = 1 ; len < (n << 1) ; len <<= 1) {
Poly<Int> b = a; b.resize(len);
res.resize(len);
Poly<Int> c = mtt(b, Inv(res), P); c.resize(len);
for (int i = 0 ; i < len ; i++) {
res[i] = add(res[i], c[i], P);
}
res = res * NTT::inv(2, P);
}
return res.resize(n), res;
}
Poly<Int> Sin(Poly<Int> &a) {
assert(a[0] == 0);
return (Exp(a * NTT::I) - Exp(a * sub(0, NTT::I))) * NTT::inv(2 * NTT::I % P);
}
Poly<Int> Sin(Poly<Int> &a, const Long I, const Long P) {
assert(a[0] == 0);
Poly<Int> b = a, c = a;
const Long II = sub(0, I, P);
const Long III = NTT::inv(2 * I % P, P);
for (int i = 0 ; i < (int) a.size() ; i++) {
b[i] = Long(b[i]) * I % P;
c[i] = Long(c[i]) * II % P;
}
b = Exp(b, P); c = Exp(c, P);
for (int i = 0 ; i < (int) a.size() ; i++) {
b[i] = sub(b[i], c[i], P);
b[i] = Long(b[i]) * III % P;
}
return b;
}
Poly<Int> Cos(Poly<Int> &a) {
assert(a[0] == 0);
return (Exp(a * NTT::I) + Exp(a * sub(0, NTT::I))) * NTT::inv(2);
}
Poly<Int> Cos(Poly<Int> &a, const Long I, const Long P) {
assert(a[0] == 0);
Poly<Int> b = a, c = a;
const Long II = sub(0, I, P);
const Long III = NTT::inv(2, P);
for (int i = 0 ; i < (int) a.size() ; i++) {
b[i] = Long(b[i]) * I % P;
c[i] = Long(c[i]) * II % P;
}
b = Exp(b, P); c = Exp(c, P);
for (int i = 0 ; i < (int) a.size() ; i++) {
b[i] = add(b[i], c[i], P);
b[i] = Long(b[i]) * III % P;
}
return b;
}
Poly<Int> Tan(Poly<Int> &a) {
assert(a[0] == 0);
Poly<Int> A = Exp(a * NTT::I), B = Exp(a * sub(0, NTT::I));
Poly<Int> C = (A - B) * NTT::inv(NTT::I) * Inv(A + B);
return C.resize(a.size()), C;
}
Poly<Int> Tan(Poly<Int> &a, const Long I, const Long P) {
assert(a[0] == 0);
Poly<Int> b = a, c = a;
const Long II = sub(0, I, P);
for (int i = 0 ; i < (int) a.size() ; i++) {
b[i] = Long(b[i]) * I % P;
c[i] = Long(c[i]) * II % P;
}
b = Exp(b, P); c = Exp(c, P);
Poly<Int> d = b, e = b;
const Long invI = NTT::inv(I, P);
for (int i = 0 ; i < (int) a.size() ; i++) {
d[i] = sub(d[i], c[i], P);
d[i] = Long(d[i]) * invI % P;
e[i] = add(e[i], c[i], P);
}
e = Inv(e, P);
d = mtt(d, e, P);
return d.resize(a.size()), d;
}
Poly<Int> Asin(Poly<Int> &a) {
assert(a[0] == 0);
int n = a.size();
Poly<Int> c = a * a; c.resize(n);
c[0] = sub(c[0], 1);
c = deriv(a) * Inv(Sqrt(-c)); c.resize(n - 1);
return integ(c);
}
Poly<Int> Asin(Poly<Int> &a, const Long P) {
assert(a[0] == 0);
int n = a.size();
Poly<Int> c = mtt(a, a, P); c.resize(n);
c[0] = sub(c[0], 1, P);
Poly<Int> d = c;
for (auto &x : d) {
x = sub(0, x, P);
}
c = deriv(a, P) * Inv(Sqrt(d, P), P); c.resize(n - 1);
return integ(c, P);
}
Poly<Int> Acos(Poly<Int> &a) {
assert(a[0] == 0);
int n = a.size();
Poly<Int> c = a * a; c.resize(n);
c[0] = sub(c[0], 1);
c = deriv(a) * Inv(Sqrt(-c)); c.resize(n - 1);
return integ(-c);
}
Poly<Int> Acos(Poly<Int> &a, const Long P) {
assert(a[0] == 0);
int n = a.size();
Poly<Int> c = mtt(a, a, P); c.resize(n);
c[0] = sub(c[0], 1, P);
Poly<Int> d = c;
for (auto &x : d) {
x = sub(0, x, P);
}
c = deriv(a) * Inv(Sqrt(d)); c.resize(n - 1);
for (auto &x : c) {
x = sub(0, x, P);
}
return integ(c, P);
}
Poly<Int> Atan(Poly<Int> &a) {
assert(a[0] == 0);
int n = a.size();
Poly<Int> c = a * a; c.resize(n);
c[0] = add(c[0], 1);
c = deriv(a) * Inv(c); c.resize(n - 1);
return integ(c);
}
Poly<Int> Atan(Poly<Int> &a, const Long P) {
assert(a[0] == 0);
int n = a.size();
Poly<Int> c = mtt(a, a, P); c.resize(n);
c[0] = add(c[0], 1, P);
c = deriv(a, P) * Inv(c, P); c.resize(n - 1);
return integ(c);
}
Poly<Int> BerlekampMassey(const Poly<Int> &a) {
// a 是 线性数列
// 返回线性递推式
int n = a.size();
Poly<Int> temp(n + 1), f(n + 1), g(n + 1);
int k = 0, last_k = 0, last_delta, last = -1;
for (int i = 0 ; i < n ; i++) {
Int delta = sub(0, a[i]);
for (int j = 1 ; j <= k ; j++) {
delta = add(delta, Long(f[j]) * a[i - j] % P);
}
if (delta == 0) continue;
if (last == -1) {
k = i + 1;
} else {
Int t = Long(delta) * qpow(last_delta, P - 2) % P;
temp[i - last] = add(temp[i - last], t);
for (int j = 1 ; j <= last_k ; j++) {
temp[i - last + j] = sub(temp[i - last + j], Long(t) * g[j] % P);
}
int p = last_k;
last_k = k;
k = max(k, i - last + p);
for (int j = 1 ; j <= last_k ; j++) {
g[j] = f[j];
}
for (int j = 1 ; j <= k ; j++) {
f[j] = temp[j];
}
}
last_delta = delta;
last = i;
}
return f.resize(k + 1), f;
}
Poly<Int> BerlekampMassey(const Poly<Int> &a, const Long P) {
// a 是 线性数列
// 返回线性递推式
int n = a.size();
Poly<Int> temp(n + 1), f(n + 1), g(n + 1);
int k = 0, last_k = 0, last_delta, last = -1;
for (int i = 0 ; i < n ; i++) {
Int delta = sub(0, a[i], P);
for (int j = 1 ; j <= k ; j++) {
delta = add(delta, Long(f[j]) * a[i - j] % P, P);
}
if (delta == 0) continue;
if (last == -1) {
k = i + 1;
} else {
Int t = Long(delta) * qpow(last_delta, P - 2, P) % P;
temp[i - last] = add(temp[i - last], t, P);
for (int j = 1 ; j <= last_k ; j++) {
temp[i - last + j] = sub(temp[i - last + j], Long(t) * g[j] % P, P);
}
int p = last_k;
last_k = k;
k = max(k, i - last + p);
for (int j = 1 ; j <= last_k ; j++) {
g[j] = f[j];
}
for (int j = 1 ; j <= k ; j++) {
f[j] = temp[j];
}
}
last_delta = delta;
last = i;
}
return f.resize(k + 1), f;
}
Int LinearRecurrence(const Poly<Int> &a, const Poly<Int> &r, i64 n) {
// // a 是线性数列
// // r 是线性递推式的系数, 可以通过 a 由 BerlekampMassey 得到
Poly<Int> F(begin(a), begin(a) + r.size() - 1), R = -r;
R[0] = 1;
F = F * R;
F.resize(R.size() - 1);
for ( ; n ; n /= 2) {
Poly<Int> rR = R;
for (int i = 0 ; i < (int) rR.size() ; i += 2) {
rR[i] = sub(0, rR[i]);
}
Poly<Int> FR = F * rR, RR = R * rR;
F.resize(0); R.resize(0);
for (int i = (n & 1) ; i < (int) FR.size() ; i+= 2) {
F.push_back(FR[i]);
}
for (int i = 0 ; i < (int) RR.size() ; i += 2) {
R.push_back(RR[i]);
}
}
return (F.empty() ? 0 : Long(F[0]) * qpow(R[0], P - 2) % P);
}
Int LinearRecurrence(const Poly<Int> &a, const Poly<Int> &r, i64 n, const Long P) {
// // a 是线性数列
// // r 是线性递推式的系数, 可以通过 a 由 BerlekampMassey 得到
Poly<Int> F(begin(a), begin(a) + r.size() - 1), R = r;
for (auto &x : R) {
x = sub(0, x, P);
}
R[0] = 1;
F = mtt(F, R, P);
F.resize(R.size() - 1);
for ( ; n ; n /= 2) {
Poly<Int> rR = R;
for (int i = 0 ; i < (int) rR.size() ; i += 2) {
rR[i] = sub(0, rR[i], P);
}
Poly<Int> FR = mtt(F, rR, P), RR = mtt(R, rR, P);
F.resize(0); R.resize(0);
for (int i = (n & 1) ; i < (int) FR.size() ; i+= 2) {
F.push_back(FR[i]);
}
for (int i = 0 ; i < (int) RR.size() ; i += 2) {
R.push_back(RR[i]);
}
}
return (F.empty() ? 0 : Long(F[0]) * qpow(R[0], P - 2, P) % P);
}
Poly<Int> ValueTrans(Poly<Int> a, int n, int m) {
vector<Int> fact(n + 1), invfact(n + 1);
fact[0] = 1;
for (int i = 1 ; i <= n ; i++) {
fact[i] = Long(fact[i - 1]) * i % P;
}
invfact[n] = qpow(fact[n], P - 2);
for (int i = n ; i >= 1 ; i--) {
invfact[i - 1] = Long(invfact[i]) * i % P;
}
vector<Int> down(n + 1);
down[0] = 1;
for (int i = m ; i >= m - n ; i--) {
down[0] = Long(down[0]) * i % P;
}
for (int i = 1 ; i <= n ; i++) {
down[i] = Long(down[i - 1]) * qpow(m - n + i - 1, P - 2) % P * add(m, i) % P;
}
for (int i = 0 ; i <= n ; i++) {
a[i] = Long(a[i]) * invfact[i] % P * invfact[n - i] % P;
if ((n ^ i) & 1) {
a[i] = sub(0, a[i]);
}
}
Poly<Int> b((n << 1) + 1);
for (int i = 0 ; i <= (n << 1) ; i++) {
b[i] = qpow(i + m - n, P - 2);
}
a = a * b;
Poly<Int> res(n + 1);
for (int i = 0 ; i <= n ; i++) {
res[i] = Long(down[i]) * a[i + n] % P;
}
return res;
}
Poly<Int> ValueTrans(Poly<Int> a, int n, int m, const Long P) {
vector<Int> fact(n + 1), invfact(n + 1);
fact[0] = 1;
for (int i = 1 ; i <= n ; i++) {
fact[i] = Long(fact[i - 1]) * i % P;
}
invfact[n] = qpow(fact[n], P - 2, P);
for (int i = n ; i >= 1 ; i--) {
invfact[i - 1] = Long(invfact[i]) * i % P;
}
vector<Int> down(n + 1);
down[0] = 1;
for (int i = m ; i >= m - n ; i--) {
down[0] = Long(down[0]) * i % P;
}
for (int i = 1 ; i <= n ; i++) {
down[i] = Long(down[i - 1]) * qpow(m - n + i - 1, P - 2, P) % P * add(m, i) % P;
}
for (int i = 0 ; i <= n ; i++) {
a[i] = Long(a[i]) * invfact[i] % P * invfact[n - i] % P;
if ((n ^ i) & 1) {
a[i] = sub(0, a[i]);
}
}
Poly<Int> b((n << 1) + 1);
for (int i = 0 ; i <= (n << 1) ; i++) {
b[i] = qpow(i + m - n, P - 2, P);
}
a = mtt(a, b, P);
Poly<Int> res(n + 1);
for (int i = 0 ; i <= n ; i++) {
res[i] = Long(down[i]) * a[i + n] % P;
}
return res;
}
}
namespace FWT {
void AND(Int *a, int n) {
for (int k = 2 ; k <= n ; k <<= 1) {
for (int i = 0 ; i < n ; i += k) {
int t = k >> 1;
for (int j = i ; j < i + t ; j++) {
a[j] = add(a[j], a[j + t]);
}
}
}
}
void IAND(Int *a, int n) {
for (int k = 2 ; k <= n ; k <<= 1) {
for (int i = 0 ; i < n ; i += k) {
int t = k >> 1;
for (int j = i ; j < i + t ; j++) {
a[j] = sub(a[j], a[j + t]);
}
}
}
}
void OR(Int *a, int n) {
for (int k = 2 ; k <= n ; k <<= 1) {
for (int i = 0 ; i < n ; i += k) {
int t = k >> 1;
for (int j = i ; j < i + t ; j++) {
a[j + t] = add(a[j + t], a[j]);
}
}
}
}
void IOR(Int *a, int n) {
for (int k = 2 ; k <= n ; k <<= 1) {
for (int i = 0 ; i < n ; i += k) {
int t = k >> 1;
for (int j = i ; j < i + t ; j++) {
a[j + t] = sub(a[j + t], a[j]);
}
}
}
}
void XOR(Int *a, int n) {
for (int k = 2 ; k <= n ; k <<= 1) {
for (int i = 0 ; i < n ; i += k) {
int t = k >> 1;
for (int j = i ; j < i + t ; j++) {
Int temp = a[j];
a[j] = add(a[j], a[j + t]);
a[j + t] = sub(temp, a[j + t]);
}
}
}
}
void IXOR(Int *a, int n) {
const Int inv2 = qpow(2, P - 2);
for (int k = 2 ; k <= n ; k <<= 1) {
for (int i = 0 ; i < n ; i += k) {
int t = k >> 1;
for (int j = i ; j < i + t ; j++) {
Int temp = a[j];
a[j] = add(a[j], a[j + t]);
a[j + t] = sub(temp, a[j + t]);
a[j] = Long(a[j]) * inv2 % P;
a[j + t] = Long(a[j + t]) * inv2 % P;
}
}
}
}
}
namespace Polynomial {
void AND(Poly<Int> &a) {
FWT::AND(a.data(), a.size());
}
void IAND(Poly<Int> &a) {
FWT::IAND(a.data(), a.size());
}
void OR(Poly<Int> &a) {
FWT::OR(a.data(), a.size());
}
void IOR(Poly<Int> &a) {
FWT::IOR(a.data(), a.size());
}
void XOR(Poly<Int> &a) {
FWT::XOR(a.data(), a.size());
}
void IXOR(Poly<Int> &a) {
FWT::IXOR(a.data(), a.size());
}
Poly<Int> operator&(Poly<Int> a, Poly<Int> b) {
int n = max(a.size(), b.size());
int L = norm(n);
if (a.size() <= 16 || b.size() <= 16) {
Poly<Int> c(L);
for (int i = 0 ; i < (int) a.size() ; i++) {
for (int j = 0; j < (int) b.size() ; j++) {
c[i & j] = add(c[i & j], Long(a[i]) * b[j] % P);
}
}
return c;
}
a.resize(L); b.resize(L);
AND(a); AND(b); dot(a, b); IAND(a);
return a.resize(n), a;
}
Poly<Int> operator|(Poly<Int> a, Poly<Int> b) {
int n = max(a.size(), b.size());
int L = norm(n);
if (a.size() <= 16 || b.size() <= 16) {
Poly<Int> c(L);
for (int i = 0 ; i < (int) a.size() ; i++) {
for (int j = 0; j < (int) b.size() ; j++) {
c[i | j] = add(c[i | j], Long(a[i]) * b[j] % P);
}
}
return c;
}
a.resize(L); b.resize(L);
OR(a); OR(b); dot(a, b); IOR(a);
return a.resize(n), a;
}
Poly<Int> operator^(Poly<Int> a, Poly<Int> b) {
int n = max(a.size(), b.size());
int L = norm(n);
if (a.size() <= 16 || b.size() <= 16) {
Poly<Int> c(L);
for (int i = 0 ; i < (int) a.size() ; i++) {
for (int j = 0; j < (int) b.size() ; j++) {
c[i ^ j] = add(c[i ^ j], Long(a[i]) * b[j] % P);
}
}
return c;
}
a.resize(L); b.resize(L);
XOR(a); XOR(b); dot(a, b); IXOR(a);
return a.resize(n), a;
}
}
namespace other {
Poly<Int> Stirling1stRow(int n) {
vector<int> fact(n + 1), invfact(n + 1);
fact[0] = 1;
for (int i = 1 ; i <= n ; i++) {
fact[i] = Long(fact[i - 1]) * i % P;
}
invfact[n] = qpow(fact[n], P - 2);
for (int i = n ; i >= 1 ; i--) {
invfact[i - 1] = Long(invfact[i]) * i % P;
}
function<Poly<Int>(int)> dfs = [&](int n) -> Poly<Int> {
if (n == 1) {
return {0, 1};
}
if (n & 1) {
Poly<Int> a = dfs(n - 1);
a.resize(n + 1);
for (int i = n ; i > 0 ; i--) {
a[i] = Long(a[i]) * sub(n, 1) % P;
a[i] = add(a[i], a[i - 1]);
}
a[0] = Long(a[0]) * sub(n, 1) % P;
return a;
} else {
Poly<Int> a = dfs(n / 2);
int len = n / 2;
a.resize(len + 1);
Poly<Int> b(len + 1), c(len + 1);
for (int i = 0 ; i <= len ; i++) {
b[i] = Long(a[len - i]) * fact[len - i] % P;
c[i] = Long(qpow(len, i)) * invfact[i] % P;
}
b = b * c;
for (int i = 0 ; i <= len ; i++) {
c[i] = Long(b[len - i]) * invfact[i] % P;
}
a = a * c;
return a;
}
};
return dfs(n);
}
Poly<Int> Stirling1Row(int n, const Long P) {
vector<int> fact(n + 1), invfact(n + 1);
fact[0] = 1;
for (int i = 1 ; i <= n ; i++) {
fact[i] = Long(fact[i - 1]) * i % P;
}
invfact[n] = qpow(fact[n], P - 2, P);
for (int i = n ; i >= 1 ; i--) {
invfact[i - 1] = Long(invfact[i]) * i % P;
}
function<Poly<Int>(int)> dfs = [&](int n) -> Poly<Int> {
if (n == 1) {
return {0, 1};
}
if (n & 1) {
Poly<Int> a = dfs(n - 1);
a.resize(n + 1);
for (int i = n ; i > 0 ; i--) {
a[i] = Long(a[i]) * sub(n, 1) % P;
a[i] = add(a[i], a[i - 1]);
}
a[0] = Long(a[0]) * sub(n, 1) % P;
return a;
} else {
Poly<Int> a = dfs(n / 2);
int len = n / 2;
a.resize(len + 1);
Poly<Int> b(len + 1), c(len + 1);
for (int i = 0 ; i <= len ; i++) {
b[i] = Long(a[len - i]) * fact[len - i] % P;
c[i] = Long(qpow(len, i, P)) * invfact[i] % P;
}
b = mtt(b, c, P);
for (int i = 0 ; i <= len ; i++) {
c[i] = Long(b[len - i]) * invfact[i] % P;
}
a = mtt(a, c, P);
return a;
}
};
return dfs(n);
}
Poly<Int> Stirling1Col(int n, int k) {
// 第 k 列的前 n 行
int mx = max(n, k);
vector<int> fact(mx + 1), invfact(mx + 1);
fact[0] = 1;
for (int i = 1 ; i <= mx ; i++) {
fact[i] = Long(fact[i - 1]) * i % P;
}
invfact[mx] = qpow(fact[mx], P - 2);
for (int i = mx ; i >= 1 ; i--) {
invfact[i - 1] = Long(invfact[i]) * i % P;
}
Poly<Int> a(n + 1);
for (int i = 1 ; i <= n ; i++) {
a[i] = Long(fact[i - 1]) * invfact[i] % P;
}
a = Pow(a, k);
for (int i = 0 ; i <= n ; i++) {
a[i] = Long(a[i]) * fact[i] % P * invfact[k] % P;
}
return a;
}
Poly<Int> Stirling1Col(int n, int k, const Long P) {
// 第 k 列的前 n 行
int mx = max(n, k);
vector<int> fact(mx + 1), invfact(mx + 1);
fact[0] = 1;
for (int i = 1 ; i <= mx ; i++) {
fact[i] = Long(fact[i - 1]) * i % P;
}
invfact[mx] = qpow(fact[mx], P - 2, P);
for (int i = mx ; i >= 1 ; i--) {
invfact[i - 1] = Long(invfact[i]) * i % P;
}
Poly<Int> a(n + 1);
for (int i = 1 ; i <= n ; i++) {
a[i] = Long(fact[i - 1]) * invfact[i] % P;
}
a = Pow(a, k, P);
for (int i = 0 ; i <= n ; i++) {
a[i] = Long(a[i]) * fact[i] % P * invfact[k] % P;
}
return a;
}
Poly<Int> Stirling2Row(int n) {
// 第 n 行
vector<int> fact(n + 1), invfact(n + 1);
fact[0] = 1;
for (int i = 1 ; i <= n ; i++) {
fact[i] = Long(fact[i - 1]) * i % P;
}
invfact[n] = qpow(fact[n], P - 2);
for (int i = n ; i >= 1 ; i--) {
invfact[i - 1] = Long(invfact[i]) * i % P;
}
Poly<Int> a(n + 1), b(n + 1);
for (int i = 0 ; i <= n ; i++) {
a[i] = Long(invfact[i]) * qpow(i, n) % P;
if (i & 1) {
b[i] = Long(invfact[i]) * (P - 1) % P;
} else {
b[i] = Long(invfact[i]) % P;
}
}
a = a * b;
return a.resize(n + 1), a;
}
Poly<Int> Stirling2Row(int n, const Long P) {
// 第 n 行
vector<int> fact(n + 1), invfact(n + 1);
fact[0] = 1;
for (int i = 1 ; i <= n ; i++) {
fact[i] = Long(fact[i - 1]) * i % P;
}
invfact[n] = qpow(fact[n], P - 2, P);
for (int i = n ; i >= 1 ; i--) {
invfact[i - 1] = Long(invfact[i]) * i % P;
}
Poly<Int> a(n + 1), b(n + 1);
for (int i = 0 ; i <= n ; i++) {
a[i] = Long(invfact[i]) * qpow(i, n, P) % P;
if (i & 1) {
b[i] = Long(invfact[i]) * (P - 1) % P;
} else {
b[i] = Long(invfact[i]) % P;
}
}
a = mtt(a, b, P);
return a.resize(n + 1), a;
}
Poly<Int> Stirling2Col(int n, int k) {
int mx = max(n, k);
vector<int> fact(mx + 1), invfact(mx + 1);
fact[0] = 1;
for (int i = 1 ; i <= mx ; i++) {
fact[i] = Long(fact[i - 1]) * i % P;
}
invfact[mx] = qpow(fact[mx], P - 2);
for (int i = mx ; i >= 1 ; i--) {
invfact[i - 1] = Long(invfact[i]) * i % P;
}
Poly<Int> a(n + 1);
for (int i = 1 ; i <= n ; i++) {
a[i] = invfact[i];
}
a = Pow(a, k);
for (int i = 0 ; i <= n ; i++) {
a[i] = Long(a[i]) * fact[i] % P * invfact[k] % P;
}
return a;
}
Poly<Int> Stirling2Col(int n, int k, const Long P) {
int mx = max(n, k);
vector<int> fact(mx + 1), invfact(mx + 1);
fact[0] = 1;
for (int i = 1 ; i <= mx ; i++) {
fact[i] = Long(fact[i - 1]) * i % P;
}
invfact[mx] = qpow(fact[mx], P - 2, P);
for (int i = mx ; i >= 1 ; i--) {
invfact[i - 1] = Long(invfact[i]) * i % P;
}
Poly<Int> a(n + 1);
for (int i = 1 ; i <= n ; i++) {
a[i] = invfact[i];
}
a = Pow(a, k, P);
for (int i = 0 ; i <= n ; i++) {
a[i] = Long(a[i]) * fact[i] % P * invfact[k] % P;
}
return a;
}
Poly<Int> Bell(int n) {
vector<int> fact(n + 1), invfact(n + 1);
fact[0] = 1;
for (int i = 1 ; i <= n ; i++) {
fact[i] = Long(fact[i - 1]) * i % P;
}
invfact[n] = qpow(fact[n], P - 2);
for (int i = n ; i >= 1 ; i--) {
invfact[i - 1] = Long(invfact[i]) * i % P;
}
Poly<Int> a(n + 1);
for (int i = 1 ; i <= n ; i++) {
a[i] = invfact[i];
}
a = Exp(a);
for (int i = 1 ; i <= n ; i++) {
a[i] = Long(a[i]) * fact[i] % P;
}
return a;
}
Poly<Int> Bell(int n, const Long P) {
vector<int> fact(n + 1), invfact(n + 1);
fact[0] = 1;
for (int i = 1 ; i <= n ; i++) {
fact[i] = Long(fact[i - 1]) * i % P;
}
invfact[n] = qpow(fact[n], P - 2, P);
for (int i = n ; i >= 1 ; i--) {
invfact[i - 1] = Long(invfact[i]) * i % P;
}
Poly<Int> a(n + 1);
for (int i = 1 ; i <= n ; i++) {
a[i] = invfact[i];
}
a = Exp(a, P);
for (int i = 1 ; i <= n ; i++) {
a[i] = Long(a[i]) * fact[i] % P;
}
return a;
}
Poly<Int> Partition(int n, int m) {
assert(n >= m);
Poly<Int> a(n + 1);
for (int i = 1 ; i <= m ; i++) {
for (int j = i ; j <= n ; j += i) {
a[j] = add(a[j], sub(0, NTT::inv(j / i)));
}
}
a = Inv(Exp(a));
return a;
}
Poly<Int> Partition(int n, int m, const Long P) {
assert(n >= m);
Poly<Int> a(n + 1);
for (int i = 1 ; i <= m ; i++) {
for (int j = i ; j <= n ; j += i) {
a[j] = add(a[j], sub(0, NTT::inv(j / i, P)), P);
}
}
a = Inv(Exp(a, P), P);
return a;
}
}