#938363 (C++17) No.2587 Random Walk on Tree

提出ソース
結果

問題	No.2587 Random Walk on Tree
ユーザー	NyaanNyaan
提出日時	2023-12-15 02:35:22
言語	C++17 (gcc 12.3.0 + boost 1.83.0)
結果	TLE
実行時間	-
コード長	63,998 bytes
コンパイル時間	5,858 ms
コンパイル使用メモリ	339,668 KB
実行使用メモリ	46,444 KB
最終ジャッジ日時	2023-12-23 23:37:24
合計ジャッジ時間	23,339 ms
ジャッジサーバーID （参考情報）	judge15 / judge11
このコードへのチャレンジ
（要ログイン）
テストケース

テストケース表示
入力	結果	実行時間実行使用メモリ
testcase_00	AC	3 ms 6,548 KB
testcase_01	AC	2 ms 6,548 KB
testcase_02	AC	2 ms 6,548 KB
testcase_03	AC	2 ms 6,548 KB
testcase_04	AC	3 ms 6,548 KB
testcase_05	AC	4 ms 6,548 KB
testcase_06	AC	4 ms 6,548 KB
testcase_07	AC	3 ms 6,548 KB
testcase_08	AC	4 ms 6,548 KB
testcase_09	AC	2 ms 6,548 KB
testcase_10	AC	4 ms 6,548 KB
testcase_11	AC	7 ms 6,548 KB
testcase_12	AC	26 ms 6,548 KB
testcase_13	AC	37 ms 6,604 KB
testcase_14	AC	13 ms 6,548 KB
testcase_15	AC	2,053 ms 39,140 KB
testcase_16	AC	1,162 ms 24,696 KB
testcase_17	AC	1,934 ms 24,580 KB
testcase_18	AC	193 ms 9,240 KB
testcase_19	TLE	-
testcase_20	--	-
testcase_21	--	-
testcase_22	--	-
testcase_23	--	-
testcase_24	--	-
testcase_25	--	-
testcase_26	--	-
testcase_27	--	-
testcase_28	--	-
testcase_29	--	-
testcase_30	--	-
testcase_31	--	-
testcase_32	--	-
testcase_33	--	-
testcase_34	--	-
testcase_35	--	-
testcase_36	--	-
testcase_37	--	-
testcase_38	--	-
testcase_39	--	-
権限があれば一括ダウンロードができます
ソースコード

raw source code




#define NDEBUG

using namespace std;


#include <immintrin.h>

#include <algorithm>
#include <array>
#include <bitset>
#include <cassert>
#include <cctype>
#include <cfenv>
#include <cfloat>
#include <chrono>
#include <cinttypes>
#include <climits>
#include <cmath>
#include <complex>
#include <cstdarg>
#include <cstddef>
#include <cstdint>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <deque>
#include <fstream>
#include <functional>
#include <initializer_list>
#include <iomanip>
#include <ios>
#include <iostream>
#include <istream>
#include <iterator>
#include <limits>
#include <list>
#include <map>
#include <memory>
#include <new>
#include <numeric>
#include <ostream>
#include <queue>
#include <random>
#include <set>
#include <sstream>
#include <stack>
#include <streambuf>
#include <string>
#include <tuple>
#include <type_traits>
#include <typeinfo>
#include <unordered_map>
#include <unordered_set>
#include <utility>
#include <vector>



namespace Nyaan {
using ll = long long;
using i64 = long long;
using u64 = unsigned long long;
using i128 = __int128_t;
using u128 = __uint128_t;

template <typename T>
using V = vector<T>;
template <typename T>
using VV = vector<vector<T>>;
using vi = vector<int>;
using vl = vector<long long>;
using vd = V<double>;
using vs = V<string>;
using vvi = vector<vector<int>>;
using vvl = vector<vector<long long>>;
template <typename T>
using minpq = priority_queue<T, vector<T>, greater<T>>;

template <typename T, typename U>
struct P : pair<T, U> {
 template <typename... Args>
 P(Args... args) : pair<T, U>(args...) {}

 using pair<T, U>::first;
 using pair<T, U>::second;

 P &operator+=(const P &r) {
 first += r.first;
 second += r.second;
 return *this;
 }
 P &operator-=(const P &r) {
 first -= r.first;
 second -= r.second;
 return *this;
 }
 P &operator*=(const P &r) {
 first *= r.first;
 second *= r.second;
 return *this;
 }
 template <typename S>
 P &operator*=(const S &r) {
 first *= r, second *= r;
 return *this;
 }
 P operator+(const P &r) const { return P(*this) += r; }
 P operator-(const P &r) const { return P(*this) -= r; }
 P operator*(const P &r) const { return P(*this) *= r; }
 template <typename S>
 P operator*(const S &r) const {
 return P(*this) *= r;
 }
 P operator-() const { return P{-first, -second}; }
};

using pl = P<ll, ll>;
using pi = P<int, int>;
using vp = V<pl>;

constexpr int inf = 1001001001;
constexpr long long infLL = 4004004004004004004LL;

template <typename T>
int sz(const T &t) {
 return t.size();
}

template <typename T, typename U>
inline bool amin(T &x, U y) {
 return (y < x) ? (x = y, true) : false;
}
template <typename T, typename U>
inline bool amax(T &x, U y) {
 return (x < y) ? (x = y, true) : false;
}

template <typename T>
inline T Max(const vector<T> &v) {
 return *max_element(begin(v), end(v));
}
template <typename T>
inline T Min(const vector<T> &v) {
 return *min_element(begin(v), end(v));
}
template <typename T>
inline long long Sum(const vector<T> &v) {
 return accumulate(begin(v), end(v), 0LL);
}

template <typename T>
int lb(const vector<T> &v, const T &a) {
 return lower_bound(begin(v), end(v), a) - begin(v);
}
template <typename T>
int ub(const vector<T> &v, const T &a) {
 return upper_bound(begin(v), end(v), a) - begin(v);
}

constexpr long long TEN(int n) {
 long long ret = 1, x = 10;
 for (; n; x *= x, n >>= 1) ret *= (n & 1 ? x : 1);
 return ret;
}

template <typename T, typename U>
pair<T, U> mkp(const T &t, const U &u) {
 return make_pair(t, u);
}

template <typename T>
vector<T> mkrui(const vector<T> &v, bool rev = false) {
 vector<T> ret(v.size() + 1);
 if (rev) {
 for (int i = int(v.size()) - 1; i >= 0; i--) ret[i] = v[i] + ret[i + 1];
 } else {
 for (int i = 0; i < int(v.size()); i++) ret[i + 1] = ret[i] + v[i];
 }
 return ret;
};

template <typename T>
vector<T> mkuni(const vector<T> &v) {
 vector<T> ret(v);
 sort(ret.begin(), ret.end());
 ret.erase(unique(ret.begin(), ret.end()), ret.end());
 return ret;
}

template <typename F>
vector<int> mkord(int N, F f) {
 vector<int> ord(N);
 iota(begin(ord), end(ord), 0);
 sort(begin(ord), end(ord), f);
 return ord;
}

template <typename T>
vector<int> mkinv(vector<T> &v) {
 int max_val = *max_element(begin(v), end(v));
 vector<int> inv(max_val + 1, -1);
 for (int i = 0; i < (int)v.size(); i++) inv[v[i]] = i;
 return inv;
}

vector<int> mkiota(int n) {
 vector<int> ret(n);
 iota(begin(ret), end(ret), 0);
 return ret;
}

template <typename T>
T mkrev(const T &v) {
 T w{v};
 reverse(begin(w), end(w));
 return w;
}

template <typename T>
bool nxp(vector<T> &v) {
 return next_permutation(begin(v), end(v));
}



template <typename T>
vector<vector<T>> product(const vector<T> &a) {
 vector<vector<T>> ret;
 vector<T> v;
 auto dfs = [&](auto rc, int i) -> void {
 if (i == (int)a.size()) {
 ret.push_back(v);
 return;
 }
 for (int j = 0; j < a[i]; j++) v.push_back(j), rc(rc, i + 1), v.pop_back();
 };
 dfs(dfs, 0);
 return ret;
}



template <typename T>
T Power(T a, long long n, const T &I, const function<void(T &)> &f) {
 T res = I;
 for (; n; f(a = a * a), n >>= 1) {
 if (n & 1) f(res = res * a);
 }
 return res;
}

template <typename T>
T Power(T a, long long n, const T &I = T{1}) {
 return Power(a, n, I, function<void(T &)>{[](T &) -> void {}});
}

template <typename T>
T Rev(const T &v) {
 T res = v;
 reverse(begin(res), end(res));
 return res;
}

template <typename T>
vector<T> Transpose(const vector<T> &v) {
 using U = typename T::value_type;
 int H = v.size(), W = v[0].size();
 vector res(W, T(H, U{}));
 for (int i = 0; i < H; i++) {
 for (int j = 0; j < W; j++) {
 res[j][i] = v[i][j];
 }
 }
 return res;
}

template <typename T>
vector<T> Rotate(const vector<T> &v, int clockwise = true) {
 using U = typename T::value_type;
 int H = v.size(), W = v[0].size();
 vector res(W, T(H, U{}));
 for (int i = 0; i < H; i++) {
 for (int j = 0; j < W; j++) {
 if (clockwise) {
 res[W - 1 - j][i] = v[i][j];
 } else {
 res[j][H - 1 - i] = v[i][j];
 }
 }
 }
 return res;
}

} 




namespace Nyaan {
__attribute__((target("popcnt"))) inline int popcnt(const u64 &a) {
 return _mm_popcnt_u64(a);
}
inline int lsb(const u64 &a) { return a ? __builtin_ctzll(a) : 64; }
inline int ctz(const u64 &a) { return a ? __builtin_ctzll(a) : 64; }
inline int msb(const u64 &a) { return a ? 63 - __builtin_clzll(a) : -1; }
template <typename T>
inline int gbit(const T &a, int i) {
 return (a >> i) & 1;
}
template <typename T>
inline void sbit(T &a, int i, bool b) {
 if (gbit(a, i) != b) a ^= T(1) << i;
}
constexpr long long PW(int n) { return 1LL << n; }
constexpr long long MSK(int n) { return (1LL << n) - 1; }
} 




namespace Nyaan {

template <typename T, typename U>
ostream &operator<<(ostream &os, const pair<T, U> &p) {
 os << p.first << " " << p.second;
 return os;
}
template <typename T, typename U>
istream &operator>>(istream &is, pair<T, U> &p) {
 is >> p.first >> p.second;
 return is;
}

template <typename T>
ostream &operator<<(ostream &os, const vector<T> &v) {
 int s = (int)v.size();
 for (int i = 0; i < s; i++) os << (i ? " " : "") << v[i];
 return os;
}
template <typename T>
istream &operator>>(istream &is, vector<T> &v) {
 for (auto &x : v) is >> x;
 return is;
}

void in() {}
template <typename T, class... U>
void in(T &t, U &...u) {
 cin >> t;
 in(u...);
}

void out() { cout << "\n"; }
template <typename T, class... U, char sep = ' '>
void out(const T &t, const U &...u) {
 cout << t;
 if (sizeof...(u)) cout << sep;
 out(u...);
}

struct IoSetupNya {
 IoSetupNya() {
 cin.tie(nullptr);
 ios::sync_with_stdio(false);
 cout << fixed << setprecision(15);
 cerr << fixed << setprecision(7);
 }
} iosetupnya;

} 





