QOJ.ac
QOJ
| ID | Problem | Submitter | Result | Time | Memory | Language | File size | Submit time | Judge time |
|---|---|---|---|---|---|---|---|---|---|
| #379492 | #8568. Expected Diameter | ucup-team180# | AC ✓ | 2046ms | 52880kb | C++17 | 86.7kb | 2024-04-06 17:40:30 | 2024-04-06 17:40:32 |
Judging History
answer
#pragma region Macros
#ifdef noimi
#pragma comment(linker, "/stack:256000000")
#include "my_template.hpp"
#else
// #pragma GCC target("avx2")
#pragma GCC optimize("Ofast")
#pragma GCC optimize("unroll-loops")
#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 <immintrin.h>
#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 <utility>
#include <variant>
#ifdef noimi
#define oj_local(a, b) b
#else
#define oj_local(a, b) a
#endif
#define LOCAL if(oj_local(0, 1))
#define OJ if(oj_local(1, 0))
using namespace std;
using ll = long long;
using ull = unsigned long long int;
using i128 = __int128_t;
using pii = pair<int, int>;
using pll = pair<ll, ll>;
using ld = long double;
template <typename T> using vc = vector<T>;
template <typename T> using vvc = vector<vc<T>>;
template <typename T> using vvvc = vector<vvc<T>>;
using vi = vc<int>;
using vl = vc<ll>;
using vpi = vc<pii>;
using vpl = vc<pll>;
template <class T> using pq = priority_queue<T>;
template <class T> using pqg = priority_queue<T, vector<T>, greater<T>>;
template <typename T> int si(const T &x) { return x.size(); }
template <class T, class S> inline bool chmax(T &a, const S &b) { return (a < b ? a = b, 1 : 0); }
template <class T, class S> inline bool chmin(T &a, const S &b) { return (a > b ? a = b, 1 : 0); }
vi iota(int n) {
vi a(n);
return iota(a.begin(), a.end(), 0), a;
}
template <typename T> vi iota(const vector<T> &a, bool greater = false) {
vi res(a.size());
iota(res.begin(), res.end(), 0);
sort(res.begin(), res.end(), [&](int i, int j) {
if(greater) return a[i] > a[j];
return a[i] < a[j];
});
return res;
}
// macros
#define overload5(a, b, c, d, e, name, ...) name
#define overload4(a, b, c, d, name, ...) name
#define endl '\n'
#define REP0(n) for(ll jidlsjf = 0; jidlsjf < n; ++jidlsjf)
#define REP1(i, n) for(ll i = 0; i < (n); ++i)
#define REP2(i, a, b) for(ll i = (a); i < (b); ++i)
#define REP3(i, a, b, c) for(ll i = (a); i < (b); i += (c))
#define rep(...) overload4(__VA_ARGS__, REP3, REP2, REP1, REP0)(__VA_ARGS__)
#define per0(n) for(int jidlsjf = 0; jidlsjf < (n); ++jidlsjf)
#define per1(i, n) for(ll i = (n)-1; i >= 0; --i)
#define per2(i, a, b) for(ll i = (a)-1; i >= b; --i)
#define per3(i, a, b, c) for(ll i = (a)-1; i >= (b); i -= (c))
#define per(...) overload4(__VA_ARGS__, per3, per2, per1, per0)(__VA_ARGS__)
#define fore0(a) rep(a.size())
#define fore1(i, a) for(auto &&i : a)
#define fore2(a, b, v) for(auto &&[a, b] : v)
#define fore3(a, b, c, v) for(auto &&[a, b, c] : v)
#define fore4(a, b, c, d, v) for(auto &&[a, b, c, d] : v)
#define fore(...) overload5(__VA_ARGS__, fore4, fore3, fore2, fore1, fore0)(__VA_ARGS__)
#define setbits(j, n) for(ll iiiii = (n), j = lowbit(iiiii); iiiii; iiiii ^= 1 << j, j = lowbit(iiiii))
#define perm(v) for(bool permrepflag = true; (permrepflag ? exchange(permrepflag, false) : next_permutation(all(v)));)
#define fi first
#define se second
#define pb push_back
#define ppb pop_back
#define ppf pop_front
#define eb emplace_back
#define drop(s) cout << #s << endl, exit(0)
#define si(c) (int)(c).size()
#define lb(c, x) distance((c).begin(), lower_bound(all(c), (x)))
#define lbg(c, x) distance((c).begin(), lower_bound(all(c), (x), greater{}))
#define ub(c, x) distance((c).begin(), upper_bound(all(c), (x)))
#define ubg(c, x) distance((c).begin(), upper_bound(all(c), (x), greater{}))
#define rng(v, l, r) v.begin() + (l), v.begin() + (r)
#define all(c) begin(c), end(c)
#define rall(c) rbegin(c), rend(c)
#define SORT(v) sort(all(v))
#define REV(v) reverse(all(v))
#define UNIQUE(x) SORT(x), x.erase(unique(all(x)), x.end())
template <typename T = ll, typename S> T SUM(const S &v) { return accumulate(all(v), T(0)); }
#define MIN(v) *min_element(all(v))
#define MAX(v) *max_element(all(v))
#define overload2(_1, _2, name, ...) name
#define vec(type, name, ...) vector<type> name(__VA_ARGS__)
#define vv(type, name, h, ...) vector<vector<type>> name(h, vector<type>(__VA_ARGS__))
#define vvv(type, name, h, w, ...) vector<vector<vector<type>>> name(h, vector<vector<type>>(w, vector<type>(__VA_ARGS__)))
#define vvvv(type, name, a, b, c, ...) \
vector<vector<vector<vector<type>>>> name(a, vector<vector<vector<type>>>(b, vector<vector<type>>(c, vector<type>(__VA_ARGS__))))
constexpr pii dx4[4] = {pii{1, 0}, pii{0, 1}, pii{-1, 0}, pii{0, -1}};
constexpr pii dx8[8] = {{1, 0}, {1, 1}, {0, 1}, {-1, 1}, {-1, 0}, {-1, -1}, {0, -1}, {1, -1}};
namespace yesno_impl {
const string YESNO[2] = {"NO", "YES"};
const string YesNo[2] = {"No", "Yes"};
const string yesno[2] = {"no", "yes"};
const string firstsecond[2] = {"second", "first"};
const string FirstSecond[2] = {"Second", "First"};
const string possiblestr[2] = {"impossible", "possible"};
const string Possiblestr[2] = {"Impossible", "Possible"};
void YES(bool t = 1) { cout << YESNO[t] << endl; }
void NO(bool t = 1) { YES(!t); }
void Yes(bool t = 1) { cout << YesNo[t] << endl; }
void No(bool t = 1) { Yes(!t); }
void yes(bool t = 1) { cout << yesno[t] << endl; }
void no(bool t = 1) { yes(!t); }
void first(bool t = 1) { cout << firstsecond[t] << endl; }
void First(bool t = 1) { cout << FirstSecond[t] << endl; }
void possible(bool t = 1) { cout << possiblestr[t] << endl; }
void Possible(bool t = 1) { cout << Possiblestr[t] << endl; }
}; // namespace yesno_impl
using namespace yesno_impl;
#define INT(...) \
int __VA_ARGS__; \
IN(__VA_ARGS__)
#define INTd(...) \
int __VA_ARGS__; \
IN2(__VA_ARGS__)
#define LL(...) \
ll __VA_ARGS__; \
IN(__VA_ARGS__)
#define LLd(...) \
ll __VA_ARGS__; \
IN2(__VA_ARGS__)
#define STR(...) \
string __VA_ARGS__; \
IN(__VA_ARGS__)
#define CHR(...) \
char __VA_ARGS__; \
IN(__VA_ARGS__)
#define DBL(...) \
double __VA_ARGS__; \
IN(__VA_ARGS__)
#define VEC(type, name, size) \
vector<type> name(size); \
IN(name)
#define VECd(type, name, size) \
vector<type> name(size); \
IN2(name)
#define VEC2(type, name1, name2, size) \
vector<type> name1(size), name2(size); \
for(int i = 0; i < size; i++) IN(name1[i], name2[i])
#define VEC2d(type, name1, name2, size) \
vector<type> name1(size), name2(size); \
for(int i = 0; i < size; i++) IN2(name1[i], name2[i])
#define VEC3(type, name1, name2, name3, size) \
vector<type> name1(size), name2(size), name3(size); \
for(int i = 0; i < size; i++) IN(name1[i], name2[i], name3[i])
#define VEC3d(type, name1, name2, name3, size) \
vector<type> name1(size), name2(size), name3(size); \
for(int i = 0; i < size; i++) IN2(name1[i], name2[i], name3[i])
#define VEC4(type, name1, name2, name3, name4, size) \
vector<type> name1(size), name2(size), name3(size), name4(size); \
for(int i = 0; i < size; i++) IN(name1[i], name2[i], name3[i], name4[i]);
#define VEC4d(type, name1, name2, name3, name4, size) \
vector<type> name1(size), name2(size), name3(size), name4(size); \
for(int i = 0; i < size; i++) IN2(name1[i], name2[i], name3[i], name4[i]);
#define VV(type, name, h, w) \
vector<vector<type>> name(h, vector<type>(w)); \
IN(name)
#define VVd(type, name, h, w) \
vector<vector<type>> name(h, vector<type>(w)); \
IN2(name)
int scan() { return getchar(); }
void scan(int &a) { cin >> a; }
void scan(long long &a) { cin >> a; }
void scan(char &a) { cin >> a; }
void scan(double &a) { cin >> a; }
void scan(string &a) { cin >> a; }
template <class T, class S> void scan(pair<T, S> &p) { scan(p.first), scan(p.second); }
template <class T> void scan(vector<T> &);
template <class T> void scan(vector<T> &a) {
for(auto &i : a) scan(i);
}
template <class T> void scan(T &a) { cin >> a; }
void IN() {}
void IN2() {}
template <class Head, class... Tail> void IN(Head &head, Tail &...tail) {
scan(head);
IN(tail...);
}
template <class Head, class... Tail> void IN2(Head &head, Tail &...tail) {
scan(head);
--head;
IN2(tail...);
}
template <int p = -1> void pat() {}
template <int p = -1, class Head, class... Tail> void pat(Head &h, Tail &...tail) {
h += p;
pat<p>(tail...);
}
template <typename T, typename S> T ceil(T x, S y) {
assert(y);
return (y < 0 ? ceil(-x, -y) : (x > 0 ? (x + y - 1) / y : x / y));
}
template <typename T, typename S> T floor(T x, S y) {
assert(y);
return (y < 0 ? floor(-x, -y) : (x > 0 ? x / y : x / y - (x % y == 0 ? 0 : 1)));
}
template <typename T, typename S, typename U> U bigmul(const T &x, const S &y, const U &lim) { // clamp(x * y, -lim, lim)
if(x < 0 and y < 0) return bigmul(-x, -y, lim);
if(x < 0) return -bigmul(-x, y, lim);
if(y < 0) return -bigmul(x, -y, lim);
return y == 0 or x <= lim / y ? x * y : lim;
}
template <class T> T POW(T x, int n) {
T res = 1;
for(; n; n >>= 1, x *= x)
if(n & 1) res *= x;
return res;
}
template <class T, class S> T POW(T x, S n, const ll &mod) {
T res = 1;
x %= mod;
for(; n; n >>= 1, x = x * x % mod)
if(n & 1) res = res * x % mod;
return res;
}
vector<pll> factor(ll x) {
vector<pll> ans;
for(ll i = 2; i * i <= x; i++)
if(x % i == 0) {
ans.push_back({i, 1});
while((x /= i) % i == 0) ans.back().second++;
}
if(x != 1) ans.push_back({x, 1});
return ans;
}
template <class T> vector<T> divisor(T x) {
vector<T> ans;
for(T i = 1; i * i <= x; i++)
if(x % i == 0) {
ans.pb(i);
if(i * i != x) ans.pb(x / i);
}
return ans;
}
template <typename T> void zip(vector<T> &x) {
vector<T> y = x;
UNIQUE(y);
for(int i = 0; i < x.size(); ++i) { x[i] = lb(y, x[i]); }
}
template <class S> void fold_in(vector<S> &v) {}
template <typename Head, typename... Tail, class S> void fold_in(vector<S> &v, Head &&a, Tail &&...tail) {
for(auto e : a) v.emplace_back(e);
fold_in(v, tail...);
}
template <class S> void renumber(vector<S> &v) {}
template <typename Head, typename... Tail, class S> void renumber(vector<S> &v, Head &&a, Tail &&...tail) {
for(auto &&e : a) e = lb(v, e);
renumber(v, tail...);
}
template <class S, class... Args> vector<S> zip(vector<S> &head, Args &&...args) {
vector<S> v;
fold_in(v, head, args...);
sort(all(v)), v.erase(unique(all(v)), v.end());
renumber(v, head, args...);
return v;
}
template <typename S> void rearrange(const vector<S> &id) {}
template <typename S, typename T> void rearrange_exec(const vector<S> &id, vector<T> &v) {
vector<T> w(v.size());
rep(i, si(id)) w[i] = v[id[i]];
v.swap(w);
}
// 並び替える順番, 並び替える vector 達
template <typename S, typename Head, typename... Tail> void rearrange(const vector<S> &id, Head &a, Tail &...tail) {
rearrange_exec(id, a);
rearrange(id, tail...);
}
template <typename T> vector<T> RUI(const vector<T> &v) {
vector<T> res(v.size() + 1);
for(int i = 0; i < v.size(); i++) res[i + 1] = res[i] + v[i];
return res;
}
template <typename T> void zeta_supersetsum(vector<T> &f) {
int n = f.size();
for(int i = 1; i < n; i <<= 1) rep(b, n) if(!(i & b)) f[b] += f[b | i];
}
template <typename T> void zeta_subsetsum(vector<T> &f) {
int n = f.size();
for(int i = 1; i < n; i <<= 1) rep(b, n) if(!(i & b)) f[b | i] += f[b];
}
template <typename T> void mobius_subset(vector<T> &f) {
int n = f.size();
for(int i = 1; i < n; i <<= 1) rep(b, n) if(!(i & b)) f[b] -= f[b | i];
}
template <typename T> void mobius_superset(vector<T> &f) {
int n = f.size();
for(int i = 1; i < n; i <<= 1) rep(b, n) if(!(i & b)) f[b | i] -= f[b];
}
// 反時計周りに 90 度回転
template <typename T> void rot(vector<vector<T>> &v) {
if(empty(v)) return;
int n = v.size(), m = v[0].size();
vector<vector<T>> res(m, vector<T>(n));
rep(i, n) rep(j, m) res[m - 1 - j][i] = v[i][j];
v.swap(res);
}
vector<int> counter(const vector<int> &v, int max_num = -1) {
if(max_num == -1) max_num = MAX(v);
vector<int> res(max_num + 1);
fore(e, v) res[e]++;
return res;
}
// x in [l, r)
template <class T, class S> bool inc(const T &x, const S &l, const S &r) { return l <= x and x < r; }
template <class T, class S> bool inc(const T &x, const pair<S, S> &p) { return p.first <= x and x < p.second; }
// 便利関数
constexpr ll ten(int n) { return n == 0 ? 1 : ten(n - 1) * 10; }
constexpr ll tri(ll n) { return n * (n + 1) / 2; }
// l + ... + r
constexpr ll tri(ll l, ll r) { return (l + r) * (r - l + 1) / 2; }
ll max(int x, ll y) { return max((ll)x, y); }
ll max(ll x, int y) { return max(x, (ll)y); }
int min(int x, ll y) { return min((ll)x, y); }
int min(ll x, int y) { return min(x, (ll)y); }
// bit 演算系
#define bit(i) (1LL << i) // (1 << i)
#define test(b, i) (b >> i & 1) // b の i bit 目が立っているか
ll pow2(int i) { return 1LL << i; }
int topbit(signed t) { return t == 0 ? -1 : 31 - __builtin_clz(t); }
int topbit(ll t) { return t == 0 ? -1 : 63 - __builtin_clzll(t); }
// int lowbit(signed a) { return a == 0 ? 32 : __builtin_ctz(a); }
int lowbit(ull a) { return a == 0 ? 64 : __builtin_ctzll(a); }
// int allbit(int n) { return (1 << n) - 1; }
constexpr ll mask(int n) { return (1LL << n) - 1; }
// int popcount(signed t) { return __builtin_popcount(t); }
// int popcount(ll t) { return __builtin_popcountll(t); }
int popcount(uint64_t t) { return __builtin_popcountll(t); }
static inline uint64_t popcount64(uint64_t x) {
uint64_t m1 = 0x5555555555555555ll;
uint64_t m2 = 0x3333333333333333ll;
uint64_t m4 = 0x0F0F0F0F0F0F0F0Fll;
uint64_t h01 = 0x0101010101010101ll;
x -= (x >> 1) & m1;
x = (x & m2) + ((x >> 2) & m2);
x = (x + (x >> 4)) & m4;
return (x * h01) >> 56;
}
bool ispow2(int i) { return i && (i & -i) == i; }
ll rnd(ll l, ll r) { //[l, r)
#ifdef noimi
static mt19937_64 gen;
#else
static mt19937_64 gen(chrono::steady_clock::now().time_since_epoch().count());
#endif
return uniform_int_distribution<ll>(l, r - 1)(gen);
}
ll rnd(ll n) { return rnd(0, n); }
template <class t> void random_shuffle(vc<t> &a) { rep(i, si(a)) swap(a[i], a[rnd(0, i + 1)]); }
int in() {
int x;
cin >> x;
return x;
}
ll lin() {
unsigned long long x;
cin >> x;
return x;
}
template <class T, class S> pair<T, S> operator-(const pair<T, S> &x) { return pair<T, S>(-x.first, -x.second); }
template <class T, class S> pair<T, S> operator-(const pair<T, S> &x, const pair<T, S> &y) { return pair<T, S>(x.fi - y.fi, x.se - y.se); }
template <class T, class S> pair<T, S> operator+(const pair<T, S> &x, const pair<T, S> &y) { return pair<T, S>(x.fi + y.fi, x.se + y.se); }
template <class T> pair<T, T> operator&(const pair<T, T> &l, const pair<T, T> &r) { return pair<T, T>(max(l.fi, r.fi), min(l.se, r.se)); }
template <class T, class S> pair<T, S> operator+=(pair<T, S> &l, const pair<T, S> &r) { return l = l + r; }
template <class T, class S> pair<T, S> operator-=(pair<T, S> &l, const pair<T, S> &r) { return l = l - r; }
template <class T> bool intersect(const pair<T, T> &l, const pair<T, T> &r) { return (l.se < r.se ? r.fi < l.se : l.fi < r.se); }
template <class T> vector<T> &operator++(vector<T> &v) {
fore(e, v) e++;
return v;
}
template <class T> vector<T> operator++(vector<T> &v, int) {
auto res = v;
fore(e, v) e++;
return res;
}
template <class T> vector<T> &operator--(vector<T> &v) {
fore(e, v) e--;
return v;
}
template <class T> vector<T> operator--(vector<T> &v, int) {
auto res = v;
fore(e, v) e--;
return res;
}
template <class T> void connect(vector<T> &l, const vector<T> &r) { fore(e, r) l.eb(e); }
template <class T> vector<T> operator+(const vector<T> &l, const vector<T> &r) {
vector<T> res(max(si(l), si(r)));
rep(i, si(l)) res[i] += l[i];
rep(i, si(r)) res[i] += r[i];
return res;
}
template <class T> vector<T> operator-(const vector<T> &l, const vector<T> &r) {
vector<T> res(max(si(l), si(r)));
rep(i, si(l)) res[i] += l[i];
rep(i, si(r)) res[i] -= r[i];
return res;
}
template <class T> vector<T> &operator+=(const vector<T> &l, const vector<T> &r) {
if(si(l) < si(r)) l.resize(si(r));
rep(i, si(r)) l[i] += r[i];
return l;
}
template <class T> vector<T> &operator-=(const vector<T> &l, const vector<T> &r) {
if(si(l) < si(r)) l.resize(si(r));
rep(i, si(r)) l[i] -= r[i];
return l;
}
template <class T> vector<T> &operator+=(vector<T> &v, const T &x) {
fore(e, v) e += x;
return v;
}
template <class T> vector<T> &operator-=(vector<T> &v, const T &x) {
fore(e, v) e -= x;
return v;
}
template <typename T> struct edge {
int from, to;
T cost;
int id;
edge(int to, T cost) : from(-1), to(to), cost(cost) {}
edge(int from, int to, T cost) : from(from), to(to), cost(cost) {}
edge(int from, int to, T cost, int id) : from(from), to(to), cost(cost), id(id) {}
constexpr bool operator<(const edge<T> &rhs) const noexcept { return cost < rhs.cost; }
edge &operator=(const int &x) {
to = x;
return *this;
}
operator int() const { return to; }
friend ostream operator<<(ostream &os, const edge &e) { return os << e.to; }
};
template <typename T> using Edges = vector<edge<T>>;
template <typename T = int> Edges<T> read_edges(int m, bool weighted = false) {
Edges<T> res;
res.reserve(m);
for(int i = 0; i < m; i++) {
int u, v, c = 0;
scan(u), scan(v), u--, v--;
if(weighted) scan(c);
res.eb(u, v, c, i);
}
return res;
}
using Tree = vector<vector<int>>;
using Graph = vector<vector<int>>;
template <class T> using Wgraph = vector<vector<edge<T>>>;
Graph getG(int n, int m = -1, bool directed = false, int margin = 1) {
Tree res(n);
if(m == -1) m = n - 1;
while(m--) {
int a, b;
cin >> a >> b;
a -= margin, b -= margin;
res[a].emplace_back(b);
if(!directed) res[b].emplace_back(a);
}
return res;
}
Graph getTreeFromPar(int n, int margin = 1) {
Graph res(n);
for(int i = 1; i < n; i++) {
int a;
cin >> a;
res[a - margin].emplace_back(i);
}
return res;
}
template <class T> Wgraph<T> getWg(int n, int m = -1, bool directed = false, int margin = 1) {
Wgraph<T> res(n);
if(m == -1) m = n - 1;
while(m--) {
int a, b;
T c;
scan(a), scan(b), scan(c);
a -= margin, b -= margin;
res[a].emplace_back(b, c);
if(!directed) res[b].emplace_back(a, c);
}
return res;
}
void add(Graph &G, int x, int y) { G[x].eb(y), G[y].eb(x); }
template <class S, class T> void add(Wgraph<S> &G, int x, int y, T c) { G[x].eb(y, c), G[y].eb(x, c); }
#define TEST \
INT(testcases); \
while(testcases--)
i128 abs(const i128 &x) { return x > 0 ? x : -x; }
istream &operator>>(istream &is, i128 &v) {
string s;
is >> s;
v = 0;
for(int i = 0; i < (int)s.size(); i++) {
if(isdigit(s[i])) { v = v * 10 + s[i] - '0'; }
}
if(s[0] == '-') { v *= -1; }
return is;
}
ostream &operator<<(ostream &os, const i128 &v) {
if(v == 0) { return (os << "0"); }
i128 num = v;
if(v < 0) {
os << '-';
num = -num;
}
string s;
for(; num > 0; num /= 10) { s.push_back((char)(num % 10) + '0'); }
reverse(s.begin(), s.end());
return (os << s);
}
namespace aux {
template <typename T, unsigned N, unsigned L> struct tp {
static void output(std::ostream &os, const T &v) {
os << std::get<N>(v) << (&os == &cerr ? ", " : " ");
tp<T, N + 1, L>::output(os, v);
}
};
template <typename T, unsigned N> struct tp<T, N, N> {
static void output(std::ostream &os, const T &v) { os << std::get<N>(v); }
};
} // namespace aux
template <typename... Ts> std::ostream &operator<<(std::ostream &os, const std::tuple<Ts...> &t) {
if(&os == &cerr) { os << '('; }
aux::tp<std::tuple<Ts...>, 0, sizeof...(Ts) - 1>::output(os, t);
if(&os == &cerr) { os << ')'; }
return os;
}
template <typename T, typename S, typename U> std::ostream &operator<<(std::ostream &os, const priority_queue<T, S, U> &_pq) {
auto pq = _pq;
vector<T> res;
while(!empty(pq)) res.emplace_back(pq.top()), pq.pop();
return os << res;
}
template <class T, class S> ostream &operator<<(ostream &os, const pair<T, S> &p) {
if(&os == &cerr) { return os << "(" << p.first << ", " << p.second << ")"; }
return os << p.first << " " << p.second;
}
template <class Ch, class Tr, class Container> std::basic_ostream<Ch, Tr> &operator<<(std::basic_ostream<Ch, Tr> &os, const Container &x) {
bool f = true;
if(&os == &cerr) os << "[";
for(auto &y : x) {
if(&os == &cerr)
os << (f ? "" : ", ") << y;
else
os << (f ? "" : " ") << y;
f = false;
}
if(&os == &cerr) os << "]";
return os;
}
#define dump(...) static_cast<void>(0)
#define dbg(...) static_cast<void>(0)
void OUT() { cout << endl; }
template <class Head, class... Tail> void OUT(const Head &head, const Tail &...tail) {
cout << head;
if(sizeof...(tail)) cout << ' ';
OUT(tail...);
}
template <typename T> static constexpr T inf = numeric_limits<T>::max() / 2;
template <class T, class S> constexpr pair<T, S> inf<pair<T, S>> = {inf<T>, inf<S>};
template <class T> void OUT2(const T &t, T INF = inf<T>, T res = -1) { OUT(t != INF ? t : res); }
template <class T> void OUT2(vector<T> &v, T INF = inf<T>, T res = -1) {
fore(e, v) if(e == INF) e = res;
OUT(v);
fore(e, v) if(e == res) e = INF;
}
template <class F> struct REC {
F f;
REC(F &&f_) : f(forward<F>(f_)) {}
template <class... Args> auto operator()(Args &&...args) const { return f(*this, forward<Args>(args)...); }
};
template <class S> vector<pair<S, int>> runLength(const vector<S> &v) {
vector<pair<S, int>> res;
for(auto &e : v) {
if(res.empty() or res.back().fi != e)
res.eb(e, 1);
else
res.back().se++;
}
return res;
}
vector<pair<char, int>> runLength(const string &v) {
vector<pair<char, int>> res;
for(auto &e : v) {
if(res.empty() or res.back().fi != e)
res.eb(e, 1);
else
res.back().se++;
}
return res;
}
struct string_converter {
char start = 0;
char type(const char &c) const { return (islower(c) ? 'a' : isupper(c) ? 'A' : isdigit(c) ? '0' : 0); }
int convert(const char &c) {
if(!start) start = type(c);
return c - start;
}
int convert(const char &c, const string &chars) { return chars.find(c); }
template <typename T> auto convert(const T &v) {
vector<decltype(convert(v[0]))> ret;
ret.reserve(size(v));
for(auto &&e : v) ret.emplace_back(convert(e));
return ret;
}
template <typename T> auto convert(const T &v, const string &chars) {
vector<decltype(convert(v[0], chars))> ret;
ret.reserve(size(v));
for(auto &&e : v) ret.emplace_back(convert(e, chars));
return ret;
}
int operator()(const char &v, char s = 0) {
start = s;
return convert(v);
}
int operator()(const char &v, const string &chars) { return convert(v, chars); }
template <typename T> auto operator()(const T &v, char s = 0) {
start = s;
return convert(v);
}
template <typename T> auto operator()(const T &v, const string &chars) { return convert(v, chars); }
} toint;
template <class T, class F> T bin_search(T ok, T ng, const F &f) {
while(abs(ok - ng) > 1) {
T mid = ok + ng >> 1;
(f(mid) ? ok : ng) = mid;
}
return ok;
}
template <class T, class F> T bin_search_double(T ok, T ng, const F &f, int iter = 80) {
while(iter--) {
T mid = (ok + ng) / 2;
(f(mid) ? ok : ng) = mid;
}
return ok;
}
struct Setup_io {
Setup_io() {
ios_base::sync_with_stdio(0), cin.tie(0), cout.tie(0);
cout << fixed << setprecision(11);
}
} setup_io;
#endif
#pragma endregion
namespace Modular998 {
#line 1 "a.cpp"
#line 2 "library/fps/berlekamp-massey.hpp"
template <typename mint> vector<mint> BerlekampMassey(const vector<mint> &s) {
const int N = (int)s.size();
vector<mint> b, c;
b.reserve(N + 1);
c.reserve(N + 1);
b.push_back(mint(1));
c.push_back(mint(1));
mint y = mint(1);
for(int ed = 1; ed <= N; ed++) {
int l = int(c.size()), m = int(b.size());
mint x = 0;
for(int i = 0; i < l; i++) x += c[i] * s[ed - l + i];
b.emplace_back(mint(0));
m++;
if(x == mint(0)) continue;
mint freq = x / y;
if(l < m) {
auto tmp = c;
c.insert(begin(c), m - l, mint(0));
for(int i = 0; i < m; i++) c[m - 1 - i] -= freq * b[m - 1 - i];
b = tmp;
y = x;
} else {
for(int i = 0; i < m; i++) c[l - 1 - i] -= freq * b[m - 1 - i];
}
}
reverse(begin(c), end(c));
return c;
}
#line 2 "library/modulo/binomial.hpp"
template <typename T> struct Binomial {
vector<T> f, g, h;
Binomial(int MAX = 0) : f(1, T(1)), g(1, T(1)), h(1, T(1)) {
while(MAX >= (int)f.size()) extend();
}
void extend() {
int n = f.size();
int m = n * 2;
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);
}
};
#line 2 "library/modint/montgomery-modint.hpp"
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(r * mod == 1, "invalid, r * mod != 1");
static_assert(mod < (1 << 30), "invalid, mod >= 2 ^ 30");
static_assert((mod & 1) == 1, "invalid, mod % 2 == 0");
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 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 { return pow(mod - 2); }
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; }
};
#line 2 "library/fps/ntt-friendly-fps.hpp"
#line 2 "library/ntt/ntt-avx2.hpp"
#line 2 "library/modint/simd-montgomery.hpp"
#include <immintrin.h>
__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);
}
#line 4 "library/ntt/ntt-avx2.hpp"
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)));
} // namespace ntt_inner
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;
}
};
#line 2 "library/fps/formal-power-series.hpp"
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 { return FPS(begin(*this), begin(*this) + min((int)this->size(), sz)); }
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)[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;
for(int i = 0; i < n; i++) {
if((*this)[i] != mint(0)) {
if(i * k > deg) return FPS(deg, mint(0));
mint rev = mint(1) / (*this)[i];
FPS ret = (((*this * rev) >> i).log(deg) * k).exp(deg) * ((*this)[i].pow(k));
ret = (ret << (i * k)).pre(deg);
if((int)ret.size() < deg) ret.resize(deg, mint(0));
return ret;
}
}
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;
/**
* @brief 多項式/形式的冪級数ライブラリ
* @docs docs/fps/formal-power-series.md
*/
#line 5 "library/fps/ntt-friendly-fps.hpp"
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};
}
/**
* @brief NTT mod用FPSライブラリ
* @docs docs/fps/ntt-friendly-fps.md
*/
#line 4 "library/fps/fast-multieval.hpp"
template <typename mint> vector<mint> FastMultiEval(const FormalPowerSeries<mint> &f, const vector<mint> &xs) {
using fps = FormalPowerSeries<mint>;
int s = xs.size();
int N = 1 << (32 - __builtin_clz((int)xs.size() - 1));
if(f.empty() || xs.empty()) return vector<mint>(s, mint(0));
vector<FormalPowerSeries<mint>> buf(2 * N);
for(int i = 0; i < N; i++) {
mint n = mint{i < s ? -xs[i] : mint(0)};
buf[i + N] = fps{n + 1, n - 1};
}
for(int i = N - 1; i > 0; i--) {
fps &g(buf[(i << 1) | 0]), &h(buf[(i << 1) | 1]);
int n = g.size();
int m = n << 1;
buf[i].reserve(m);
buf[i].resize(n);
for(int j = 0; j < n; j++) buf[i][j] = g[j] * h[j] - mint(1);
if(i != 1) {
buf[i].ntt_doubling();
for(int j = 0; j < m; j++) buf[i][j] += j < n ? mint(1) : -mint(1);
}
}
int fs = f.size();
fps root = buf[1];
root.intt();
root.push_back(1);
reverse(begin(root), end(root));
root = root.inv(fs).rev() * f;
root.erase(begin(root), begin(root) + fs - 1);
root.resize(N, mint(0));
vector<mint> ans(s);
auto calc = [&](auto rec, int i, int l, int r, fps g) -> void {
if(i >= N) {
ans[i - N] = g[0];
return;
}
int len = g.size(), m = (l + r) >> 1;
g.ntt();
fps tmp = buf[i * 2 + 1];
for(int j = 0; j < len; j++) tmp[j] *= g[j];
tmp.intt();
rec(rec, i * 2 + 0, l, m, fps{begin(tmp) + (len >> 1), end(tmp)});
if(m >= s) return;
tmp = buf[i * 2 + 0];
for(int j = 0; j < len; j++) tmp[j] *= g[j];
tmp.intt();
rec(rec, i * 2 + 1, m, r, fps{begin(tmp) + (len >> 1), end(tmp)});
};
calc(calc, 1, 0, N, root);
return ans;
}
/**
* @brief Multipoint Evaluation(高速化版)
*/
#line 2 "library/modint/arbitrary-prime-modint.hpp"
struct ArbitraryLazyMontgomeryModInt {
using mint = ArbitraryLazyMontgomeryModInt;
using i32 = int32_t;
using u32 = uint32_t;
using u64 = uint64_t;
static u32 mod;
static u32 r;
static u32 n2;
static u32 get_r() {
u32 ret = mod;
for(i32 i = 0; i < 4; ++i) ret *= 2 - mod * ret;
return ret;
}
static void set_mod(u32 m) {
assert(m < (1 << 30));
assert((m & 1) == 1);
mod = m;
n2 = -u64(m) % m;
r = get_r();
assert(r * mod == 1);
}
u32 a;
ArbitraryLazyMontgomeryModInt() : a(0) {}
ArbitraryLazyMontgomeryModInt(const int64_t &b) : a(reduce(u64(b % mod + mod) * n2)){};
static u32 reduce(const u64 &b) { return (b + u64(u32(b) * u32(-r)) * mod) >> 32; }
mint &operator+=(const mint &b) {
if(i32(a += b.a - 2 * mod) < 0) a += 2 * mod;
return *this;
}
mint &operator-=(const mint &b) {
if(i32(a -= b.a) < 0) a += 2 * mod;
return *this;
}
mint &operator*=(const mint &b) {
a = reduce(u64(a) * b.a);
return *this;
}
mint &operator/=(const mint &b) {
*this *= b.inverse();
return *this;
}
mint operator+(const mint &b) const { return mint(*this) += b; }
mint operator-(const mint &b) const { return mint(*this) -= b; }
mint operator*(const mint &b) const { return mint(*this) *= b; }
mint operator/(const mint &b) const { return mint(*this) /= b; }
bool operator==(const mint &b) const { return (a >= mod ? a - mod : a) == (b.a >= mod ? b.a - mod : b.a); }
bool operator!=(const mint &b) const { return (a >= mod ? a - mod : a) != (b.a >= mod ? b.a - mod : b.a); }
mint operator-() const { return mint() - mint(*this); }
mint pow(u64 n) const {
mint ret(1), mul(*this);
while(n > 0) {
if(n & 1) ret *= mul;
mul *= mul;
n >>= 1;
}
return ret;
}
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 = ArbitraryLazyMontgomeryModInt(t);
return (is);
}
mint inverse() const { return pow(mod - 2); }
u32 get() const {
u32 ret = reduce(a);
return ret >= mod ? ret - mod : ret;
}
static u32 get_mod() { return mod; }
};
typename ArbitraryLazyMontgomeryModInt::u32 ArbitraryLazyMontgomeryModInt::mod;
typename ArbitraryLazyMontgomeryModInt::u32 ArbitraryLazyMontgomeryModInt::r;
typename ArbitraryLazyMontgomeryModInt::u32 ArbitraryLazyMontgomeryModInt::n2;
#line 3 "library/modulo/mod-sqrt.hpp"
int64_t mod_sqrt(const int64_t &a, const int64_t &p) {
assert(0 <= a && a < p);
if(a < 2) return a;
using Mint = ArbitraryLazyMontgomeryModInt;
Mint::set_mod(p);
if(Mint(a).pow((p - 1) >> 1) != 1) return -1;
Mint b = 1, one = 1;
while(b.pow((p - 1) >> 1) == 1) b += one;
int64_t m = p - 1, e = 0;
while(m % 2 == 0) m >>= 1, e += 1;
Mint x = Mint(a).pow((m - 1) >> 1);
Mint y = Mint(a) * x * x;
x *= a;
Mint z = Mint(b).pow(m);
while(y != 1) {
int64_t j = 0;
Mint t = y;
while(t != one) {
j += 1;
t *= t;
}
z = z.pow(int64_t(1) << (e - j - 1));
x *= z;
z *= z;
y *= z;
e = j;
}
return x.get();
}
/**
* @brief mod sqrt(Tonelli-Shanks algorithm)
* @docs docs/modulo/mod-sqrt.md
*/
#line 4 "library/fps/fps-sqrt.hpp"
template <typename mint> FormalPowerSeries<mint> sqrt(const FormalPowerSeries<mint> &f, int deg = -1) {
if(deg == -1) deg = (int)f.size();
if((int)f.size() == 0) return FormalPowerSeries<mint>(deg, 0);
if(f[0] == mint(0)) {
for(int i = 1; i < (int)f.size(); i++) {
if(f[i] != mint(0)) {
if(i & 1) return {};
if(deg - i / 2 <= 0) break;
auto ret = sqrt(f >> i, deg - i / 2);
if(ret.empty()) return {};
ret = ret << (i / 2);
if((int)ret.size() < deg) ret.resize(deg, mint(0));
return ret;
}
}
return FormalPowerSeries<mint>(deg, 0);
}
int64_t sqr = mod_sqrt(f[0].get(), mint::get_mod());
if(sqr == -1) return {};
assert(sqr * sqr % mint::get_mod() == f[0].get());
FormalPowerSeries<mint> ret = {mint(sqr)};
mint inv2 = mint(2).inverse();
for(int i = 1; i < deg; i <<= 1) { ret = (ret + f.pre(i << 1) * ret.inv(i << 1)) * inv2; }
return ret.pre(deg);
}
/**
* @brief 平方根
* @docs docs/fps/fps-sqrt.md
*/
#line 2 "library/fps/kitamasa.hpp"
#line 4 "library/fps/kitamasa.hpp"
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();
// even degree of Q(x)Q(-x)
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) {
// odd degree of P(x)Q(-x)
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 {
// even degree of P(x)Q(-x)
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);
}
/**
* @brief 線形漸化式の高速計算
* @docs docs/fps/kitamasa.md
*/
#line 2 "library/fps/lagrange-interpolation-point.hpp"
#line 4 "library/fps/lagrange-interpolation-point.hpp"
// given : y(x=0) , y(x=1) , ... , y(k)
// return : y(x)
template <typename mint> mint lagrange_interpolation(const vector<mint> &y, long long x, Binomial<mint> &C) {
int N = (int)y.size() - 1;
if(x <= N) return y[x];
mint ret = 0;
vector<mint> dp(N + 1, 1), pd(N + 1, 1);
mint a = x, one = 1;
for(int i = 0; i < N; i++) dp[i + 1] = dp[i] * a, a -= one;
for(int i = N; i > 0; i--) pd[i - 1] = pd[i] * a, a += one;
for(int i = 0; i <= N; i++) {
mint tmp = y[i] * dp[i] * pd[i] * C.finv(i) * C.finv(N - i);
ret += ((N - i) & 1) ? -tmp : tmp;
}
return ret;
}
#line 2 "library/fps/nth-term.hpp"
#line 5 "library/fps/nth-term.hpp"
template <typename mint> mint nth_term(long long n, const vector<mint> &s) {
using fps = FormalPowerSeries<mint>;
auto bm = BerlekampMassey<mint>(s);
return kitamasa(n, fps{begin(bm), end(bm)}, fps{begin(s), end(s)});
}
/**
* @brief 線形回帰数列の高速計算(Berlekamp-Massey/Bostan-Mori)
* @docs docs/fps/nth-term.md
*/
#line 12 "a.cpp"
using mint = LazyMontgomeryModInt<998244353>;
using fps = FormalPowerSeries<mint>;
using vmint = vector<mint>;
Binomial<mint> binomial;
mint inv(int i) { return binomial.inv(i); }
mint C(int r, int c) { return binomial.C(r, c); }
mint P(int r, int c) { return binomial.P(r, c); }
mint fact(int r) { return binomial.fac(r); }
mint ifact(int r) { return binomial.finv(r); }
} // namespace Modular998
using namespace Modular998;
int main() {
INT(n, x, y);
y = y - x;
vector<fps> dp(n + 5);
dp[0] = fps{0, 1};
rep(i, 1, n + 5) {
fps f = dp[i - 1] * x;
if(i > 1) { f += dp[i - 2] * y; }
dp[i] = (f.exp(n) << 1);
}
// vv(mint, f, 0);
vector<fps> f(n + 5);
swap(f, dp);
vector<fps> g(n + 5);
rep(i, 1, n + 5) g[i] = f[i] - f[i - 1];
g[0] = fps{0, 1};
// rep(i, n + 5) dump(i, f[i]);
// rep(i, n + 5) dump(i, g[i]);
mint ans;
// 中心が頂点
rep(d, 1, n + 5) {
auto a = (d >= 2 ? f[d - 2] : fps()) * x + (d >= 3 ? f[d - 3] : fps()) * y;
auto b = (d >= 1 ? g[d - 1] : fps()) * x + (d >= 2 ? g[d - 2] : fps()) * y;
a.resize(n), b.resize(n);
// dump(a, b);
// dump(d, ((a + b).exp(n)[n - 1] - a.exp(n)[n - 1] - (a.exp(n) * b)[n - 1]));
ans += ((a + b).exp(n)[n - 1] - a.exp(n)[n - 1] - (a.exp(n) * b)[n - 1]) * (d * 2);
}
// 中心が辺、d + d
rep(d, n + 5) {
auto a = g[d];
a = (a * a) * inv(2);
if(si(a) > n) {
ans += mint(x * (d * 2 + 1) + y * (d * 2 + 2)) * a[n];
// dump(d, a[n]);
}
}
// 中心が辺、d + (d + 1)
rep(d, n + 4) {
auto a = g[d];
auto b = g[d + 1];
a = a * b;
if(si(a) > n) {
ans += mint(y * (d * 2 + 3)) * a[n];
// dump(d, a[n]);
}
}
dump(ans);
OUT(ans * fact(n) / mint(n).pow(n - 2) / mint(x + y).pow(n - 1));
}
Details
Tip: Click on the bar to expand more detailed information
Test #1:
score: 100
Accepted
time: 0ms
memory: 3684kb
input:
2 1 3
output:
665496237
result:
ok 1 number(s): "665496237"
Test #2:
score: 0
Accepted
time: 1ms
memory: 5964kb
input:
3 2 3
output:
665496238
result:
ok 1 number(s): "665496238"
Test #3:
score: 0
Accepted
time: 2046ms
memory: 52880kb
input:
2000 1 2
output:
254870088
result:
ok 1 number(s): "254870088"
Test #4:
score: 0
Accepted
time: 2026ms
memory: 52740kb
input:
2000 1 3
output:
193693601
result:
ok 1 number(s): "193693601"
Test #5:
score: 0
Accepted
time: 2040ms
memory: 52680kb
input:
1999 188 211
output:
463395288
result:
ok 1 number(s): "463395288"
Test #6:
score: 0
Accepted
time: 2033ms
memory: 52620kb
input:
1990 470 818
output:
479264654
result:
ok 1 number(s): "479264654"
Test #7:
score: 0
Accepted
time: 508ms
memory: 17324kb
input:
1000 407 783
output:
20089106
result:
ok 1 number(s): "20089106"
Test #8:
score: 0
Accepted
time: 503ms
memory: 17128kb
input:
990 884 901
output:
94051884
result:
ok 1 number(s): "94051884"
Test #9:
score: 0
Accepted
time: 500ms
memory: 17164kb
input:
995 873 988
output:
209191626
result:
ok 1 number(s): "209191626"
Test #10:
score: 0
Accepted
time: 131ms
memory: 8448kb
input:
500 307 938
output:
603465152
result:
ok 1 number(s): "603465152"
Test #11:
score: 0
Accepted
time: 124ms
memory: 8736kb
input:
490 237 732
output:
402554558
result:
ok 1 number(s): "402554558"
Test #12:
score: 0
Accepted
time: 124ms
memory: 9000kb
input:
495 473 511
output:
833418554
result:
ok 1 number(s): "833418554"
Test #13:
score: 0
Accepted
time: 35ms
memory: 6800kb
input:
250 69 207
output:
786182422
result:
ok 1 number(s): "786182422"
Test #14:
score: 0
Accepted
time: 29ms
memory: 6512kb
input:
240 184 259
output:
473414786
result:
ok 1 number(s): "473414786"
Test #15:
score: 0
Accepted
time: 29ms
memory: 6564kb
input:
245 478 807
output:
312847415
result:
ok 1 number(s): "312847415"
Test #16:
score: 0
Accepted
time: 10ms
memory: 6148kb
input:
125 112 253
output:
31497383
result:
ok 1 number(s): "31497383"
Test #17:
score: 0
Accepted
time: 10ms
memory: 5980kb
input:
120 137 498
output:
923043504
result:
ok 1 number(s): "923043504"
Test #18:
score: 0
Accepted
time: 8ms
memory: 6156kb
input:
100 230 792
output:
203877027
result:
ok 1 number(s): "203877027"
Extra Test:
score: 0
Extra Test Passed