···cpp
include<bits/stdc++.h>
typedef std::vector
const int Mo = 998244353;
inline void add(int &x, int y) { (x += y) >= Mo && (x -= Mo); }
inline void dec(int &x, int y) { (x -= y) < 0 && (x += Mo); }
inline int add(int x) { return x >= Mo ? x - Mo : x; }
inline int dec(int x) { return x < 0 ? x + Mo : x; }
inline void Mod(int &x) { x >= Mo && (x -= Mo), x < 0 && (x += Mo); }
inline int mul(int x, int y, int Mo = Mo) {return 1ll * x * y % Mo; }
inline int pow(int x, int k, int Mo = Mo) {
int ans = 1;
for (; k; k >>= 1, x = 1ll * x * x % Mo)
if (k & 1)
ans = 1ll * ans * x % Mo;
return ans;
}
inline int inv(int x, int Mo = Mo) { return pow(x, Mo - 2); }
inline int BSGS(int a, int b, int p) {
std::unordered_map<int, int> mp;
int m = ceil(sqrt(p));
for (int i = 0; i <= m; b = 1ll * b * a % p, i++)
mp[b] = i;
a = pow(a, m);
for (int i = 0, j = 1; i < m; j = 1ll * j * a % p, i++)
if (mp.count(j) && i * m >= mp[j])
return i * m - mp[j];
return -1;
}
inline int root(int p) {
Vec fac;
int x = p - 1;
for (int i = 2; i * i <= x; i++)
if (x % i == 0) {
fac.emplace_back(i);
while (x % i == 0) x /= i;
}
if (x > 1) fac.emplace_back(x);
for (int i = 2; i <= (p - 1); i++) {
bool tf = true;
for (auto j : fac)
if (pow(i, (p - 1) / j, p) == 1) {
tf = false;
break;
}
if (tf) return i;
}
return -1;
}
inline int degree(int a, int b, int p) {
int y = BSGS(3, a, p);
assert(y >= 0 && y % b == 0);
int z = pow(3, y / b);
return std::min(z, p - z);
}
namespace FFT {
int extend(int x);
void NTT(Vec &A, bool x);
void DFT(Vec &A);
void IDFT(Vec &A);
int extend(int x) {
int n = 1;
for (; n < x; n <<= 1) ;
return n;
}
void NTT(Vec &A, bool x) {
int n = A.size(), k = 0;
for (; (1 << k) < n; k++);
Vec rev(n);
for (int i = 0; i < n; i++)
rev[i] = rev[i >> 1] >> 1 | (i & 1) << (k - 1);
for (int i = 0; i < n; i++)
if (i < rev[i])
std::swap(A[i], A[rev[i]]);
for (int i = 2; i <= n; i <<= 1) {
int i1 = i >> 1, w = pow(3, (Mo - 1) / i);
if (x)
w = inv(w);
for (int j = 0; j < n; j += i) {
int x1, x2, x3 = 1;
for (int k = 0; k < i1; ++k, x3 = 1LL * x3 * w % Mo) {
x1 = A[j + k];
x2 = 1LL * x3 * A[j + k + i1] % Mo;
A[j + k] = (x1 + x2) % Mo;
A[j + k + i1] = (x1 - x2 + Mo) % Mo;
}
}
}
}
void DFT(Vec &A) {
NTT(A, false);
}
void IDFT(Vec &A) {
NTT(A, true);
int t = inv(A.size());
for (auto &x : A) x = 1ll * x * t % Mo;
}
}
using namespace FFT;
namespace Poly {
Vec operator +(Vec a, Vec b);
Vec operator +(Vec a, int b);
Vec operator +(int b, Vec a);
Vec operator -(Vec a, Vec b);
Vec operator -(Vec a, int b);
Vec operator -(int b, Vec a);
Vec operator -(Vec a);
Vec operator *(Vec a, Vec b);
Vec operator *(Vec a, int b);
Vec operator *(int b, Vec a);
Vec operator /(Vec a, Vec b);
Vec operator /(Vec a, int b);
Vec operator %(Vec a, Vec b);
Vec operator ~(Vec a);
Vec operator ^(Vec a, int b);
Vec operator <<(Vec a, int b);
Vec operator >>(Vec a, int b);
Vec fix(Vec a, int b);
Vec der(Vec a, bool tf);
Vec Inte(Vec a, bool tf);
Vec Sqrt(Vec a);
Vec root(Vec a, int b);
Vec Ln(Vec a);
Vec Exp(Vec a);
Vec sin(Vec a);
Vec cos(Vec a);
Vec tan(Vec a);
Vec asin(Vec a);
Vec acos(Vec a);
Vec atan(Vec a);
void print(Vec a);
Vec operator +(Vec a, Vec b) {
int n = std::max(a.size(), b.size());
a.resize(n), b.resize(n);
for (int i = 0; i < n; i++) add(a[i], b[i]);
return a;
}
Vec operator +(Vec a, int b) {
add(a[0], b);
return a;
}
Vec operator +(int b, Vec a) {
add(a[0], b);
return a;
}
Vec operator -(Vec a, Vec b) {
int n = std::max(a.size(), b.size());
a.resize(n), b.resize(n);
for (int i = 0; i < n; i++) dec(a[i], b[i]);
return a;
}
Vec operator -(Vec a, int b) {
dec(a[0], b);
return a;
}
Vec operator -(int b, Vec a) {
a = -a;
add(a[0], b);
return a;
}
Vec operator -(Vec a) {
for (auto &x : a)
x = dec(-x);
return a;
}
Vec operator *(Vec a, Vec b) {
int n = a.size() + b.size() - 1;
int n1 = extend(n);
a.resize(n1), DFT(a), b.resize(n1), DFT(b);
for (int i = 0; i < n1; i++)
a[i] = 1ll * a[i] * b[i] % Mo;
IDFT(a), a.resize(n);
return a;
}
Vec operator *(Vec a, int b) {
for (auto &x : a) x = 1ll * x * b % Mo;
return a;
}
Vec operator *(int b, Vec a) {
for (auto &x : a) x = 1ll * x * b % Mo;
return a;
}
Vec operator /(Vec a, Vec b) {
int n = a.size() - b.size() + 1;
if (n <= 0)
return Vec(1, 0);
std::reverse(a.begin(), a.end());
std::reverse(b.begin(), b.end());
a.resize(n), b.resize(n);
a = fix(a * ~b, n);
std::reverse(a.begin(), a.end());
return a;
}
Vec operator /(Vec a, int b) {
return a * inv(b);
}
Vec operator %(Vec a, Vec b) {
int n = b.size() - 1;
return fix(a - a / b * b, n);
}
Vec operator ~(Vec a) {
int n = a.size();
int n1 = extend(n), x1 = 0;
a.resize(n1);
Vec b(n1, 0);
b[0] = inv(a[0]);
for (int i = 2; i <= n1; i <<= 1) {
Vec x(i), y(i);
std::copy(a.begin(), a.begin() + i, x.begin());
std::copy(b.begin(), b.begin() + i, y.begin());
x1 = i << 1;
x.resize(x1), DFT(x);
y.resize(x1), DFT(y);
for (int ii = 0; ii < x1; ii++)
x[ii] = 1ll * y[ii] * (2 - 1ll * x[ii] * y[ii] % Mo + Mo) % Mo;
IDFT(x);
x.resize(i);
std::copy(x.begin(), x.begin() + i, b.begin());
}
b.resize(n);
return b;
}
Vec operator ^(Vec a, int b) {
int n = a.size();
int x = 0, y;
while (x < n && a[x] == 0)
x++;
if (1ll * x * b >= n)
return Vec(n, 0);
a = a >> x, y = a[0];
return (Exp(Ln(a) * b) * pow(y, b)) << (x * b);
}
Vec operator <<(Vec a, int b) {
int n = a.size();
Vec B(n, 0);
for (int i = 0; i < n - b; i++)
B[i + b] = a[i];
return B;
}
Vec operator >>(Vec a, int b) {
int n = a.size();
Vec B(n, 0);
for (int i = 0; i < n - b; i++)
B[i] = a[i + b];
return B;
}
Vec fix(Vec a, int b) {
a.resize(b);
return a;
}
Vec der(Vec a, bool tf = true) {
int n = a.size();
if (n == 1)
return Vec(1, 0);
Vec b(n - 1, 0);
for (int i = 1; i < n; i++)
b[i - 1] = 1ll * i * a[i] % Mo;
if (tf)
b.resize(n);
return b;
}
Vec Inte(Vec a, bool tf = true) {
int n = a.size();
// if (n == 1)
// return Vec(1, 0);
Vec b(n + 1, 0);
for (int i = 1; i <= n; i++)
b[i] = 1ll * inv(i) * a[i - 1] % Mo;
if (tf)
b.resize(n);
return b;
}
Vec Sqrt(Vec a) {
int n = a.size();
int n1 = extend(n);
a.resize(n1);
Vec b(n1, 0);
b[0] = degree(a[0], 2, Mo);
int t = inv(2);
for (int i = 2; i <= n1; i <<= 1) {
Vec x(i), y(i);
std::copy(a.begin(), a.begin() + i, x.begin());
std::copy(b.begin(), b.begin() + i, y.begin());
Vec z = ~y;
int x1 = i << 1;
x.resize(x1), DFT(x);
y.resize(x1), DFT(y);
z.resize(x1), DFT(z);
for (int ii = 0; ii < x1; ii++)
y[ii] = 1ll * y[ii] * y[ii] % Mo, x[ii] = 1ll * (x[ii] + y[ii]) * t % Mo * z[ii] % Mo;
IDFT(x);
x.resize(i);
std::copy(x.begin(), x.begin() + i, b.begin());
}
b.resize(n);
return b;
}
Vec root(Vec a, int b) { return b == 1 ? a : b == 2 ? Sqrt(a) : Exp(Ln(a) / b); }
Vec Ln(Vec a) {
assert(a[0] == 1);
int n = a.size();
return Inte(fix(der(a) * ~a, n));
}
Vec Exp(Vec a) {
assert(a[0] == 0);
int n = a.size();
int n1 = extend(n);
a.resize(n1);
Vec b(n1, 0);
b[0] = 1;
for (int i = 2; i <= n1; i <<= 1) {
Vec x = (-Ln(fix(b, i)) + fix(a, i) + 1) * fix(b, i);
std::copy(x.begin(), x.begin() + i, b.begin());
}
b.resize(n);
return b;
}
Vec sin(Vec a) {
assert(a[0] == 0);
int x = degree(Mo - 1, 2, Mo);
Vec b = Exp(x * a);
return (b - ~b) / ((long long)2 * x % Mo);
}
Vec cos(Vec a) {
assert(a[0] == 0);
int x = degree(Mo - 1, 2, Mo);
Vec b = Exp(x * a);
return (b + ~b) / 2;
}
Vec tan(Vec a) {
assert(a[0] == 0);
// int x = degree(Mo - 1, 2, Mo);
// Vec b = Exp(x * a);
int n = a.size();
return fix(sin(a) * ~cos(a), n);
}
Vec asin(Vec a) {
assert(a[0] == 0);
int n = a.size();
return Inte(fix(der(a) * ~Sqrt(1 - fix(a * a, n)), n));
}
Vec acos(Vec a) {
assert(a[0] == 0);
return -asin(a);
}
Vec atan(Vec a) {
assert(a[0] == 0);
int n = a.size();
return Inte(fix(der(a) * ~(1 + fix(a * a , n)), n));
}
void print(Vec a) {
int n = a.size();
for (int i = 0; i < n - 1; i++) printf("%d ", a[i]);
printf("%d
", a[n - 1]);
}
}
using namespace Poly;
···