#ifdef NyaanDebug
#define trc(...) (void(0))
#else
#define trc(...) (void(0))
#endif

#ifdef NyaanLocal
#define trc2(...) (void(0))
#else
#define trc2(...) (void(0))
#endif




#define each(x, v) for (auto&& x : v)
#define each2(x, y, v) for (auto&& [x, y] : v)
#define all(v) (v).begin(), (v).end()
#define rep(i, N) for (long long i = 0; i < (long long)(N); i++)
#define repr(i, N) for (long long i = (long long)(N)-1; i >= 0; i--)
#define rep1(i, N) for (long long i = 1; i <= (long long)(N); i++)
#define repr1(i, N) for (long long i = (N); (long long)(i) > 0; i--)
#define reg(i, a, b) for (long long i = (a); i < (b); i++)
#define regr(i, a, b) for (long long i = (b)-1; i >= (a); i--)
#define fi first
#define se second
#define ini(...) \
 int __VA_ARGS__; \
 in(__VA_ARGS__)
#define inl(...) \
 long long __VA_ARGS__; \
 in(__VA_ARGS__)
#define ins(...) \
 string __VA_ARGS__; \
 in(__VA_ARGS__)
#define in2(s, t) \
 for (int i = 0; i < (int)s.size(); i++) { \
 in(s[i], t[i]); \
 }
#define in3(s, t, u) \
 for (int i = 0; i < (int)s.size(); i++) { \
 in(s[i], t[i], u[i]); \
 }
#define in4(s, t, u, v) \
 for (int i = 0; i < (int)s.size(); i++) { \
 in(s[i], t[i], u[i], v[i]); \
 }
#define die(...) \
 do { \
 Nyaan::out(__VA_ARGS__); \
 return; \
 } while (0)


namespace Nyaan {
void solve();
}
int main() { Nyaan::solve(); }







template <typename T>
struct edge {
 int src, to;
 T cost;

 edge(int _to, T _cost) : src(-1), to(_to), cost(_cost) {}
 edge(int _src, int _to, T _cost) : src(_src), to(_to), cost(_cost) {}

 edge &operator=(const int &x) {
 to = x;
 return *this;
 }

 operator int() const { return to; }
};
template <typename T>
using Edges = vector<edge<T>>;
template <typename T>
using WeightedGraph = vector<Edges<T>>;
using UnweightedGraph = vector<vector<int>>;


UnweightedGraph graph(int N, int M = -1, bool is_directed = false,
 bool is_1origin = true) {
 UnweightedGraph g(N);
 if (M == -1) M = N - 1;
 for (int _ = 0; _ < M; _++) {
 int x, y;
 cin >> x >> y;
 if (is_1origin) x--, y--;
 g[x].push_back(y);
 if (!is_directed) g[y].push_back(x);
 }
 return g;
}


template <typename T>
WeightedGraph<T> wgraph(int N, int M = -1, bool is_directed = false,
 bool is_1origin = true) {
 WeightedGraph<T> g(N);
 if (M == -1) M = N - 1;
 for (int _ = 0; _ < M; _++) {
 int x, y;
 cin >> x >> y;
 T c;
 cin >> c;
 if (is_1origin) x--, y--;
 g[x].emplace_back(x, y, c);
 if (!is_directed) g[y].emplace_back(y, x, c);
 }
 return g;
}


template <typename T>
Edges<T> esgraph(int N, int M, int is_weighted = true, bool is_1origin = true) {
 Edges<T> es;
 for (int _ = 0; _ < M; _++) {
 int x, y;
 cin >> x >> y;
 T c;
 if (is_weighted)
 cin >> c;
 else
 c = 1;
 if (is_1origin) x--, y--;
 es.emplace_back(x, y, c);
 }
 return es;
}


template <typename T>
vector<vector<T>> adjgraph(int N, int M, T INF, int is_weighted = true,
 bool is_directed = false, bool is_1origin = true) {
 vector<vector<T>> d(N, vector<T>(N, INF));
 for (int _ = 0; _ < M; _++) {
 int x, y;
 cin >> x >> y;
 T c;
 if (is_weighted)
 cin >> c;
 else
 c = 1;
 if (is_1origin) x--, y--;
 d[x][y] = c;
 if (!is_directed) d[y][x] = c;
 }
 return d;
}













vector<int> Depth(const UnweightedGraph &g, int start = 0) {
 int n = g.size();
 vector<int> ds(n, -1);
 ds[start] = 0;
 queue<int> q;
 q.push(start);
 while (!q.empty()) {
 int c = q.front();
 q.pop();
 int dc = ds[c];
 for (auto &d : g[c]) {
 if (ds[d] == -1) {
 ds[d] = dc + 1;
 q.push(d);
 }
 }
 }
 return ds;
}



template <typename T>
vector<T> Depth(const WeightedGraph<T> &g, int start = 0) {
 vector<T> d(g.size(), -1);
 auto dfs = [&](auto rec, int cur, T val, int par = -1) -> void {
 d[cur] = val;
 for (auto &dst : g[cur]) {
 if (dst == par) continue;
 rec(rec, dst, val + dst.cost, cur);
 }
 };
 dfs(dfs, start, 0);
 return d;
}



pair<pair<int, int>, int> Diameter(const UnweightedGraph &g) {
 auto d = Depth(g, 0);
 int u = max_element(begin(d), end(d)) - begin(d);
 d = Depth(g, u);
 int v = max_element(begin(d), end(d)) - begin(d);
 return make_pair(make_pair(u, v), d[v]);
}



template <typename T>
pair<pair<int, int>, T> Diameter(const WeightedGraph<T> &g) {
 auto d = Depth(g, 0);
 int u = max_element(begin(d), end(d)) - begin(d);
 d = Depth(g, u);
 int v = max_element(begin(d), end(d)) - begin(d);
 return make_pair(make_pair(u, v), d[v]);
}


template <typename G>
vector<int> Path(G &g, int u, int v) {
 vector<int> ret;
 int end = 0;
 auto dfs = [&](auto rec, int cur, int par = -1) -> void {
 ret.push_back(cur);
 if (cur == v) {
 end = 1;
 return;
 }
 for (int dst : g[cur]) {
 if (dst == par) continue;
 rec(rec, dst, cur);
 if (end) return;
 }
 if (end) return;
 ret.pop_back();
 };
 dfs(dfs, u);
 return ret;
}











template <typename T>
struct has_cost {
 private:
 template <typename U>
 static auto confirm(U u) -> decltype(u.cost, std::true_type());
 static auto confirm(...) -> std::false_type;

 public:
 enum : bool { value = decltype(confirm(std::declval<T>()))::value };
};

template <typename T>
vector<vector<T>> inverse_tree(const vector<vector<T>>& g) {
 int N = (int)g.size();
 vector<vector<T>> rg(N);
 for (int i = 0; i < N; i++) {
 for (auto& e : g[i]) {
 if constexpr (is_same<T, int>::value) {
 rg[e].push_back(i);
 } else if constexpr (has_cost<T>::value) {
 rg[e].emplace_back(e.to, i, e.cost);
 } else {
 assert(0);
 }
 }
 }
 return rg;
}

template <typename T>
vector<vector<T>> rooted_tree(const vector<vector<T>>& g, int root = 0) {
 int N = (int)g.size();
 vector<vector<T>> rg(N);
 vector<char> v(N, false);
 v[root] = true;
 queue<int> que;
 que.emplace(root);
 while (!que.empty()) {
 auto p = que.front();
 que.pop();
 for (auto& e : g[p]) {
 if (v[e] == false) {
 v[e] = true;
 que.push(e);
 rg[p].push_back(e);
 }
 }
 }
 return rg;
}










template <typename G>
struct HeavyLightDecomposition {
 private:
 void dfs_sz(int cur) {
 size[cur] = 1;
 for (auto& dst : g[cur]) {
 if (dst == par[cur]) {
 if (g[cur].size() >= 2 && int(dst) == int(g[cur][0]))
 swap(g[cur][0], g[cur][1]);
 else
 continue;
 }
 depth[dst] = depth[cur] + 1;
 par[dst] = cur;
 dfs_sz(dst);
 size[cur] += size[dst];
 if (size[dst] > size[g[cur][0]]) {
 swap(dst, g[cur][0]);
 }
 }
 }

 void dfs_hld(int cur) {
 down[cur] = id++;
 for (auto dst : g[cur]) {
 if (dst == par[cur]) continue;
 nxt[dst] = (int(dst) == int(g[cur][0]) ? nxt[cur] : int(dst));
 dfs_hld(dst);
 }
 up[cur] = id;
 }

 
 vector<pair<int, int>> ascend(int u, int v) const {
 vector<pair<int, int>> res;
 while (nxt[u] != nxt[v]) {
 res.emplace_back(down[u], down[nxt[u]]);
 u = par[nxt[u]];
 }
 if (u != v) res.emplace_back(down[u], down[v] + 1);
 return res;
 }

 
 vector<pair<int, int>> descend(int u, int v) const {
 if (u == v) return {};
 if (nxt[u] == nxt[v]) return {{down[u] + 1, down[v]}};
 auto res = descend(u, par[nxt[v]]);
 res.emplace_back(down[nxt[v]], down[v]);
 return res;
 }

 public:
 G& g;
 int id;
 vector<int> size, depth, down, up, nxt, par;
 HeavyLightDecomposition(G& _g, int root = 0)
 : g(_g),
 id(0),
 size(g.size(), 0),
 depth(g.size(), 0),
 down(g.size(), -1),
 up(g.size(), -1),
 nxt(g.size(), root),
 par(g.size(), root) {
 dfs_sz(root);
 dfs_hld(root);
 }

 void build(int root) {
 dfs_sz(root);
 dfs_hld(root);
 }

 pair<int, int> idx(int i) const { return make_pair(down[i], up[i]); }

 template <typename F>
 void path_query(int u, int v, bool vertex, const F& f) {
 int l = lca(u, v);
 for (auto&& [a, b] : ascend(u, l)) {
 int s = a + 1, t = b;
 s > t ? f(t, s) : f(s, t);
 }
 if (vertex) f(down[l], down[l] + 1);
 for (auto&& [a, b] : descend(l, v)) {
 int s = a, t = b + 1;
 s > t ? f(t, s) : f(s, t);
 }
 }

 template <typename F>
 void path_noncommutative_query(int u, int v, bool vertex, const F& f) {
 int l = lca(u, v);
 for (auto&& [a, b] : ascend(u, l)) f(a + 1, b);
 if (vertex) f(down[l], down[l] + 1);
 for (auto&& [a, b] : descend(l, v)) f(a, b + 1);
 }

 template <typename F>
 void subtree_query(int u, bool vertex, const F& f) {
 f(down[u] + int(!vertex), up[u]);
 }

 int lca(int a, int b) {
 while (nxt[a] != nxt[b]) {
 if (down[a] < down[b]) swap(a, b);
 a = par[nxt[a]];
 }
 return depth[a] < depth[b] ? a : b;
 }

 int dist(int a, int b) { return depth[a] + depth[b] - depth[lca(a, b)] * 2; }
};










template <typename fps>
struct fps_fraction {
 using frac = fps_fraction;
 using mint = typename fps::value_type;

 fps p, q;
 fps_fraction(const fps& numerator = fps{0}, const fps& denominator = fps{1})
 : p(numerator), q(denominator) {}

 friend frac operator+(const frac& l, const frac& r) {
 return frac{l.p * r.q + r.p * l.q, l.q * r.q};
 }
 friend frac operator-(const frac& l, const frac& r) {
 return frac{l.p * r.q - r.p * l.q, l.q * r.q};
 }
 friend frac operator*(const frac& l, const frac& r) {
 return frac{l.p * r.p, l.q * r.q};
 }
 friend frac operator/(const frac& l, const frac& r) {
 return frac{l.p * r.q, l.q * r.p};
 }

 frac& operator+=(const mint& r) {
 (*this).p += (*this).q * r;
 return *this;
 }
 frac& operator-=(const mint& r) {
 (*this).p -= (*this).q * r;
 return *this;
 }
 frac& operator*=(const mint& r) {
 (*this).p *= r;
 return *this;
 }

 frac operator+(const mint& r) { return frac{*this} += r; }
 frac operator-(const mint& r) { return frac{*this} -= r; }
 frac operator*(const mint& r) { return frac{*this} *= r; }
 frac operator/(const mint& r) { return frac{*this} /= r; }
 frac& operator+=(const frac& r) { return *this = (*this) + r; }
 frac& operator-=(const frac& r) { return *this = (*this) - r; }
 frac& operator*=(const frac& r) { return *this = (*this) * r; }
 frac operator-() const { return frac{-p, q}; }
 frac inverse() const { return frac{q, p}; };

 void shrink() { p.shrink(), q.shrink(); }
 friend bool operator==(const frac& l, const frac& r) {
 return l.p == r.p && l.q == r.q;
 }
 friend bool operator!=(const frac& l, const frac& r) {
 return l.p != r.p || l.q != r.q;
 }
 friend ostream& operator<<(ostream& os, const frac& r) {
 os << "[ ";
 for (auto& x : r.p) os << x << ", ";
 os << "], ";
 os << "[ ";
 for (auto& x : r.q) os << x << ", ";
 os << " ]";
 return os;
 }
};





using namespace std;




template <typename T>
struct Binomial {
 vector<T> f, g, h;
 Binomial(int MAX = 0) {
 assert(T::get_mod() != 0 && "Binomial<mint>()");
 f.resize(1, T{1});
 g.resize(1, T{1});
 h.resize(1, T{1});
 if (MAX > 0) extend(MAX + 1);
 }

 void extend(int m = -1) {
 int n = f.size();
 if (m == -1) m = n * 2;
 m = min<int>(m, T::get_mod());
 if (n >= m) return;
 f.resize(m);
 g.resize(m);
 h.resize(m);
 for (int i = n; i < m; i++) f[i] = f[i - 1] * T(i);
 g[m - 1] = f[m - 1].inverse();
 h[m - 1] = g[m - 1] * f[m - 2];
 for (int i = m - 2; i >= n; i--) {
 g[i] = g[i + 1] * T(i + 1);
 h[i] = g[i] * f[i - 1];
 }
 }

 T fac(int i) {
 if (i < 0) return T(0);
 while (i >= (int)f.size()) extend();
 return f[i];
 }

 T finv(int i) {
 if (i < 0) return T(0);
 while (i >= (int)g.size()) extend();
 return g[i];
 }

 T inv(int i) {
 if (i < 0) return -inv(-i);
 while (i >= (int)h.size()) extend();
 return h[i];
 }

 T C(int n, int r) {
 if (n < 0 || n < r || r < 0) return T(0);
 return fac(n) * finv(n - r) * finv(r);
 }

 inline T operator()(int n, int r) { return C(n, r); }

 template <typename I>
 T multinomial(const vector<I>& r) {
 static_assert(is_integral<I>::value == true);
 int n = 0;
 for (auto& x : r) {
 if (x < 0) return T(0);
 n += x;
 }
 T res = fac(n);
 for (auto& x : r) res *= finv(x);
 return res;
 }

 template <typename I>
 T operator()(const vector<I>& r) {
 return multinomial(r);
 }

 T C_naive(int n, int r) {
 if (n < 0 || n < r || r < 0) return T(0);
 T ret = T(1);
 r = min(r, n - r);
 for (int i = 1; i <= r; ++i) ret *= inv(i) * (n--);
 return ret;
 }

 T P(int n, int r) {
 if (n < 0 || n < r || r < 0) return T(0);
 return fac(n) * finv(n - r);
 }

 
 T H(int n, int r) {
 if (n < 0 || r < 0) return T(0);
 return r == 0 ? 1 : C(n + r - 1, r);
 }
};


template <typename mint>
struct FormalPowerSeries : vector<mint> {
 using vector<mint>::vector;
 using FPS = FormalPowerSeries;

 FPS &operator+=(const FPS &r) {
 if (r.size() > this->size()) this->resize(r.size());
 for (int i = 0; i < (int)r.size(); i++) (*this)[i] += r[i];
 return *this;
 }

 FPS &operator+=(const mint &r) {
 if (this->empty()) this->resize(1);
 (*this)[0] += r;
 return *this;
 }

 FPS &operator-=(const FPS &r) {
 if (r.size() > this->size()) this->resize(r.size());
 for (int i = 0; i < (int)r.size(); i++) (*this)[i] -= r[i];
 return *this;
 }

 FPS &operator-=(const mint &r) {
 if (this->empty()) this->resize(1);
 (*this)[0] -= r;
 return *this;
 }

 FPS &operator*=(const mint &v) {
 for (int k = 0; k < (int)this->size(); k++) (*this)[k] *= v;
 return *this;
 }

 FPS &operator/=(const FPS &r) {
 if (this->size() < r.size()) {
 this->clear();
 return *this;
 }
 int n = this->size() - r.size() + 1;
 if ((int)r.size() <= 64) {
 FPS f(*this), g(r);
 g.shrink();
 mint coeff = g.back().inverse();
 for (auto &x : g) x *= coeff;
 int deg = (int)f.size() - (int)g.size() + 1;
 int gs = g.size();
 FPS quo(deg);
 for (int i = deg - 1; i >= 0; i--) {
 quo[i] = f[i + gs - 1];
 for (int j = 0; j < gs; j++) f[i + j] -= quo[i] * g[j];
 }
 *this = quo * coeff;
 this->resize(n, mint(0));
 return *this;
 }
 return *this = ((*this).rev().pre(n) * r.rev().inv(n)).pre(n).rev();
 }

 FPS &operator%=(const FPS &r) {
 *this -= *this / r * r;
 shrink();
 return *this;
 }

 FPS operator+(const FPS &r) const { return FPS(*this) += r; }
 FPS operator+(const mint &v) const { return FPS(*this) += v; }
 FPS operator-(const FPS &r) const { return FPS(*this) -= r; }
 FPS operator-(const mint &v) const { return FPS(*this) -= v; }
 FPS operator*(const FPS &r) const { return FPS(*this) *= r; }
 FPS operator*(const mint &v) const { return FPS(*this) *= v; }
 FPS operator/(const FPS &r) const { return FPS(*this) /= r; }
 FPS operator%(const FPS &r) const { return FPS(*this) %= r; }
 FPS operator-() const {
 FPS ret(this->size());
 for (int i = 0; i < (int)this->size(); i++) ret[i] = -(*this)[i];
 return ret;
 }

 void shrink() {
 while (this->size() && this->back() == mint(0)) this->pop_back();
 }

 FPS rev() const {
 FPS ret(*this);
 reverse(begin(ret), end(ret));
 return ret;
 }

 FPS dot(FPS r) const {
 FPS ret(min(this->size(), r.size()));
 for (int i = 0; i < (int)ret.size(); i++) ret[i] = (*this)[i] * r[i];
 return ret;
 }

 
 FPS pre(int sz) const {
 FPS ret(begin(*this), begin(*this) + min((int)this->size(), sz));
 if ((int)ret.size() < sz) ret.resize(sz);
 return ret;
 }

 FPS operator>>(int sz) const {
 if ((int)this->size() <= sz) return {};
 FPS ret(*this);
 ret.erase(ret.begin(), ret.begin() + sz);
 return ret;
 }

 FPS operator<<(int sz) const {
 FPS ret(*this);
 ret.insert(ret.begin(), sz, mint(0));
 return ret;
 }

 FPS diff() const {
 const int n = (int)this->size();
 FPS ret(max(0, n - 1));
 mint one(1), coeff(1);
 for (int i = 1; i < n; i++) {
 ret[i - 1] = (*this)[i] * coeff;
 coeff += one;
 }
 return ret;
 }

 FPS integral() const {
 const int n = (int)this->size();
 FPS ret(n + 1);
 ret[0] = mint(0);
 if (n > 0) ret[1] = mint(1);
 auto mod = mint::get_mod();
 for (int i = 2; i <= n; i++) ret[i] = (-ret[mod % i]) * (mod / i);
 for (int i = 0; i < n; i++) ret[i + 1] *= (*this)[i];
 return ret;
 }

 mint eval(mint x) const {
 mint r = 0, w = 1;
 for (auto &v : *this) r += w * v, w *= x;
 return r;
 }

 FPS log(int deg = -1) const {
 assert(!(*this).empty() && (*this)[0] == mint(1));
 if (deg == -1) deg = (int)this->size();
 return (this->diff() * this->inv(deg)).pre(deg - 1).integral();
 }

 FPS pow(int64_t k, int deg = -1) const {
 const int n = (int)this->size();
 if (deg == -1) deg = n;
 if (k == 0) {
 FPS ret(deg);
 if (deg) ret[0] = 1;
 return ret;
 }
 for (int i = 0; i < n; i++) {
 if ((*this)[i] != mint(0)) {
 mint rev = mint(1) / (*this)[i];
 FPS ret = (((*this * rev) >> i).log(deg) * k).exp(deg);
 ret *= (*this)[i].pow(k);
 ret = (ret << (i * k)).pre(deg);
 if ((int)ret.size() < deg) ret.resize(deg, mint(0));
 return ret;
 }
 if (__int128_t(i + 1) * k >= deg) return FPS(deg, mint(0));
 }
 return FPS(deg, mint(0));
 }

 static void *ntt_ptr;
 static void set_fft();
 FPS &operator*=(const FPS &r);
 void ntt();
 void intt();
 void ntt_doubling();
 static int ntt_pr();
 FPS inv(int deg = -1) const;
 FPS exp(int deg = -1) const;
};
template <typename mint>
void *FormalPowerSeries<mint>::ntt_ptr = nullptr;







template <typename fps>
fps Pi(vector<fps> v) {
 if ((int)v.size() == 0) return fps{1};
 while ((int)v.size() >= 2) {
 vector<fps> nx;
 for (int i = 0; i + 1 < (int)v.size(); i += 2)
 nx.push_back(v[i] * v[i + 1]);
 if (v.size() % 2) nx.push_back(v.back());
 v = nx;
 }
 return v.back();
}






template <typename mint>
mint LinearRecurrence(long long k, FormalPowerSeries<mint> Q,
 FormalPowerSeries<mint> P) {
 Q.shrink();
 mint ret = 0;
 if (P.size() >= Q.size()) {
 auto R = P / Q;
 P -= R * Q;
 P.shrink();
 if (k < (int)R.size()) ret += R[k];
 }
 if ((int)P.size() == 0) return ret;

 FormalPowerSeries<mint>::set_fft();
 if (FormalPowerSeries<mint>::ntt_ptr == nullptr) {
 P.resize((int)Q.size() - 1);
 while (k) {
 auto Q2 = Q;
 for (int i = 1; i < (int)Q2.size(); i += 2) Q2[i] = -Q2[i];
 auto S = P * Q2;
 auto T = Q * Q2;
 if (k & 1) {
 for (int i = 1; i < (int)S.size(); i += 2) P[i >> 1] = S[i];
 for (int i = 0; i < (int)T.size(); i += 2) Q[i >> 1] = T[i];
 } else {
 for (int i = 0; i < (int)S.size(); i += 2) P[i >> 1] = S[i];
 for (int i = 0; i < (int)T.size(); i += 2) Q[i >> 1] = T[i];
 }
 k >>= 1;
 }
 return ret + P[0];
 } else {
 int N = 1;
 while (N < (int)Q.size()) N <<= 1;

 P.resize(2 * N);
 Q.resize(2 * N);
 P.ntt();
 Q.ntt();
 vector<mint> S(2 * N), T(2 * N);

 vector<int> btr(N);
 for (int i = 0, logn = __builtin_ctz(N); i < (1 << logn); i++) {
 btr[i] = (btr[i >> 1] >> 1) + ((i & 1) << (logn - 1));
 }
 mint dw = mint(FormalPowerSeries<mint>::ntt_pr())
 .inverse()
 .pow((mint::get_mod() - 1) / (2 * N));

 while (k) {
 mint inv2 = mint(2).inverse();

 
 T.resize(N);
 for (int i = 0; i < N; i++) T[i] = Q[(i << 1) | 0] * Q[(i << 1) | 1];

 S.resize(N);
 if (k & 1) {
 
 for (auto &i : btr) {
 S[i] = (P[(i << 1) | 0] * Q[(i << 1) | 1] -
 P[(i << 1) | 1] * Q[(i << 1) | 0]) *
 inv2;
 inv2 *= dw;
 }
 } else {
 
 for (int i = 0; i < N; i++) {
 S[i] = (P[(i << 1) | 0] * Q[(i << 1) | 1] +
 P[(i << 1) | 1] * Q[(i << 1) | 0]) *
 inv2;
 }
 }
 swap(P, S);
 swap(Q, T);
 k >>= 1;
 if (k < N) break;
 P.ntt_doubling();
 Q.ntt_doubling();
 }
 P.intt();
 Q.intt();
 return ret + (P * (Q.inv()))[k];
 }
}

template <typename mint>
mint kitamasa(long long N, FormalPowerSeries<mint> Q,
 FormalPowerSeries<mint> a) {
 assert(!Q.empty() && Q[0] != 0);
 if (N < (int)a.size()) return a[N];
 assert((int)a.size() >= int(Q.size()) - 1);
 auto P = a.pre((int)Q.size() - 1) * Q;
 P.resize(Q.size() - 1);
 return LinearRecurrence<mint>(N, Q, P);
}














using namespace std;





template <typename T>
std::pair<int, T> GaussElimination(vector<vector<T>> &a, int pivot_end = -1,
 bool diagonalize = false) {
 int H = a.size(), W = a[0].size(), rank = 0;
 if (pivot_end == -1) pivot_end = W;
 T det = 1;
 for (int j = 0; j < pivot_end; j++) {
 int idx = -1;
 for (int i = rank; i < H; i++) {
 if (a[i][j] != T(0)) {
 idx = i;
 break;
 }
 }
 if (idx == -1) {
 det = 0;
 continue;
 }
 if (rank != idx) det = -det, swap(a[rank], a[idx]);
 det *= a[rank][j];
 if (diagonalize && a[rank][j] != T(1)) {
 T coeff = T(1) / a[rank][j];
 for (int k = j; k < W; k++) a[rank][k] *= coeff;
 }
 int is = diagonalize ? 0 : rank + 1;
 for (int i = is; i < H; i++) {
 if (i == rank) continue;
 if (a[i][j] != T(0)) {
 T coeff = a[i][j] / a[rank][j];
 for (int k = j; k < W; k++) a[i][k] -= a[rank][k] * coeff;
 }
 }
 rank++;
 }
 return make_pair(rank, det);
}


template <typename mint>
vector<vector<mint>> inverse_matrix(const vector<vector<mint>>& a) {
 int N = a.size();
 assert(N > 0);
 assert(N == (int)a[0].size());

 vector<vector<mint>> m(N, vector<mint>(2 * N));
 for (int i = 0; i < N; i++) {
 copy(begin(a[i]), end(a[i]), begin(m[i]));
 m[i][N + i] = 1;
 }

 auto [rank, det] = GaussElimination(m, N, true);
 if (rank != N) return {};

 vector<vector<mint>> b(N);
 for (int i = 0; i < N; i++) {
 copy(begin(m[i]) + N, end(m[i]), back_inserter(b[i]));
 }
 return b;
}


template <class T>
struct Matrix {
 vector<vector<T> > A;

 Matrix() = default;
 Matrix(int n, int m) : A(n, vector<T>(m, T())) {}
 Matrix(int n) : A(n, vector<T>(n, T())){};

 int H() const { return A.size(); }

 int W() const { return A[0].size(); }

 int size() const { return A.size(); }

 inline const vector<T> &operator[](int k) const { return A[k]; }

 inline vector<T> &operator[](int k) { return A[k]; }

 static Matrix I(int n) {
 Matrix mat(n);
 for (int i = 0; i < n; i++) mat[i][i] = 1;
 return (mat);
 }

 Matrix &operator+=(const Matrix &B) {
 int n = H(), m = W();
 assert(n == B.H() && m == B.W());
 for (int i = 0; i < n; i++)
 for (int j = 0; j < m; j++) (*this)[i][j] += B[i][j];
 return (*this);
 }

 Matrix &operator-=(const Matrix &B) {
 int n = H(), m = W();
 assert(n == B.H() && m == B.W());
 for (int i = 0; i < n; i++)
 for (int j = 0; j < m; j++) (*this)[i][j] -= B[i][j];
 return (*this);
 }

 Matrix &operator*=(const Matrix &B) {
 int n = H(), m = B.W(), p = W();
 assert(p == B.H());
 vector<vector<T> > C(n, vector<T>(m, T{}));
 for (int i = 0; i < n; i++)
 for (int k = 0; k < p; k++)
 for (int j = 0; j < m; j++) C[i][j] += (*this)[i][k] * B[k][j];
 A.swap(C);
 return (*this);
 }

 Matrix &operator^=(long long k) {
 Matrix B = Matrix::I(H());
 while (k > 0) {
 if (k & 1) B *= *this;
 *this *= *this;
 k >>= 1LL;
 }
 A.swap(B.A);
 return (*this);
 }

 Matrix operator+(const Matrix &B) const { return (Matrix(*this) += B); }

 Matrix operator-(const Matrix &B) const { return (Matrix(*this) -= B); }

 Matrix operator*(const Matrix &B) const { return (Matrix(*this) *= B); }

 Matrix operator^(const long long k) const { return (Matrix(*this) ^= k); }

 bool operator==(const Matrix &B) const {
 assert(H() == B.H() && W() == B.W());
 for (int i = 0; i < H(); i++)
 for (int j = 0; j < W(); j++)
 if (A[i][j] != B[i][j]) return false;
 return true;
 }

 bool operator!=(const Matrix &B) const {
 assert(H() == B.H() && W() == B.W());
 for (int i = 0; i < H(); i++)
 for (int j = 0; j < W(); j++)
 if (A[i][j] != B[i][j]) return true;
 return false;
 }

 Matrix inverse() const {
 assert(H() == W());
 Matrix B(H());
 B.A = inverse_matrix(A);
 return B;
 }

 friend ostream &operator<<(ostream &os, const Matrix &p) {
 int n = p.H(), m = p.W();
 for (int i = 0; i < n; i++) {
 os << (i ? "   " : "") << "[";
 for (int j = 0; j < m; j++) {
 os << p[i][j] << (j + 1 == m ? "]\n" : ",");
 }
 }
 return (os);
 }

 T determinant() const {
 Matrix B(*this);
 assert(H() == W());
 T ret = 1;
 for (int i = 0; i < H(); i++) {
 int idx = -1;
 for (int j = i; j < W(); j++) {
 if (B[j][i] != 0) {
 idx = j;
 break;
 }
 }
 if (idx == -1) return 0;
 if (i != idx) {
 ret *= T(-1);
 swap(B[i], B[idx]);
 }
 ret *= B[i][i];
 T inv = T(1) / B[i][i];
 for (int j = 0; j < W(); j++) {
 B[i][j] *= inv;
 }
 for (int j = i + 1; j < H(); j++) {
 T a = B[j][i];
 if (a == 0) continue;
 for (int k = i; k < W(); k++) {
 B[j][k] -= B[i][k] * a;
 }
 }
 }
 return ret;
 }
};















__attribute__((target("sse4.2"))) inline __m128i my128_mullo_epu32(
 const __m128i &a, const __m128i &b) {
 return _mm_mullo_epi32(a, b);
}

__attribute__((target("sse4.2"))) inline __m128i my128_mulhi_epu32(
 const __m128i &a, const __m128i &b) {
 __m128i a13 = _mm_shuffle_epi32(a, 0xF5);
 __m128i b13 = _mm_shuffle_epi32(b, 0xF5);
 __m128i prod02 = _mm_mul_epu32(a, b);
 __m128i prod13 = _mm_mul_epu32(a13, b13);
 __m128i prod = _mm_unpackhi_epi64(_mm_unpacklo_epi32(prod02, prod13),
 _mm_unpackhi_epi32(prod02, prod13));
 return prod;
}

__attribute__((target("sse4.2"))) inline __m128i montgomery_mul_128(
 const __m128i &a, const __m128i &b, const __m128i &r, const __m128i &m1) {
 return _mm_sub_epi32(
 _mm_add_epi32(my128_mulhi_epu32(a, b), m1),
 my128_mulhi_epu32(my128_mullo_epu32(my128_mullo_epu32(a, b), r), m1));
}

__attribute__((target("sse4.2"))) inline __m128i montgomery_add_128(
 const __m128i &a, const __m128i &b, const __m128i &m2, const __m128i &m0) {
 __m128i ret = _mm_sub_epi32(_mm_add_epi32(a, b), m2);
 return _mm_add_epi32(_mm_and_si128(_mm_cmpgt_epi32(m0, ret), m2), ret);
}

__attribute__((target("sse4.2"))) inline __m128i montgomery_sub_128(
 const __m128i &a, const __m128i &b, const __m128i &m2, const __m128i &m0) {
 __m128i ret = _mm_sub_epi32(a, b);
 return _mm_add_epi32(_mm_and_si128(_mm_cmpgt_epi32(m0, ret), m2), ret);
}

__attribute__((target("avx2"))) inline __m256i my256_mullo_epu32(
 const __m256i &a, const __m256i &b) {
 return _mm256_mullo_epi32(a, b);
}

__attribute__((target("avx2"))) inline __m256i my256_mulhi_epu32(
 const __m256i &a, const __m256i &b) {
 __m256i a13 = _mm256_shuffle_epi32(a, 0xF5);
 __m256i b13 = _mm256_shuffle_epi32(b, 0xF5);
 __m256i prod02 = _mm256_mul_epu32(a, b);
 __m256i prod13 = _mm256_mul_epu32(a13, b13);
 __m256i prod = _mm256_unpackhi_epi64(_mm256_unpacklo_epi32(prod02, prod13),
 _mm256_unpackhi_epi32(prod02, prod13));
 return prod;
}

__attribute__((target("avx2"))) inline __m256i montgomery_mul_256(
 const __m256i &a, const __m256i &b, const __m256i &r, const __m256i &m1) {
 return _mm256_sub_epi32(
 _mm256_add_epi32(my256_mulhi_epu32(a, b), m1),
 my256_mulhi_epu32(my256_mullo_epu32(my256_mullo_epu32(a, b), r), m1));
}

__attribute__((target("avx2"))) inline __m256i montgomery_add_256(
 const __m256i &a, const __m256i &b, const __m256i &m2, const __m256i &m0) {
 __m256i ret = _mm256_sub_epi32(_mm256_add_epi32(a, b), m2);
 return _mm256_add_epi32(_mm256_and_si256(_mm256_cmpgt_epi32(m0, ret), m2),
 ret);
}

__attribute__((target("avx2"))) inline __m256i montgomery_sub_256(
 const __m256i &a, const __m256i &b, const __m256i &m2, const __m256i &m0) {
 __m256i ret = _mm256_sub_epi32(a, b);
 return _mm256_add_epi32(_mm256_and_si256(_mm256_cmpgt_epi32(m0, ret), m2),
 ret);
}

namespace ntt_inner {
using u64 = uint64_t;
constexpr uint32_t get_pr(uint32_t mod) {
 if (mod == 2) return 1;
 u64 ds[32] = {};
 int idx = 0;
 u64 m = mod - 1;
 for (u64 i = 2; i * i <= m; ++i) {
 if (m % i == 0) {
 ds[idx++] = i;
 while (m % i == 0) m /= i;
 }
 }
 if (m != 1) ds[idx++] = m;

 uint32_t pr = 2;
 while (1) {
 int flg = 1;
 for (int i = 0; i < idx; ++i) {
 u64 a = pr, b = (mod - 1) / ds[i], r = 1;
 while (b) {
 if (b & 1) r = r * a % mod;
 a = a * a % mod;
 b >>= 1;
 }
 if (r == 1) {
 flg = 0;
 break;
 }
 }
 if (flg == 1) break;
 ++pr;
 }
 return pr;
}

constexpr int SZ_FFT_BUF = 1 << 23;
uint32_t _buf1[SZ_FFT_BUF] __attribute__((aligned(64)));
uint32_t _buf2[SZ_FFT_BUF] __attribute__((aligned(64)));
} 

template <typename mint>
struct NTT {
 static constexpr uint32_t mod = mint::get_mod();
 static constexpr uint32_t pr = ntt_inner::get_pr(mint::get_mod());
 static constexpr int level = __builtin_ctzll(mod - 1);
 mint dw[level], dy[level];
 mint *buf1, *buf2;

 constexpr NTT() {
 setwy(level);
 union raw_cast {
 mint dat;
 uint32_t _;
 };
 buf1 = &(((raw_cast *)(ntt_inner::_buf1))->dat);
 buf2 = &(((raw_cast *)(ntt_inner::_buf2))->dat);
 }

 constexpr void setwy(int k) {
 mint w[level], y[level];
 w[k - 1] = mint(pr).pow((mod - 1) / (1 << k));
 y[k - 1] = w[k - 1].inverse();
 for (int i = k - 2; i > 0; --i)
 w[i] = w[i + 1] * w[i + 1], y[i] = y[i + 1] * y[i + 1];
 dw[0] = dy[0] = w[1] * w[1];
 dw[1] = w[1], dy[1] = y[1], dw[2] = w[2], dy[2] = y[2];
 for (int i = 3; i < k; ++i) {
 dw[i] = dw[i - 1] * y[i - 2] * w[i];
 dy[i] = dy[i - 1] * w[i - 2] * y[i];
 }
 }

 __attribute__((target("avx2"))) void ntt(mint *a, int n) {
 int k = n ? __builtin_ctz(n) : 0;
 if (k == 0) return;
 if (k == 1) {
 mint a1 = a[1];
 a[1] = a[0] - a[1];
 a[0] = a[0] + a1;
 return;
 }
 if (k & 1) {
 int v = 1 << (k - 1);
 if (v < 8) {
 for (int j = 0; j < v; ++j) {
 mint ajv = a[j + v];
 a[j + v] = a[j] - ajv;
 a[j] += ajv;
 }
 } else {
 const __m256i m0 = _mm256_set1_epi32(0);
 const __m256i m2 = _mm256_set1_epi32(mod + mod);
 int j0 = 0;
 int j1 = v;
 for (; j0 < v; j0 += 8, j1 += 8) {
 __m256i T0 = _mm256_loadu_si256((__m256i *)(a + j0));
 __m256i T1 = _mm256_loadu_si256((__m256i *)(a + j1));
 __m256i naj = montgomery_add_256(T0, T1, m2, m0);
 __m256i najv = montgomery_sub_256(T0, T1, m2, m0);
 _mm256_storeu_si256((__m256i *)(a + j0), naj);
 _mm256_storeu_si256((__m256i *)(a + j1), najv);
 }
 }
 }
 int u = 1 << (2 + (k & 1));
 int v = 1 << (k - 2 - (k & 1));
 mint one = mint(1);
 mint imag = dw[1];
 while (v) {
 if (v == 1) {
 mint ww = one, xx = one, wx = one;
 for (int jh = 0; jh < u;) {
 ww = xx * xx, wx = ww * xx;
 mint t0 = a[jh + 0], t1 = a[jh + 1] * xx;
 mint t2 = a[jh + 2] * ww, t3 = a[jh + 3] * wx;
 mint t0p2 = t0 + t2, t1p3 = t1 + t3;
 mint t0m2 = t0 - t2, t1m3 = (t1 - t3) * imag;
 a[jh + 0] = t0p2 + t1p3, a[jh + 1] = t0p2 - t1p3;
 a[jh + 2] = t0m2 + t1m3, a[jh + 3] = t0m2 - t1m3;
 xx *= dw[__builtin_ctz((jh += 4))];
 }
 } else if (v == 4) {
 const __m128i m0 = _mm_set1_epi32(0);
 const __m128i m1 = _mm_set1_epi32(mod);
 const __m128i m2 = _mm_set1_epi32(mod + mod);
 const __m128i r = _mm_set1_epi32(mint::r);
 const __m128i Imag = _mm_set1_epi32(imag.a);
 mint ww = one, xx = one, wx = one;
 for (int jh = 0; jh < u;) {
 if (jh == 0) {
 int j0 = 0;
 int j1 = v;
 int j2 = j1 + v;
 int j3 = j2 + v;
 int je = v;
 for (; j0 < je; j0 += 4, j1 += 4, j2 += 4, j3 += 4) {
 const __m128i T0 = _mm_loadu_si128((__m128i *)(a + j0));
 const __m128i T1 = _mm_loadu_si128((__m128i *)(a + j1));
 const __m128i T2 = _mm_loadu_si128((__m128i *)(a + j2));
 const __m128i T3 = _mm_loadu_si128((__m128i *)(a + j3));
 const __m128i T0P2 = montgomery_add_128(T0, T2, m2, m0);
 const __m128i T1P3 = montgomery_add_128(T1, T3, m2, m0);
 const __m128i T0M2 = montgomery_sub_128(T0, T2, m2, m0);
 const __m128i T1M3 = montgomery_mul_128(
 montgomery_sub_128(T1, T3, m2, m0), Imag, r, m1);
 _mm_storeu_si128((__m128i *)(a + j0),
 montgomery_add_128(T0P2, T1P3, m2, m0));
 _mm_storeu_si128((__m128i *)(a + j1),
 montgomery_sub_128(T0P2, T1P3, m2, m0));
 _mm_storeu_si128((__m128i *)(a + j2),
 montgomery_add_128(T0M2, T1M3, m2, m0));
 _mm_storeu_si128((__m128i *)(a + j3),
 montgomery_sub_128(T0M2, T1M3, m2, m0));
 }
 } else {
 ww = xx * xx, wx = ww * xx;
 const __m128i WW = _mm_set1_epi32(ww.a);
 const __m128i WX = _mm_set1_epi32(wx.a);
 const __m128i XX = _mm_set1_epi32(xx.a);
 int j0 = jh * v;
 int j1 = j0 + v;
 int j2 = j1 + v;
 int j3 = j2 + v;
 int je = j1;
 for (; j0 < je; j0 += 4, j1 += 4, j2 += 4, j3 += 4) {
 const __m128i T0 = _mm_loadu_si128((__m128i *)(a + j0));
 const __m128i T1 = _mm_loadu_si128((__m128i *)(a + j1));
 const __m128i T2 = _mm_loadu_si128((__m128i *)(a + j2));
 const __m128i T3 = _mm_loadu_si128((__m128i *)(a + j3));
 const __m128i MT1 = montgomery_mul_128(T1, XX, r, m1);
 const __m128i MT2 = montgomery_mul_128(T2, WW, r, m1);
 const __m128i MT3 = montgomery_mul_128(T3, WX, r, m1);
 const __m128i T0P2 = montgomery_add_128(T0, MT2, m2, m0);
 const __m128i T1P3 = montgomery_add_128(MT1, MT3, m2, m0);
 const __m128i T0M2 = montgomery_sub_128(T0, MT2, m2, m0);
 const __m128i T1M3 = montgomery_mul_128(
 montgomery_sub_128(MT1, MT3, m2, m0), Imag, r, m1);
 _mm_storeu_si128((__m128i *)(a + j0),
 montgomery_add_128(T0P2, T1P3, m2, m0));
 _mm_storeu_si128((__m128i *)(a + j1),
 montgomery_sub_128(T0P2, T1P3, m2, m0));
 _mm_storeu_si128((__m128i *)(a + j2),
 montgomery_add_128(T0M2, T1M3, m2, m0));
 _mm_storeu_si128((__m128i *)(a + j3),
 montgomery_sub_128(T0M2, T1M3, m2, m0));
 }
 }
 xx *= dw[__builtin_ctz((jh += 4))];
 }
 } else {
 const __m256i m0 = _mm256_set1_epi32(0);
 const __m256i m1 = _mm256_set1_epi32(mod);
 const __m256i m2 = _mm256_set1_epi32(mod + mod);
 const __m256i r = _mm256_set1_epi32(mint::r);
 const __m256i Imag = _mm256_set1_epi32(imag.a);
 mint ww = one, xx = one, wx = one;
 for (int jh = 0; jh < u;) {
 if (jh == 0) {
 int j0 = 0;
 int j1 = v;
 int j2 = j1 + v;
 int j3 = j2 + v;
 int je = v;
 for (; j0 < je; j0 += 8, j1 += 8, j2 += 8, j3 += 8) {
 const __m256i T0 = _mm256_loadu_si256((__m256i *)(a + j0));
 const __m256i T1 = _mm256_loadu_si256((__m256i *)(a + j1));
 const __m256i T2 = _mm256_loadu_si256((__m256i *)(a + j2));
 const __m256i T3 = _mm256_loadu_si256((__m256i *)(a + j3));
 const __m256i T0P2 = montgomery_add_256(T0, T2, m2, m0);
 const __m256i T1P3 = montgomery_add_256(T1, T3, m2, m0);
 const __m256i T0M2 = montgomery_sub_256(T0, T2, m2, m0);
 const __m256i T1M3 = montgomery_mul_256(
 montgomery_sub_256(T1, T3, m2, m0), Imag, r, m1);
 _mm256_storeu_si256((__m256i *)(a + j0),
 montgomery_add_256(T0P2, T1P3, m2, m0));
 _mm256_storeu_si256((__m256i *)(a + j1),
 montgomery_sub_256(T0P2, T1P3, m2, m0));
 _mm256_storeu_si256((__m256i *)(a + j2),
 montgomery_add_256(T0M2, T1M3, m2, m0));
 _mm256_storeu_si256((__m256i *)(a + j3),
 montgomery_sub_256(T0M2, T1M3, m2, m0));
 }
 } else {
 ww = xx * xx, wx = ww * xx;
 const __m256i WW = _mm256_set1_epi32(ww.a);
 const __m256i WX = _mm256_set1_epi32(wx.a);
 const __m256i XX = _mm256_set1_epi32(xx.a);
 int j0 = jh * v;
 int j1 = j0 + v;
 int j2 = j1 + v;
 int j3 = j2 + v;
 int je = j1;
 for (; j0 < je; j0 += 8, j1 += 8, j2 += 8, j3 += 8) {
 const __m256i T0 = _mm256_loadu_si256((__m256i *)(a + j0));
 const __m256i T1 = _mm256_loadu_si256((__m256i *)(a + j1));
 const __m256i T2 = _mm256_loadu_si256((__m256i *)(a + j2));
 const __m256i T3 = _mm256_loadu_si256((__m256i *)(a + j3));
 const __m256i MT1 = montgomery_mul_256(T1, XX, r, m1);
 const __m256i MT2 = montgomery_mul_256(T2, WW, r, m1);
 const __m256i MT3 = montgomery_mul_256(T3, WX, r, m1);
 const __m256i T0P2 = montgomery_add_256(T0, MT2, m2, m0);
 const __m256i T1P3 = montgomery_add_256(MT1, MT3, m2, m0);
 const __m256i T0M2 = montgomery_sub_256(T0, MT2, m2, m0);
 const __m256i T1M3 = montgomery_mul_256(
 montgomery_sub_256(MT1, MT3, m2, m0), Imag, r, m1);
 _mm256_storeu_si256((__m256i *)(a + j0),
 montgomery_add_256(T0P2, T1P3, m2, m0));
 _mm256_storeu_si256((__m256i *)(a + j1),
 montgomery_sub_256(T0P2, T1P3, m2, m0));
 _mm256_storeu_si256((__m256i *)(a + j2),
 montgomery_add_256(T0M2, T1M3, m2, m0));
 _mm256_storeu_si256((__m256i *)(a + j3),
 montgomery_sub_256(T0M2, T1M3, m2, m0));
 }
 }
 xx *= dw[__builtin_ctz((jh += 4))];
 }
 }
 u <<= 2;
 v >>= 2;
 }
 }

 __attribute__((target("avx2"))) void intt(mint *a, int n,
 int normalize = true) {
 int k = n ? __builtin_ctz(n) : 0;
 if (k == 0) return;
 if (k == 1) {
 mint a1 = a[1];
 a[1] = a[0] - a[1];
 a[0] = a[0] + a1;
 if (normalize) {
 a[0] *= mint(2).inverse();
 a[1] *= mint(2).inverse();
 }
 return;
 }
 int u = 1 << (k - 2);
 int v = 1;
 mint one = mint(1);
 mint imag = dy[1];
 while (u) {
 if (v == 1) {
 mint ww = one, xx = one, yy = one;
 u <<= 2;
 for (int jh = 0; jh < u;) {
 ww = xx * xx, yy = xx * imag;
 mint t0 = a[jh + 0], t1 = a[jh + 1];
 mint t2 = a[jh + 2], t3 = a[jh + 3];
 mint t0p1 = t0 + t1, t2p3 = t2 + t3;
 mint t0m1 = (t0 - t1) * xx, t2m3 = (t2 - t3) * yy;
 a[jh + 0] = t0p1 + t2p3, a[jh + 2] = (t0p1 - t2p3) * ww;
 a[jh + 1] = t0m1 + t2m3, a[jh + 3] = (t0m1 - t2m3) * ww;
 xx *= dy[__builtin_ctz(jh += 4)];
 }
 } else if (v == 4) {
 const __m128i m0 = _mm_set1_epi32(0);
 const __m128i m1 = _mm_set1_epi32(mod);
 const __m128i m2 = _mm_set1_epi32(mod + mod);
 const __m128i r = _mm_set1_epi32(mint::r);
 const __m128i Imag = _mm_set1_epi32(imag.a);
 mint ww = one, xx = one, yy = one;
 u <<= 2;
 for (int jh = 0; jh < u;) {
 if (jh == 0) {
 int j0 = 0;
 int j1 = v;
 int j2 = v + v;
 int j3 = j2 + v;
 for (; j0 < v; j0 += 4, j1 += 4, j2 += 4, j3 += 4) {
 const __m128i T0 = _mm_loadu_si128((__m128i *)(a + j0));
 const __m128i T1 = _mm_loadu_si128((__m128i *)(a + j1));
 const __m128i T2 = _mm_loadu_si128((__m128i *)(a + j2));
 const __m128i T3 = _mm_loadu_si128((__m128i *)(a + j3));
 const __m128i T0P1 = montgomery_add_128(T0, T1, m2, m0);
 const __m128i T2P3 = montgomery_add_128(T2, T3, m2, m0);
 const __m128i T0M1 = montgomery_sub_128(T0, T1, m2, m0);
 const __m128i T2M3 = montgomery_mul_128(
 montgomery_sub_128(T2, T3, m2, m0), Imag, r, m1);
 _mm_storeu_si128((__m128i *)(a + j0),
 montgomery_add_128(T0P1, T2P3, m2, m0));
 _mm_storeu_si128((__m128i *)(a + j2),
 montgomery_sub_128(T0P1, T2P3, m2, m0));
 _mm_storeu_si128((__m128i *)(a + j1),
 montgomery_add_128(T0M1, T2M3, m2, m0));
 _mm_storeu_si128((__m128i *)(a + j3),
 montgomery_sub_128(T0M1, T2M3, m2, m0));
 }
 } else {
 ww = xx * xx, yy = xx * imag;
 const __m128i WW = _mm_set1_epi32(ww.a);
 const __m128i XX = _mm_set1_epi32(xx.a);
 const __m128i YY = _mm_set1_epi32(yy.a);
 int j0 = jh * v;
 int j1 = j0 + v;
 int j2 = j1 + v;
 int j3 = j2 + v;
 int je = j1;
 for (; j0 < je; j0 += 4, j1 += 4, j2 += 4, j3 += 4) {
 const __m128i T0 = _mm_loadu_si128((__m128i *)(a + j0));
 const __m128i T1 = _mm_loadu_si128((__m128i *)(a + j1));
 const __m128i T2 = _mm_loadu_si128((__m128i *)(a + j2));
 const __m128i T3 = _mm_loadu_si128((__m128i *)(a + j3));
 const __m128i T0P1 = montgomery_add_128(T0, T1, m2, m0);
 const __m128i T2P3 = montgomery_add_128(T2, T3, m2, m0);
 const __m128i T0M1 = montgomery_mul_128(
 montgomery_sub_128(T0, T1, m2, m0), XX, r, m1);
 __m128i T2M3 = montgomery_mul_128(
 montgomery_sub_128(T2, T3, m2, m0), YY, r, m1);
 _mm_storeu_si128((__m128i *)(a + j0),
 montgomery_add_128(T0P1, T2P3, m2, m0));
 _mm_storeu_si128(
 (__m128i *)(a + j2),
 montgomery_mul_128(montgomery_sub_128(T0P1, T2P3, m2, m0), WW,
 r, m1));
 _mm_storeu_si128((__m128i *)(a + j1),
 montgomery_add_128(T0M1, T2M3, m2, m0));
 _mm_storeu_si128(
 (__m128i *)(a + j3),
 montgomery_mul_128(montgomery_sub_128(T0M1, T2M3, m2, m0), WW,
 r, m1));
 }
 }
 xx *= dy[__builtin_ctz(jh += 4)];
 }
 } else {
 const __m256i m0 = _mm256_set1_epi32(0);
 const __m256i m1 = _mm256_set1_epi32(mod);
 const __m256i m2 = _mm256_set1_epi32(mod + mod);
 const __m256i r = _mm256_set1_epi32(mint::r);
 const __m256i Imag = _mm256_set1_epi32(imag.a);
 mint ww = one, xx = one, yy = one;
 u <<= 2;
 for (int jh = 0; jh < u;) {
 if (jh == 0) {
 int j0 = 0;
 int j1 = v;
 int j2 = v + v;
 int j3 = j2 + v;
 for (; j0 < v; j0 += 8, j1 += 8, j2 += 8, j3 += 8) {
 const __m256i T0 = _mm256_loadu_si256((__m256i *)(a + j0));
 const __m256i T1 = _mm256_loadu_si256((__m256i *)(a + j1));
 const __m256i T2 = _mm256_loadu_si256((__m256i *)(a + j2));
 const __m256i T3 = _mm256_loadu_si256((__m256i *)(a + j3));
 const __m256i T0P1 = montgomery_add_256(T0, T1, m2, m0);
 const __m256i T2P3 = montgomery_add_256(T2, T3, m2, m0);
 const __m256i T0M1 = montgomery_sub_256(T0, T1, m2, m0);
 const __m256i T2M3 = montgomery_mul_256(
 montgomery_sub_256(T2, T3, m2, m0), Imag, r, m1);
 _mm256_storeu_si256((__m256i *)(a + j0),
 montgomery_add_256(T0P1, T2P3, m2, m0));
 _mm256_storeu_si256((__m256i *)(a + j2),
 montgomery_sub_256(T0P1, T2P3, m2, m0));
 _mm256_storeu_si256((__m256i *)(a + j1),
 montgomery_add_256(T0M1, T2M3, m2, m0));
 _mm256_storeu_si256((__m256i *)(a + j3),
 montgomery_sub_256(T0M1, T2M3, m2, m0));
 }
 } else {
 ww = xx * xx, yy = xx * imag;
 const __m256i WW = _mm256_set1_epi32(ww.a);
 const __m256i XX = _mm256_set1_epi32(xx.a);
 const __m256i YY = _mm256_set1_epi32(yy.a);
 int j0 = jh * v;
 int j1 = j0 + v;
 int j2 = j1 + v;
 int j3 = j2 + v;
 int je = j1;
 for (; j0 < je; j0 += 8, j1 += 8, j2 += 8, j3 += 8) {
 const __m256i T0 = _mm256_loadu_si256((__m256i *)(a + j0));
 const __m256i T1 = _mm256_loadu_si256((__m256i *)(a + j1));
 const __m256i T2 = _mm256_loadu_si256((__m256i *)(a + j2));
 const __m256i T3 = _mm256_loadu_si256((__m256i *)(a + j3));
 const __m256i T0P1 = montgomery_add_256(T0, T1, m2, m0);
 const __m256i T2P3 = montgomery_add_256(T2, T3, m2, m0);
 const __m256i T0M1 = montgomery_mul_256(
 montgomery_sub_256(T0, T1, m2, m0), XX, r, m1);
 const __m256i T2M3 = montgomery_mul_256(
 montgomery_sub_256(T2, T3, m2, m0), YY, r, m1);
 _mm256_storeu_si256((__m256i *)(a + j0),
 montgomery_add_256(T0P1, T2P3, m2, m0));
 _mm256_storeu_si256(
 (__m256i *)(a + j2),
 montgomery_mul_256(montgomery_sub_256(T0P1, T2P3, m2, m0), WW,
 r, m1));
 _mm256_storeu_si256((__m256i *)(a + j1),
 montgomery_add_256(T0M1, T2M3, m2, m0));
 _mm256_storeu_si256(
 (__m256i *)(a + j3),
 montgomery_mul_256(montgomery_sub_256(T0M1, T2M3, m2, m0), WW,
 r, m1));
 }
 }
 xx *= dy[__builtin_ctz(jh += 4)];
 }
 }
 u >>= 4;
 v <<= 2;
 }
 if (k & 1) {
 v = 1 << (k - 1);
 if (v < 8) {
 for (int j = 0; j < v; ++j) {
 mint ajv = a[j] - a[j + v];
 a[j] += a[j + v];
 a[j + v] = ajv;
 }
 } else {
 const __m256i m0 = _mm256_set1_epi32(0);
 const __m256i m2 = _mm256_set1_epi32(mod + mod);
 int j0 = 0;
 int j1 = v;
 for (; j0 < v; j0 += 8, j1 += 8) {
 const __m256i T0 = _mm256_loadu_si256((__m256i *)(a + j0));
 const __m256i T1 = _mm256_loadu_si256((__m256i *)(a + j1));
 __m256i naj = montgomery_add_256(T0, T1, m2, m0);
 __m256i najv = montgomery_sub_256(T0, T1, m2, m0);
 _mm256_storeu_si256((__m256i *)(a + j0), naj);
 _mm256_storeu_si256((__m256i *)(a + j1), najv);
 }
 }
 }
 if (normalize) {
 mint invn = mint(n).inverse();
 for (int i = 0; i < n; i++) a[i] *= invn;
 }
 }

 __attribute__((target("avx2"))) void inplace_multiply(
 int l1, int l2, int zero_padding = true) {
 int l = l1 + l2 - 1;
 int M = 4;
 while (M < l) M <<= 1;
 if (zero_padding) {
 for (int i = l1; i < M; i++) ntt_inner::_buf1[i] = 0;
 for (int i = l2; i < M; i++) ntt_inner::_buf2[i] = 0;
 }
 const __m256i m0 = _mm256_set1_epi32(0);
 const __m256i m1 = _mm256_set1_epi32(mod);
 const __m256i r = _mm256_set1_epi32(mint::r);
 const __m256i N2 = _mm256_set1_epi32(mint::n2);
 for (int i = 0; i < l1; i += 8) {
 __m256i a = _mm256_loadu_si256((__m256i *)(ntt_inner::_buf1 + i));
 __m256i b = montgomery_mul_256(a, N2, r, m1);
 _mm256_storeu_si256((__m256i *)(ntt_inner::_buf1 + i), b);
 }
 for (int i = 0; i < l2; i += 8) {
 __m256i a = _mm256_loadu_si256((__m256i *)(ntt_inner::_buf2 + i));
 __m256i b = montgomery_mul_256(a, N2, r, m1);
 _mm256_storeu_si256((__m256i *)(ntt_inner::_buf2 + i), b);
 }
 ntt(buf1, M);
 ntt(buf2, M);
 for (int i = 0; i < M; i += 8) {
 __m256i a = _mm256_loadu_si256((__m256i *)(ntt_inner::_buf1 + i));
 __m256i b = _mm256_loadu_si256((__m256i *)(ntt_inner::_buf2 + i));
 __m256i c = montgomery_mul_256(a, b, r, m1);
 _mm256_storeu_si256((__m256i *)(ntt_inner::_buf1 + i), c);
 }
 intt(buf1, M, false);
 const __m256i INVM = _mm256_set1_epi32((mint(M).inverse()).a);
 for (int i = 0; i < l; i += 8) {
 __m256i a = _mm256_loadu_si256((__m256i *)(ntt_inner::_buf1 + i));
 __m256i b = montgomery_mul_256(a, INVM, r, m1);
 __m256i c = my256_mulhi_epu32(my256_mullo_epu32(b, r), m1);
 __m256i d = _mm256_and_si256(_mm256_cmpgt_epi32(c, m0), m1);
 __m256i e = _mm256_sub_epi32(d, c);
 _mm256_storeu_si256((__m256i *)(ntt_inner::_buf1 + i), e);
 }
 }

 void ntt(vector<mint> &a) {
 int M = (int)a.size();
 for (int i = 0; i < M; i++) buf1[i].a = a[i].a;
 ntt(buf1, M);
 for (int i = 0; i < M; i++) a[i].a = buf1[i].a;
 }

 void intt(vector<mint> &a) {
 int M = (int)a.size();
 for (int i = 0; i < M; i++) buf1[i].a = a[i].a;
 intt(buf1, M, true);
 for (int i = 0; i < M; i++) a[i].a = buf1[i].a;
 }

 vector<mint> multiply(const vector<mint> &a, const vector<mint> &b) {
 if (a.size() == 0 && b.size() == 0) return vector<mint>{};
 int l = a.size() + b.size() - 1;
 if (min<int>(a.size(), b.size()) <= 40) {
 vector<mint> s(l);
 for (int i = 0; i < (int)a.size(); ++i)
 for (int j = 0; j < (int)b.size(); ++j) s[i + j] += a[i] * b[j];
 return s;
 }
 assert(l <= ntt_inner::SZ_FFT_BUF);
 int M = 4;
 while (M < l) M <<= 1;
 for (int i = 0; i < (int)a.size(); ++i) buf1[i].a = a[i].a;
 for (int i = (int)a.size(); i < M; ++i) buf1[i].a = 0;
 for (int i = 0; i < (int)b.size(); ++i) buf2[i].a = b[i].a;
 for (int i = (int)b.size(); i < M; ++i) buf2[i].a = 0;
 ntt(buf1, M);
 ntt(buf2, M);
 for (int i = 0; i < M; ++i)
 buf1[i].a = mint::reduce(uint64_t(buf1[i].a) * buf2[i].a);
 intt(buf1, M, false);
 vector<mint> s(l);
 mint invm = mint(M).inverse();
 for (int i = 0; i < l; ++i) s[i] = buf1[i] * invm;
 return s;
 }

 void ntt_doubling(vector<mint> &a) {
 int M = (int)a.size();
 for (int i = 0; i < M; i++) buf1[i].a = a[i].a;
 intt(buf1, M);
 mint r = 1, zeta = mint(pr).pow((mint::get_mod() - 1) / (M << 1));
 for (int i = 0; i < M; i++) buf1[i] *= r, r *= zeta;
 ntt(buf1, M);
 a.resize(2 * M);
 for (int i = 0; i < M; i++) a[M + i].a = buf1[i].a;
 }
};


template <typename mint>
void FormalPowerSeries<mint>::set_fft() {
 if (!ntt_ptr) ntt_ptr = new NTT<mint>;
}

template <typename mint>
FormalPowerSeries<mint>& FormalPowerSeries<mint>::operator*=(
 const FormalPowerSeries<mint>& r) {
 if (this->empty() || r.empty()) {
 this->clear();
 return *this;
 }
 set_fft();
 auto ret = static_cast<NTT<mint>*>(ntt_ptr)->multiply(*this, r);
 return *this = FormalPowerSeries<mint>(ret.begin(), ret.end());
}

template <typename mint>
void FormalPowerSeries<mint>::ntt() {
 set_fft();
 static_cast<NTT<mint>*>(ntt_ptr)->ntt(*this);
}

template <typename mint>
void FormalPowerSeries<mint>::intt() {
 set_fft();
 static_cast<NTT<mint>*>(ntt_ptr)->intt(*this);
}

template <typename mint>
void FormalPowerSeries<mint>::ntt_doubling() {
 set_fft();
 static_cast<NTT<mint>*>(ntt_ptr)->ntt_doubling(*this);
}

template <typename mint>
int FormalPowerSeries<mint>::ntt_pr() {
 set_fft();
 return static_cast<NTT<mint>*>(ntt_ptr)->pr;
}

template <typename mint>
FormalPowerSeries<mint> FormalPowerSeries<mint>::inv(int deg) const {
 assert((*this)[0] != mint(0));
 if (deg == -1) deg = (int)this->size();
 FormalPowerSeries<mint> res(deg);
 res[0] = {mint(1) / (*this)[0]};
 for (int d = 1; d < deg; d <<= 1) {
 FormalPowerSeries<mint> f(2 * d), g(2 * d);
 for (int j = 0; j < min((int)this->size(), 2 * d); j++) f[j] = (*this)[j];
 for (int j = 0; j < d; j++) g[j] = res[j];
 f.ntt();
 g.ntt();
 for (int j = 0; j < 2 * d; j++) f[j] *= g[j];
 f.intt();
 for (int j = 0; j < d; j++) f[j] = 0;
 f.ntt();
 for (int j = 0; j < 2 * d; j++) f[j] *= g[j];
 f.intt();
 for (int j = d; j < min(2 * d, deg); j++) res[j] = -f[j];
 }
 return res.pre(deg);
}

template <typename mint>
FormalPowerSeries<mint> FormalPowerSeries<mint>::exp(int deg) const {
 using fps = FormalPowerSeries<mint>;
 assert((*this).size() == 0 || (*this)[0] == mint(0));
 if (deg == -1) deg = this->size();

 fps inv;
 inv.reserve(deg + 1);
 inv.push_back(mint(0));
 inv.push_back(mint(1));

 auto inplace_integral = [&](fps& F) -> void {
 const int n = (int)F.size();
 auto mod = mint::get_mod();
 while ((int)inv.size() <= n) {
 int i = inv.size();
 inv.push_back((-inv[mod % i]) * (mod / i));
 }
 F.insert(begin(F), mint(0));
 for (int i = 1; i <= n; i++) F[i] *= inv[i];
 };

 auto inplace_diff = [](fps& F) -> void {
 if (F.empty()) return;
 F.erase(begin(F));
 mint coeff = 1, one = 1;
 for (int i = 0; i < (int)F.size(); i++) {
 F[i] *= coeff;
 coeff += one;
 }
 };

 fps b{1, 1 < (int)this->size() ? (*this)[1] : 0}, c{1}, z1, z2{1, 1};
 for (int m = 2; m < deg; m *= 2) {
 auto y = b;
 y.resize(2 * m);
 y.ntt();
 z1 = z2;
 fps z(m);
 for (int i = 0; i < m; ++i) z[i] = y[i] * z1[i];
 z.intt();
 fill(begin(z), begin(z) + m / 2, mint(0));
 z.ntt();
 for (int i = 0; i < m; ++i) z[i] *= -z1[i];
 z.intt();
 c.insert(end(c), begin(z) + m / 2, end(z));
 z2 = c;
 z2.resize(2 * m);
 z2.ntt();
 fps x(begin(*this), begin(*this) + min<int>(this->size(), m));
 x.resize(m);
 inplace_diff(x);
 x.push_back(mint(0));
 x.ntt();
 for (int i = 0; i < m; ++i) x[i] *= y[i];
 x.intt();
 x -= b.diff();
 x.resize(2 * m);
 for (int i = 0; i < m - 1; ++i) x[m + i] = x[i], x[i] = mint(0);
 x.ntt();
 for (int i = 0; i < 2 * m; ++i) x[i] *= z2[i];
 x.intt();
 x.pop_back();
 inplace_integral(x);
 for (int i = m; i < min<int>(this->size(), 2 * m); ++i) x[i] += (*this)[i];
 fill(begin(x), begin(x) + m, mint(0));
 x.ntt();
 for (int i = 0; i < 2 * m; ++i) x[i] *= y[i];
 x.intt();
 b.insert(end(b), begin(x) + m, end(x));
 }
 return fps{begin(b), begin(b) + deg};
}









template <uint32_t mod>
struct LazyMontgomeryModInt {
 using mint = LazyMontgomeryModInt;
 using i32 = int32_t;
 using u32 = uint32_t;
 using u64 = uint64_t;

 static constexpr u32 get_r() {
 u32 ret = mod;
 for (i32 i = 0; i < 4; ++i) ret *= 2 - mod * ret;
 return ret;
 }

 static constexpr u32 r = get_r();
 static constexpr u32 n2 = -u64(mod) % mod;
 static_assert(mod < (1 << 30), "invalid, mod >= 2 ^ 30");
 static_assert((mod & 1) == 1, "invalid, mod % 2 == 0");
 static_assert(r * mod == 1, "this code has bugs.");

 u32 a;

 constexpr LazyMontgomeryModInt() : a(0) {}
 constexpr LazyMontgomeryModInt(const int64_t &b)
 : a(reduce(u64(b % mod + mod) * n2)){};

 static constexpr u32 reduce(const u64 &b) {
 return (b + u64(u32(b) * u32(-r)) * mod) >> 32;
 }

 constexpr mint &operator+=(const mint &b) {
 if (i32(a += b.a - 2 * mod) < 0) a += 2 * mod;
 return *this;
 }

 constexpr mint &operator-=(const mint &b) {
 if (i32(a -= b.a) < 0) a += 2 * mod;
 return *this;
 }

 constexpr mint &operator*=(const mint &b) {
 a = reduce(u64(a) * b.a);
 return *this;
 }

 constexpr mint &operator/=(const mint &b) {
 *this *= b.inverse();
 return *this;
 }

 constexpr mint operator+(const mint &b) const { return mint(*this) += b; }
 constexpr mint operator-(const mint &b) const { return mint(*this) -= b; }
 constexpr mint operator*(const mint &b) const { return mint(*this) *= b; }
 constexpr mint operator/(const mint &b) const { return mint(*this) /= b; }
 constexpr bool operator==(const mint &b) const {
 return (a >= mod ? a - mod : a) == (b.a >= mod ? b.a - mod : b.a);
 }
 constexpr bool operator!=(const mint &b) const {
 return (a >= mod ? a - mod : a) != (b.a >= mod ? b.a - mod : b.a);
 }
 constexpr mint operator-() const { return mint() - mint(*this); }
 constexpr mint operator+() const { return mint(*this); }

 constexpr mint pow(u64 n) const {
 mint ret(1), mul(*this);
 while (n > 0) {
 if (n & 1) ret *= mul;
 mul *= mul;
 n >>= 1;
 }
 return ret;
 }

 constexpr mint inverse() const {
 int x = get(), y = mod, u = 1, v = 0, t = 0, tmp = 0;
 while (y > 0) {
 t = x / y;
 x -= t * y, u -= t * v;
 tmp = x, x = y, y = tmp;
 tmp = u, u = v, v = tmp;
 }
 return mint{u};
 }

 friend ostream &operator<<(ostream &os, const mint &b) {
 return os << b.get();
 }

 friend istream &operator>>(istream &is, mint &b) {
 int64_t t;
 is >> t;
 b = LazyMontgomeryModInt<mod>(t);
 return (is);
 }

 constexpr u32 get() const {
 u32 ret = reduce(a);
 return ret >= mod ? ret - mod : ret;
 }

 static constexpr u32 get_mod() { return mod; }
};






using namespace Nyaan;
using mint = LazyMontgomeryModInt<998244353>;

using vm = vector<mint>;
using vvm = vector<vm>;
Binomial<mint> C;
using fps = FormalPowerSeries<mint>;
using namespace Nyaan;

using frac = fps_fraction<fps>;
using Mat = Matrix<fps>;

void shrink(frac& f) { f.p.shrink(), f.q.shrink(); }
void shrink(Mat& m) { rep(i, 2) rep(j, 2) m[i][j].shrink(); }

vvi g2;

frac light(int c);

V<frac> heavy(int c) {
 if (g2[c].empty()) return V<frac>{frac{fps{0}, fps{1}}};
 vector<frac> v, res = heavy(g2[c][0]);
 v.push_back(frac{fps{0}, fps{1}});
 rep1(i, sz(g2[c]) - 1) v.push_back(light(g2[c][i]));
 while (v.size() >= 2u) {
 vector<frac> w;
 for (int i = 0; i + 1 < (int)v.size(); i += 2) w.push_back(v[i] + v[i + 1]);
 if (v.size() & 1) w.push_back(v.back());
 each(x, w) shrink(x);
 swap(v, w);
 }
 res.push_back(v.back());
 return res;
}

frac light(int c) {
 auto v = heavy(c);

 Mat I(2);
 I[0][0] = I[1][1] = fps{1};
 Mat oya(2);
 oya[0][1] = fps{1};
 oya[1][0] = fps{0, 0, -1};
 oya[1][1] = fps{1, -1};

 V<Mat> ms{I, I};
 each(f, v) {
 Mat m(2);
 m[0][0] = m[1][1] = f.q;
 m[0][1] = f.p;
 ms.push_back(m);
 ms.push_back(oya);
 }
 reverse(all(ms));
 while (sz(ms) >= 2) {
 V<Mat> nx;
 for (int i = 0; i + 1 < sz(ms); i += 2) {
 nx.push_back(ms[i] * ms[i + 1]);
 }
 if (sz(ms) % 2) nx.push_back(ms.back());
 each(x, nx) shrink(x);
 ms = nx;
 }
 auto m = ms[0];
 fps p = m[0][1];
 fps q = m[1][1];
 return frac{p, q};
}

mint calc(int N, int M, int S, int T, vvi g) {
 auto path = Path(g, S, T);
 g = rooted_tree(g, T);
 HeavyLightDecomposition hld{g};
 if (sz(path) - 1 > M) return 0;

 {
 set<int> pathset;
 each(p, path) pathset.insert(p);
 g2.resize(N);
 rep(i, N) each(j, g[i]) {
 if (pathset.count(j) == 0) g2[i].push_back(j);
 }
 }

 
 V<frac> dp(sz(path));
 rep(j, sz(path)) {
 int c = path[j];
 V<frac> v;
 v.push_back(frac{fps{0}, fps{1}});
 rep(i, sz(g2[c])) v.push_back(light(g2[c][i]));
 while (v.size() >= 2u) {
 vector<frac> w;
 for (int i = 0; i + 1 < (int)v.size(); i += 2)
 w.push_back(v[i] + v[i + 1]);
 if (v.size() & 1) w.push_back(v.back());
 each(x, w) shrink(x);
 swap(v, w);
 }
 dp[j] = v[0];
 shrink(dp[j]);
 }

 














 trc2("dp ok");
 {
 ll dsum = 0;
 each(x, dp) dsum += sz(x.p) + sz(x.q);
 trc2(dsum);
 }

 frac fs = fps{1, -1} - dp[0] * fps{0, 0, 1};
 swap(fs.p, fs.q);

 V<Mat> ms;
 {
 Mat m(2);
 m[0][0] = m[1][1] = fps{1};
 ms.push_back(m);
 }
 rep1(i, sz(path) - 1) {
 Mat m(2);
 m[0][1] = dp[i].q;
 m[1][0] = dp[i].q * fps{0, 0, -1};
 m[1][1] = dp[i].q * fps{1, -1} - dp[i].p * fps{0, 0, 1};
 ms.push_back(m);
 }
 reverse(all(ms));
 while (sz(ms) >= 2) {
 V<Mat> nx;
 for (int i = 0; i + 1 < sz(ms); i += 2) nx.push_back(ms[i] * ms[i + 1]);
 if (sz(ms) % 2) nx.push_back(ms.back());
 each(x, nx) shrink(x);
 ms = nx;
 }

 auto m = ms[0];
 V<fps> ftps;
 ftps.push_back(fs.p);
 rep1(i, sz(path) - 1) ftps.push_back(dp[i].q);
 fps ftp = Pi(ftps);
 fps ftq = fs.p * m[1][0] + fs.q * m[1][1];
 trc2("ft ok");
 trc2(sz(ftp), sz(ftq));
 return LinearRecurrence(M - (sz(path) - 1), ftq, ftp);
}

void q() {
 ini(N, M, S, T);
 auto g = graph(N);
 S--, T--;
 
 out(calc(N, M, S, T, g));
}

void Nyaan::solve() {
 int t = 1;
 
 while (t--) q();
}
yukicoder

結果

テストケース

ソースコード
