結果
| 問題 | No.3538 Not First Place |
| コンテスト | |
| ユーザー |
 jastaway
|
| 提出日時 | 2026-05-08 23:46:39 |
| 言語 | C++23 (gcc 15.2.0 + boost 1.89.0) |
| 結果 |
AC
|
| 実行時間 | 411 ms / 5,000 ms |
| コード長 | 55,463 bytes |
| 記録 | |
| コンパイル時間 | 7,133 ms |
| コンパイル使用メモリ | 438,476 KB |
| 実行使用メモリ | 67,328 KB |
| 最終ジャッジ日時 | 2026-05-08 23:46:55 |
| 合計ジャッジ時間 | 10,950 ms |
|
ジャッジサーバーID （参考情報） |
judge1_1 / judge2_1 |
（要ログイン）
| ファイルパターン | 結果 |
|---|---|
| sample | AC * 2 |
| other | AC * 26 |
ソースコード
#ifdef _DEBUG
#define _GLIBCXX_DEBUG
#endif
#include <bits/stdc++.h>
using namespace std;
#include <atcoder/all>
using namespace atcoder;
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>
#include <ext/pb_ds/tag_and_trait.hpp>
#include <ext/rope>
using namespace __gnu_pbds;
template<class T> using gset = tree<T,null_type,less<T>,rb_tree_tag,tree_order_statistics_node_update>;
template<class T, class U> using gmap = tree<T,U,less<T>,rb_tree_tag,tree_order_statistics_node_update>;
template<class T, class U> using ghmap = gp_hash_table<T, U>;
// #include <boost/rational.hpp>
// using namespace boost;
// using rat = rational<long long int>;
using mint = modint998244353;
// using mint = modint1000000007;
// using mint = modint;
using ll = long long;
using lll = __int128_t;
using ld = long double;
using ull = uint64_t;
using pll = pair<ll, ll>;
using vll = vector<ll>;
using vvll = vector<vll>;
using vvvll = vector<vvll>;
using vpll = vector<pll>;
using vvpll = vector<vpll>;
using vm = vector<mint>;
using vvm = vector<vm>;
using vvvm = vector<vvm>;
using vstr = vector<string>;
#define v(T) vector<T>
#define vv(T) vector<vector<T>>
#define vvv(T) vector<vector<vector<T>>>
#define vvvv(T) vector<vector<vector<vector<T>>>>
template<int M> istream &operator>>(istream &is, static_modint<M> &a){ll tmp; is >> tmp; a = tmp; return is;}
template<int M> ostream &operator<<(ostream &os, const static_modint<M> &a){ os << a.val(); return os; }
template<int M> istream &operator>>(istream &is, dynamic_modint<M> &a){ll tmp; is >> tmp; a = tmp; return is;}
template<int M> ostream &operator<<(ostream &os, const dynamic_modint<M> &a){ os << a.val(); return os; }
string to_string(const __int128_t &a) { if (a == 0) return "0"; string s = ""; __int128_t num = a; bool is_negative = false; if (num < 0) { is_negative = true; num = -num; } while (num > 0) { s += '0' + (num % 10); num /= 10; } if (is_negative) s += '-'; reverse(s.begin(), s.end()); return s; }
istream &operator>>(istream &is, __int128_t &a){ string s; is >> s; a = 0; for(char c : s) { if(isdigit(c)) {a = a*10 + (c - '0'); } } if(s[0]=='-'){ a *= -1; } return is; }
ostream &operator<<(ostream &os, const __int128_t &a){ os << to_string(a); return os; }
template<class T, class U> istream &operator>>(istream &is, pair<T, U> &p) { is >> p.first >> p.second; return is; }
template<class T, class U> ostream &operator<<(ostream &os, const pair<T, U> &p) { if(&os == &std::cerr) { os << "(" << p.first << ", " << p.second << ")"; } else { os << p.first << " " << p.second; } return os; }
template<class T> istream &operator>>(istream &is, vector<T> &vec){ for(T &e : vec){is >> e;} return is; }
template<class T> ostream &operator<<(ostream &os, const vector<T> &vec) { for(int i = 0; i < (int)vec.size(); i++) { os << vec[i] << (i + 1 != (int)vec.size() ? " " : ""); } return os; }
template<class T> istream &operator>>(istream &is, deque<T> &vec){ for(T &e : vec){is >> e;} return is; }
template<class T> ostream &operator<<(ostream &os, const deque<T> &vec) { for(int i = 0; i < (int)vec.size(); i++) { os << vec[i] << (i + 1 != (int)vec.size() ? " " : ""); } return os; }
template<class T> ostream &operator<<(ostream &os, const set<T> &vec) { for(auto itr = vec.begin(); itr != vec.end(); itr++) { os << (itr != vec.begin() ? " " : "") << (*itr); } return os; }
template<class T> ostream &operator<<(ostream &os, const unordered_set<T> &vec) { for(auto itr = vec.begin(); itr != vec.end(); itr++) { os << (itr != vec.begin() ? " " : "") << (*itr); } return os; }
template<class T, class U> ostream &operator<<(ostream &os, const map<T, U> &vec) { for(auto itr = vec.begin(); itr != vec.end(); itr++) { os << (itr != vec.begin() ? " " : "") << (*itr); } return os; }
template<class T, class U> ostream &operator<<(ostream &os, const unordered_map<T, U> &vec) { for(auto itr = vec.begin(); itr != vec.end(); itr++) { os << (itr != vec.begin() ? " " : "") << (*itr); } return os; }
template <class Tuple, size_t... Is> void read_tuple(istream& is, Tuple& t, index_sequence<Is...>) { ((is >> std::get<Is>(t)), ...); }
template <class... Ts> istream& operator>>(istream& is, tuple<Ts...>& t) { read_tuple(is, t, index_sequence_for<Ts...>{}); return is; }
template <class Tuple, size_t... Is> void print_tuple(ostream& os, const Tuple& t, index_sequence<Is...>) { ((os << (Is == 0 ? "" : (&os == &std::cerr ? ", " : " ")) << std::get<Is>(t)), ...); }
template <class... Ts> ostream& operator<<(ostream& os, const tuple<Ts...>& t) { if(&os == &std::cerr) { os << "("; } print_tuple(os, t, index_sequence_for<Ts...>{}); if(&os == &std::cerr) { os << ")"; } return os; }
#ifdef _DEBUG
template <typename... Args>
void debug_out(Args... args) { ((std::cerr << " " << args), ...); std::cerr << "\n"; }
#define debug(...) do { \
std::cerr << "[" << #__VA_ARGS__ << "]:"; \
debug_out(__VA_ARGS__); \
} while(0)
#else
#define debug(...) (void)0
#endif
template <class... T> constexpr auto min (T... a) { return min(initializer_list<common_type_t<T...>>{a...}); }
template <class... T> constexpr auto max (T... a) { return max(initializer_list<common_type_t<T...>>{a...}); }
template<class T> using pqg = priority_queue<T, vector<T>, greater<T>>;
template<class T> T opmin(T x, T y) { return min(x, y); }
template<class T> T einf() { return numeric_limits<T>::max()/2; }
template<class T> T opmax(T x, T y) { return max(x, y); }
template<class T> T eminf() { return numeric_limits<T>::min()/2; }
template<class T> T opsum(T x, T y) { return x + y; }
template<class T> T ezero() { return (T)0; }
template<class T> using minseg = segtree<T, opmin<T>, einf<T>>;
template<class T> using maxseg = segtree<T, opmax<T>, eminf<T>>;
template<class T> using sumseg = segtree<T, opsum<T>, ezero<T>>;
// template<class T> struct v : vector<T> { using vector<T> :: vector; };
// template<class T> struct vv : vector<v<T>> { using vector<v<T>> :: vector; };
// template<class T> struct vvv : vector<vv<T>> { using vector<vv<T>> :: vector; };
template<class T, class U> inline bool chmin(T& a, U b) {if(a > b){a = b; return true;} else {return false;}};
template<class T, class U> inline bool chmax(T& a, U b) {if(a < b){a = b; return true;} else {return false;}};
template<class T> T dist_sq(const pair<T, T> &x, const pair<T, T> &y){return (x.first-y.first)*(x.first-y.first)+(x.second-y.second)*(x.second-y.second);}
template<class T> T cross_product(const pair<T, T> &x, const pair<T, T> &y){return x.first * y.second - x.second * y.first;}
template <typename T> struct Reverser {
T& t;
auto begin() const { return std::rbegin(t); }
auto end() const { return std::rend(t); }
};
template <typename T> Reverser<T> reversed(T& t) { return {t}; }
template <typename T> Reverser<const T> reversed(const T& t) { return {t}; }
// #define rep(i,n) for(ll i = 0; i < (ll)(n); i++)
#define _GET5(_1,_2,_3,_4,NAME,...) NAME
#define _rep1(i, n) for(long long i=0;i<(long long)(n);i++)
#define _rep2(i,k,n) for(long long i=(long long)(k);i<(long long)(n);i++)
#define _rep3(i,k,n,s) for(long long i=(long long)(k);i<(long long)(n);i+=(s))
#define rep(...) _GET5(__VA_ARGS__, _rep3, _rep2, _rep1)(__VA_ARGS__)
#define _rrep1(i, n) for(long long i=(long long)(n) - 1;i>=0;i--)
#define _rrep2(i,k,n) for(long long i=(long long)(n)-1;i>=(long long)(k);i--)
#define _rrep3(i,k,n,s) for(long long i=(long long)(n);i>=(long long)(k);i-=(s))
#define rrep(...) _GET5(__VA_ARGS__, _rrep3, _rrep2, _rrep1)(__VA_ARGS__)
#define repr(i,n) for(ll i = (ll)(n) - 1; i >= 0; i--)
#define REP(i, l, r) for(ll i = (ll)l; i <= (ll)(r); i++)
#define REPR(i, l, r) for(ll i = (ll)r; i >= (ll)(l); i--)
const ll inf = (1 << 30);
const ll INF = (1LL << 60);
const vector<pair<ll, ll>> DIJ = {{1, 0}, {0, -1}, {-1, 0}, {0, 1}};
void out(){cout<<'\n';}
template<class T, class... Ts> void out(const T& a, const Ts&... b){ cout<<a; (cout<<... << (cout << ' ', b)); cout << '\n';}
void outf(){cout<<endl;}
template<class T, class... Ts> void outf(const T& a, const Ts&... b){ cout<<a; (cout<<... << (cout << ' ', b)); cout << endl;}
template<class T, class U> void outp(pair<T, U> a){ out((a).first, (a).second); }
template<class T, class U> void outpf(pair<T, U> a){ outf((a).first, (a).second); }
template<class T> void outv(T a){rep(i, (a).size()){ cout << (a)[i] << " "; } cout << endl;}
template<class T> void outvL(T a){rep(i, (a).size()){out((a)[i]);} cout << flush; }
// template<class T> void outvv(T a){rep(i, a.size()){ rep(j, a.at(i).size()){cout << a.at(i).at(j) << " "; } cout << endl; }}
// template<class T> void outvp(T a){rep(i, a.size()){ out2(a.at(i).first, a.at(i).second); }}
void setpre(int a){cout << fixed << setprecision(a);}
#define outN cout << "No\n"
#define outY cout << "Yes\n"
void outYN(bool flag) { cout << (flag ? "Yes\n" : "No\n"); }
#define dame(...) {outf(__VA_ARGS__);return 0;}
template<class... T> void read(T&... a){(cin >> ... >> a);}
#define readll(...) ll __VA_ARGS__; read(__VA_ARGS__)
#define readvll(a, n) vector<ll> a(n); read(a)
#define readvvll(a, n, m) vector<vector<ll>> a(n, vector<ll>(m)); read(a)
#define readvstr(a, n) vector<string> a(n); read(a)
#define readvt(type, a, n) vector<type> a(n); read(a)
#define readvll2(a, b, n) vector<ll> a(n), b(n); for(int lopi = 0; lopi < (int)(n); lopi++) cin >> (a)[lopi] >> (b)[lopi]
#define readvll3(a, b, c, n) vector<ll> a(n), b(n), c(n); for(int lopi = 0; lopi < (int)(n); lopi++) cin >> (a)[lopi] >> (b)[lopi] >> (c)[lopi]
#define readstr(...) string __VA_ARGS__; read(__VA_ARGS__)
#define readundirG(G, N, M) G = vvll(N); rep(lopi, M) {ll a, b; cin >> a >> b; G[a-1].push_back(b-1); G[b-1].push_back(a-1);}
#define readdirG(G, N, M) G = vvll(N); rep(lopi, M) {ll a, b; cin >> a >> b; G[a-1].push_back(b-1);}
#define readundirwghG(G, N, M) G = vv(pll)(N); rep(lopi, M) {ll a, b, c; cin >> a >> b >> c; G[a-1].emplace_back(b-1,c); G[b-1].emplace_back(a-1, c);}
#define readdirwghG(G, N, M) G = vv(pll)(N); rep(lopi, M) {ll a, b, c; cin >> a >> b >> c; G[a-1].emplace_back(b-1, c);}
#define All(a) (a).begin(), (a).end()
#define all(a) begin(a), end(a)
template<class T> inline void sortr(T& a){ sort((a).rbegin(), (a).rend()); }
template<class T> inline vector<int> argsort(T V, bool rev = false){vector<int> res(V.size()); iota(res.begin(), res.end(), 0); sort(res.begin(), res.end(), [&](int x, int y){if(!rev){return V[x] < V[y];}else{return V[x] > V[y];}}); return res;}
template<class T, class U> inline void sort_by_idx(T& V, vector<U>& I){assert(V.size() == I.size()); T tmpv = V; for(int loopi = 0; loopi < (int)I.size(); loopi++){V[loopi] = tmpv[I.at(loopi)];}}
template<class T, class U> inline void sortp(vector<T>& v1, vector<U>& v2, bool rev1 = false, int rev2 = false){assert(v1.size() == v2.size()); vector<int> I(v1.size()); iota(I.begin(), I.end(), 0); sort(I.begin(), I.end(), [&](const int x, const int y){if(v1[x] != v1[y]){return (bool)(rev1 ^ (v1[x] < v1[y]));}else{if(v2[x]==v2[y]){return false;} return (bool)(rev2 ^ (v2[x] < v2[y]));}}); sort_by_idx(v1, I); sort_by_idx(v2, I);}
template<class T> T POW(T x, ll n) {T ret = 1; while(n > 0){if(n & 1) ret *= x; x *= x; n >>= 1;} return ret;}
ll powll(ll x, ll n){ll ret = 1; while(n > 0){if(n & 1) ret *= x; x *= x; n >>= 1;} return ret;}
ll modpow(ll a,ll b,ll modd){ll res=1;a%=modd;while(b){if(b&1)res=res*a%modd;a=a*a%modd;b>>=1;}return res;}
ll sqrtll(ll a){assert(a >= 0); ll r = (ll)sqrtl((ld)a)-1; while(r < 0 || (r+1)*(r+1) <= a) r++; return r;}
ll cbrtll(ll a){assert(a >= 0); ll r = (ll)cbrtl((ld)a)-1; while(r < 0 || (r+1)*(r+1)*(r+1) <= a) r++; return r;}
ll modinv(ll a, ll b){assert(a); if(a == 1) return 1; if(b == 1) return 0; ll ret = (1ll-b*modinv(b%a, a))/a; ret %= b; if(ret < 0) ret += b; return ret;}
inline ll divceil(ll x, ll y) { if(x >= 0) {return(x / y + (ll)(x % y != 0)); } else { return -((-x) / y); } }
inline ll divfloor(ll x, ll y) { if(x >= 0) { return x/y; } else { return -((-x)/y + (ll)((-x) % y != 0)); } }
inline bool inLR(ll x, ll L, ll R){ return (L <= x && x < R); }
inline bool inRect(ll pos_x, ll pos_y, ll rect_H, ll rect_W, ll rect_h = 0, ll rect_w = 0){ return (rect_h <= pos_x && pos_x < rect_H && rect_w <= pos_y && pos_y < rect_W); }
template<class T> vector<T> &operator++(vector<T>& v){for(T& x : v) x++; return v;}
template<class T> vector<T> &operator--(vector<T>& v){for(T& x : v) x--; return v;}
template<class T> vector<T> operator++(vector<T>& v, signed){auto res = v; for(T& x : v) x++; return res;}
template<class T> vector<T> operator--(vector<T>& v, signed){auto res = v; for(T& x : v) x--; return res;}
template<class T> vector<T> operator+=(vector<T>& v, const vector<T>& w){if(v.size() < w.size()) v.resize(w.size()); for(int i = 0; i < (int)w.size(); i++) v[i] += w[i]; return v;}
template<class T> vector<T> operator-=(vector<T>& v, const vector<T>& w){if(v.size() < w.size()) v.resize(w.size()); for(int i = 0; i < (int)w.size(); i++) v[i] -= w[i]; return v;}
template<class T> vector<T> operator*=(vector<T>& v, const vector<T>& w){if(v.size() < w.size()) v.resize(w.size()); for(int i = 0; i < (int)w.size(); i++) v[i] *= w[i]; return v;}
template<class T> vector<T> operator/=(vector<T>& v, const vector<T>& w){if(v.size() < w.size()) v.resize(w.size()); for(int i = 0; i < (int)w.size(); i++) v[i] /= w[i]; return v;}
template<class T> vector<T> operator+(vector<T> v, const vector<T>& w){return (v += w);}
template<class T> vector<T> operator-(vector<T> v, const vector<T>& w){return (v -= w);}
template<class T> vector<T> operator*(vector<T> v, const vector<T>& w){return (v *= w);}
template<class T> vector<T> operator/(vector<T> v, const vector<T>& w){return (v /= w);}
template<class T> vector<T> operator*=(vector<T>& v, T w){for(T& x : v) x*=w; return v;}
template<class T> vector<T> operator/=(vector<T>& v, T w){for(T& x : v) x/=w; return v;}
template<class T> vector<T> operator*(vector<T> v, T x){return (v *= x);}
template<class T> vector<T> operator/(vector<T> v, T x){return (v /= x);}
template<class S, class T> pair<S,T> operator+=(pair<S,T>& p, const pair<S,T>& q){p.fi+=q.fi, p.se+=q.se; return p;}
template<class S, class T> pair<S,T> operator+(pair<S,T> p, const pair<S,T>& q){return (p += q);}
template<class S, class T> pair<S,T> operator-=(pair<S,T>& p, const pair<S,T>& q){p.fi-=q.fi, p.se-=q.se; return p;}
template<class S, class T> pair<S,T> operator-(pair<S,T> p, const pair<S,T>& q){return (p -= q);}
template<class T> size_t HashCombine(const size_t seed,const T &v){ return seed^(std::hash<T>()(v)+0x9e3779b9+(seed<<6)+(seed>>2)); }
template<class T,class S> struct std::hash<std::pair<T,S>>{ size_t operator()(const std::pair<T,S> &keyval) const noexcept { return HashCombine(std::hash<T>()(keyval.first), keyval.second); } };
template<class T> struct std::hash<std::vector<T>>{ size_t operator()(const std::vector<T> &keyval) const noexcept { size_t s=0; for (auto&& v: keyval) s=HashCombine(s,v); return s; } };
template<int N> struct HashTupleCore{ template<class Tuple> size_t operator()(const Tuple &keyval) const noexcept{ size_t s=HashTupleCore<N-1>()(keyval); return HashCombine(s,std::get<N-1>(keyval)); } };
template <> struct HashTupleCore<0>{ template<class Tuple> size_t operator()(const Tuple &keyval) const noexcept{ return 0; } };
template<class... Args> struct std::hash<std::tuple<Args...>>{ size_t operator()(const tuple<Args...> &keyval) const noexcept { return HashTupleCore<tuple_size<tuple<Args...>>::value>()(keyval); } };
template <class T>
class modint_factorials
{
private:
vector<T> invs;
vector<T> fact;
vector<T> inv_fact;
int max_calculated = -1;
modint_factorials() = default;
void calculate_until(int n)
{
assert(n >= 0);
long long mod = T::mod();
if (n < 1)
n = 1;
invs.resize(n + 1);
fact.resize(n + 1);
inv_fact.resize(n + 1);
if (max_calculated == -1)
{
invs[0] = 0;
invs[1] = 1;
fact[0] = 1;
fact[1] = 1;
inv_fact[0] = 1;
inv_fact[1] = 1;
max_calculated = 1;
}
if (n <= max_calculated)
return;
for (int i = max_calculated + 1; i <= n; i++)
{
invs[i] = invs[mod % i] * -(mod / i);
fact[i] = fact[i - 1] * i;
inv_fact[i] = inv_fact[i - 1] * invs[i];
}
max_calculated = n;
}
public:
static modint_factorials &instance()
{
static modint_factorials inst;
return inst;
}
static void precompute(int n)
{
assert(n >= 0);
instance().calculate_until(n);
}
static T mint_inv(int x)
{
assert(x != 0);
if (x < 0)
return -mint_inv(-x);
modint_factorials &inst = instance();
if (x > inst.max_calculated)
inst.calculate_until(x);
return inst.invs[x];
}
static T mint_factorial(int x)
{
assert(x >= 0);
modint_factorials &inst = instance();
if (x > inst.max_calculated)
inst.calculate_until(x);
return inst.fact[x];
}
static T mint_inv_factorial(int x)
{
assert(x >= 0);
modint_factorials &inst = instance();
if (x > inst.max_calculated)
inst.calculate_until(x);
return inst.inv_fact[x];
}
static T nPk(int n, int k)
{
if (k < 0)
return 0;
if (n < 0)
return ((k & 1) ? -nPk(-n + k - 1, k) : nPk(-n + k - 1, k));
if (n >= 0 && k > n)
return 0;
modint_factorials &inst = instance();
if (n > inst.max_calculated)
inst.calculate_until(n);
return inst.fact[n] * inst.inv_fact[n - k];
}
static T nCk(int n, int k)
{
if (k < 0)
return 0;
if (n < 0)
return ((k & 1) ? -nCk(-n + k - 1, k) : nCk(-n + k - 1, k));
if (n >= 0 && k > n)
return 0;
modint_factorials &inst = instance();
if (n > inst.max_calculated)
inst.calculate_until(n);
return inst.fact[n] * inst.inv_fact[k] * inst.inv_fact[n - k];
;
}
};
using mint_factorials = modint_factorials<mint>;
struct FPS : vector<mint>
{
using vector<mint>::vector;
using vm = vector<mint>;
using vvm = vector<vector<mint>>;
using pm = pair<mint, mint>;
using vpm = vector<pair<mint, mint>>;
using ll = long long;
using pi = pair<ll, ll>;
using vi = vector<ll>;
using vvi = vector<vector<ll>>;
FPS(const vm &v)
{
FPS ret;
ret.resize(v.size());
rep(i, v.size()) ret[i] = v[i];
*this = ret;
}
FPS(const vi &v)
{
FPS ret;
ret.resize(v.size());
rep(i, v.size()) ret[i] = v[i];
*this = ret;
}
FPS operator+=(const FPS &other)
{
if (other.size() > size())
resize(other.size());
rep(i, other.size()) at(i) += other[i];
return *this;
}
FPS operator-=(const FPS &other)
{
if (other.size() > size())
resize(other.size());
rep(i, other.size()) at(i) -= other[i];
return *this;
}
FPS operator+=(const mint &m)
{
if (empty())
resize(1);
at(0) += m;
return *this;
}
FPS operator-=(const mint &m)
{
if (empty())
resize(1);
at(0) -= m;
return *this;
}
FPS operator*=(const mint &m)
{
rep(i, size()) at(i) *= m;
return *this;
}
FPS operator/=(const mint &m)
{
mint mi = m.inv();
rep(i, size()) at(i) *= mi;
return *this;
}
FPS operator+(const FPS &other) const { return FPS(*this) += other; }
FPS operator-(const FPS &other) const { return FPS(*this) -= other; }
FPS operator+(const mint &m) const { return FPS(*this) += m; }
FPS operator-(const mint &m) const { return FPS(*this) -= m; }
FPS operator*(const mint &m) const { return FPS(*this) *= m; }
FPS operator/(const mint &m) const
{
mint mi = m.inv();
return FPS(*this) *= mi;
}
FPS operator-() const
{
FPS ret(size());
rep(i, size()) ret[i] = -at(i);
return ret;
}
FPS operator*(const FPS &other) const
{
return FPS(convolution(*this, other));
}
FPS operator*=(const FPS &other)
{
*this = FPS(convolution(*this, other));
return *this;
}
void shrink()
{
while (size() > 0 && back() == 0)
pop_back();
}
void reversed()
{
reverse(begin(), end());
}
FPS rev() const
{
FPS ret(*this);
reverse(ret.begin(), ret.end());
return ret;
}
FPS pre(ll k) const
{
return FPS(begin(), begin() + min((ll)size(), k));
}
FPS dot(FPS g) const
{
FPS ret(*this);
if (ret.size() < g.size())
ret.resize(g.size());
rep(i, g.size()) ret[i] *= g[i];
return ret;
}
FPS inv(ll N = -1) const
{
if (N == -1)
N = size();
FPS g;
g += at(0).inv();
ll k = 1;
while (k < N)
{
k *= 2;
g = (-pre(k) * g + 2) * g;
g.resize(k);
}
g.resize(N);
return g;
}
FPS operator/(const FPS &other) const
{
FPS ret = (*this) * other.inv(size());
ret.resize(size());
return ret;
}
FPS operator/=(const FPS &other)
{
*this = (*this) / other;
return *this;
}
FPS operator>>(ll k) const
{
if ((ll)size() < k)
return {};
FPS ret(*this);
ret.erase(ret.begin(), ret.begin() + k);
return ret;
}
FPS operator<<(ll k) const
{
FPS ret(*this);
ret.insert(ret.begin(), k, mint(0));
return ret;
}
FPS diff() const
{
FPS ret = (*this) >> 1;
rep(i, ret.size()) ret[i] *= i + 1;
return ret;
}
FPS integral() const
{
FPS ret = (*this) << 1;
rep(i, ret.size() - 1) ret[i + 1] *= mint_factorials::mint_inv(i + 1);
return ret;
}
FPS log(ll N = -1) const
{
FPS ret(*this);
if (N > -1)
ret.resize(N);
return (ret.diff() / ret).integral();
}
FPS exp(ll N = -1) const
{
if (N == -1)
N = size();
FPS g;
g += 1;
ll k = 1;
while (k < N)
{
k *= 2;
g = (pre(k) - g.log(k) + 1) * g;
g.resize(k);
}
g.resize(N);
return g;
}
FPS pow(ll k, ll N = -1) const
{
if (N == -1)
N = size();
if (k == 0)
{
FPS ret(N);
if (N > 0)
ret += 1;
return ret;
}
FPS ret(*(this));
ll n = ret.size();
ll c = 0;
while (c < n && ret[c] == 0)
c++;
if (c == n || c >= (N + k - 1) / k)
return FPS(N);
ret = ret >> c;
mint a = ret[0];
ret /= a;
return ((ret.log(N - c * k) * k).exp(N - c * k) << (c * k)) * a.pow(k);
}
static mint modsqrt(mint a)
{
static constexpr int MOD = mint::mod();
if (MOD == 2)
return a;
static constexpr int M = __builtin_ctzll(MOD - 1);
static constexpr long long Q = (MOD - 1) >> M;
static constexpr long long Q2 = (Q + 1) >> 1;
static constexpr long long R = (MOD - 1) / 2;
static constexpr long long G = internal::pow_mod_constexpr(internal::primitive_root<MOD>, Q, MOD);
if (a == 0)
return 0;
if (a.pow(R) == -1)
return -1;
int k = 0;
mint b = a.pow(Q);
for (int i = 0; i < M; i++)
{
if (b.pow(1ll << (M - i - 1)) == -1)
{
b *= mint(G).pow(1ll << i);
k |= 1ll << i;
}
}
return a.pow(Q2) * mint(G).pow(k >> 1);
}
FPS _sqrt(ll N = -1) const
{
if (N == -1)
N = size();
FPS g;
g += 1;
ll k = 1;
while (k < N)
{
k *= 2;
g = (g + pre(k) * g.inv(k)) / 2;
g.resize(k);
}
g.resize(N);
return g;
}
FPS sqrt(ll N = -1) const
{
if (N == -1)
N = size();
FPS ret(*(this));
ll n = ret.size();
ll c = 0;
while (c < n && ret[c] == 0)
c++;
if (c == n)
return FPS(N);
if (c & 1)
return {};
ret = ret >> c;
mint a = ret[0];
mint b = modsqrt(a);
if (b == -1)
return {};
ret /= a;
return (ret._sqrt(N - c / 2) << (c / 2)) * b;
}
pair<FPS, FPS> div(FPS g) const
{
FPS f(*this);
f.shrink();
g.shrink();
ll n = f.size();
ll m = g.size();
if (n < m)
return {{}, f};
FPS q = f.rev();
g.reversed();
q /= g;
q.resize(n - m + 1);
q.reversed();
g.reversed();
FPS r = f - g * q;
r.shrink();
return {q, r};
}
mint eval(mint x) const
{
mint ans = 0;
rrep(i, size())
{
ans *= x;
ans += at(i);
}
return ans;
}
vm mult_eval(vm p) const
{
ll n = p.size();
if (n == 0)
return {};
ll N = 1ll << (64 - __builtin_clzll(n - 1));
vector<FPS> g(2 * N);
vector<FPS> r(N);
vm ans(n);
rep(i, n) g[i + N] = {-p[i], 1};
rep(i, n, N) g[i + N] = {1};
rrep(i, 1, N) g[i] = g[2 * i] * g[2 * i + 1];
r[0] = (*this);
rep(i, 1, N) r[i] = r[i / 2].div(g[i]).second;
rep(i, n) ans[i] = r[(i + N) / 2].eval(p[i]);
return ans;
}
vm to_newton(vm p) const
{
assert(p.size() == size());
ll n = size();
if (n == 0)
return {};
ll N = 1ll << (64 - __builtin_clzll(n - 1));
vector<FPS> g(2 * N);
vector<FPS> r(2 * N);
vm ans(n);
rep(i, n) g[i + N] = {-p[i], 1};
rep(i, n, N) g[i + N] = {1};
rrep(i, 1, N) g[i] = g[2 * i] * g[2 * i + 1];
r[1] = (*this);
rep(i, 1, N)
{
auto [p, q] = r[i].div(g[2 * i]);
r[2 * i] = q, r[2 * i + 1] = p;
}
rep(i, n) ans[i] = r[i + N].empty() ? 0 : r[i + N][0];
return ans;
}
static FPS prod(vector<FPS> F)
{
ll n = F.size();
if (n == 0)
return {1};
priority_queue<pi, vector<pi>, greater<pi>> q;
rep(i, n) q.push({F[i].size(), i});
while (q.size() > 1)
{
ll i = q.top().second;
q.pop();
ll j = q.top().second;
q.pop();
F[i] *= F[j];
q.push({F[i].size(), i});
}
return F[q.top().second];
}
static FPS quotient_sum(vector<pair<FPS, FPS>> QR, ll d = -1)
{
FPS ans;
ll n = QR.size();
rep(i, n)
{
if (QR[i].first.size() >= QR[i].second.size())
{
pair<FPS, FPS> p = QR[i].first.div(QR[i].second);
ans += p.first;
QR[i].first = p.second;
}
}
ll k = 1;
while (k < n)
{
for (ll i = 0; i + k < n; i += (k << 1))
{
QR[i] = {
QR[i].second * QR[i + k].first + QR[i].first * QR[i + k].second,
QR[i].second * QR[i + k].second};
}
k = k << 1;
}
if (d == -1)
d = max(ans.size(), QR[0].second.size());
ans += QR[0].first * QR[0].second.inv(d);
ans.resize(d);
return ans;
}
vm geo_eval(mint a, mint r, ll m)
{
if (m == 0)
return {};
ll n = size();
if (n == 0)
return {};
if (r == 0)
{
vm ret(m, at(0));
ret[0] = eval(a);
return ret;
}
FPS x(n + m, 1), y(n, 1);
mint e = 1;
mint ri = r.inv();
rep(i, n + m - 1)
{
x[i + 1] = e * x[i];
e *= r;
}
e = 1;
rrep(i, n - 1)
{
y[i] = a * e * y[i + 1];
e *= ri;
}
rep(i, n) y[i] *= at(n - i - 1);
FPS z = x * y;
vm ret(m);
e = 1;
mint f = 1;
rep(i, m)
{
ret[i] = f * z[n + i - 1];
f *= e;
e *= ri;
}
return ret;
}
static FPS interpolation(const vm &x, const vm &y)
{
ll n = x.size();
if (n == 0)
return {};
vector<FPS> g0(n);
rep(i, n) g0[i] = {-x[i], 1};
FPS G = prod(g0);
vm g = G.diff().mult_eval(x);
vector<pair<FPS, FPS>> qr(n);
ll flag = -1;
rep(i, n) qr[i].first = {y[i]};
rep(i, n)
{
if (x[i] == 0)
{
qr[i] = {{0}, {1}};
flag = i;
}
else
qr[i].second = {-g[i] * x[i], g[i]};
}
FPS ret = quotient_sum(qr, n) * G;
if (flag >= 0)
ret += (G >> 1) * (y[flag] / g[flag]);
ret.resize(n);
return ret;
}
static FPS geo_interpolation(const vm &y, mint a, mint r)
{
ll n = y.size();
if (!n)
return {};
if (n == 1)
return {y[0]};
FPS g(n + 1);
g[n] = 1;
mint rn = r.pow(n);
mint ri = r.inv();
mint ai = a.inv();
mint e = 1;
if (rn == 1)
g[0] = -a.pow(n);
else
{
rrep(i, n)
{
g[i] = g[i + 1] * a * (e - rn) / (e * r - 1);
e *= r;
}
}
vm s = g.diff().geo_eval(a, r, n);
g.pop_back();
FPS c(n);
rep(i, n) c[i] = y[i] / s[i];
FPS f(c.geo_eval(ri, ri, n + 1));
e = -ai;
rep(i, n + 1)
{
f[i] *= e;
e *= ai;
}
return (g * f).pre(n);
}
FPS Taylor_shift(mint a) const
{
if (a == 0)
return *this;
FPS ret(*this);
ll n = size();
rep(i, n) ret[i] *= mint_factorials::mint_factorial(i);
FPS g(n, 1);
rep(i, n - 1) g[i + 1] = g[i] * a;
rep(i, n) g[i] *= mint_factorials::mint_inv_factorial(i);
ret = ret.rev() * g;
ret.resize(n);
ret.reversed();
rep(i, n) ret[i] *= mint_factorials::mint_inv_factorial(i);
return ret;
}
// return (cx+d)^n f((ax+b)/(cx+d))
FPS Mobius_Transform(mint a, mint b, mint c, mint d) const
{
mint tmp;
if (empty())
return {};
int n = size() - 1;
if (c == 0)
{
if (d == 0)
{
FPS res(n + 1);
tmp = back();
rep(i, n + 1)
{
res[i] = mint_factorials::nCk(n, i) * tmp;
tmp *= a;
}
tmp = 1;
rrep(i, n + 1)
{
res[i] *= tmp;
tmp *= b;
}
return res;
}
FPS res(*this);
tmp = 1;
rrep(i, n + 1)
{
res[i] *= tmp;
tmp *= d;
}
rep(i, n + 1) res[i] *= mint_factorials::mint_factorial(i);
FPS g(n + 1, 1);
rep(i, n) g[i + 1] = g[i] * b * mint_factorials::mint_inv(i + 1);
reverse(res.begin(), res.end());
res = convolution(res, g);
res.resize(n + 1);
reverse(res.begin(), res.end());
rep(i, n + 1) res[i] *= mint_factorials::mint_inv_factorial(i);
tmp = 1;
rep(i, n + 1)
{
res[i] *= tmp;
tmp *= a;
}
return res;
}
a /= c;
FPS res(*this);
rep(i, n + 1) res[i] *= mint_factorials::mint_factorial(i);
FPS g(n + 1, 1);
rep(i, n) g[i + 1] = g[i] * a * mint_factorials::mint_inv(i + 1);
reverse(res.begin(), res.end());
res = convolution(res, g);
res.resize(n + 1);
b -= a * d;
tmp = 1;
rrep(i, n + 1)
{
res[i] *= tmp * mint_factorials::mint_factorial(i) * mint_factorials::mint_inv_factorial(n - i);
tmp *= b;
}
rep(i, n) g[i + 1] = g[i] * d * mint_factorials::mint_inv(i + 1);
reverse(res.begin(), res.end());
res = convolution(res, g);
res.resize(n + 1);
reverse(res.begin(), res.end());
rep(i, n + 1) res[i] *= mint_factorials::mint_inv_factorial(i);
tmp = 1;
rep(i, n + 1)
{
res[i] *= tmp;
tmp *= c;
}
return res;
}
// FPS composition(FPS G, ll deg = -1) const {
// if(deg == -1) deg = size();
// ll n = size();
// if((ll)G.size() > deg) G.resize(deg);
// if(n == 0) return FPS(deg);
// ll k = max((ll)sqrtl(n), 1ll);
// ll l = (n+k-1)/k;
// vector<FPS> Pi(l);
// FPS Qb = {1};
// rep(i,k){
// rep(j,l){
// if(i + j*k < n){
// Pi[j] += Qb * at(i + j*k);
// }
// }
// Qb *= G;
// if((ll)Qb.size() > deg) Qb.resize(deg);
// }
// FPS ret;
// FPS Qg = {1};
// rep(i,l){
// ret += (Pi[i] * Qg).pre(deg);
// Qg *= Qb;
// if((ll)Qg.size() > deg) Qg.resize(deg);
// }
// return ret;
// }
static vm shift_sampling(vm f, mint c, ll m = -1)
{
ll n = f.size() - 1;
if (m == -1)
m = n;
if (m == 1)
{
if (c.val() <= n)
return {f[c.val()]};
mint t = 1;
mint s1 = 1;
mint s2 = 1;
rep(i, n + 1)
{
f[i] *= s1 / t;
f[n - i] *= s2 / t;
t *= i + 1;
s1 *= c - i;
s2 *= c - n + i;
if ((i + n) & 1)
f[i] *= -1;
}
mint ans = 0;
rep(i, n + 1) ans += f[i];
return {ans};
}
n++;
if (f.size() == 0)
return vm(0, m);
FPS F(n);
rep(i, n) F[i] = f[i] * mint_factorials::mint_inv_factorial(i);
FPS G(n);
rep(i, n) G[i] = mint_factorials::mint_inv_factorial(i) * (i & 1 ? -1 : 1);
F *= G;
F.resize(n);
rep(i, n) F[i] *= mint_factorials::mint_factorial(i);
FPS H(n, 1);
rep(i, n - 1) H[i + 1] = H[i] * (c - i);
rep(i, n) H[i] *= mint_factorials::mint_inv_factorial(i);
F = ((F * H.rev()) >> (n - 1));
rep(i, n) F[i] *= mint_factorials::mint_inv_factorial(i);
FPS I(m);
rep(i, m) I[i] = mint_factorials::mint_inv_factorial(i);
F = (F * I).pre(m);
rep(i, m) F[i] *= mint_factorials::mint_factorial(i);
return F;
}
FPS times_sparse(const vector<pair<ll, mint>> &x, ll deg = -1) const
{
if (deg == -1)
deg = size();
FPS ret(deg);
for (auto [i, a] : x)
{
ret += ((*this) << i).pre(deg) * a;
}
return ret;
}
static FPS inv_sparse(const vector<pair<ll, mint>> &x, ll deg)
{
if (deg == 0)
return {};
FPS ret(deg);
mint a = x[0].second.inv();
ret[0] = a;
rep(i, 1, deg)
{
rep(j, 1, x.size())
{
if (x[j].first > i)
break;
ret[i] -= ret[i - x[j].first] * x[j].second;
}
ret[i] *= a;
}
return ret;
}
static FPS log_sparse(const vector<pair<ll, mint>> &x, ll deg)
{
if (deg == 0)
return {};
FPS ret = inv_sparse(x, deg - 1);
vector<pair<ll, mint>> y = {x.begin() + 1, x.end()};
rep(i, y.size())
{
y[i].second *= y[i].first;
y[i].first--;
}
return ret.times_sparse(y).integral();
}
static FPS exp_sparse(const vector<pair<ll, mint>> &x, ll deg)
{
if (deg == 0)
return {};
FPS ret(deg);
ret[0] = 1;
rep(i, 1, deg)
{
rep(j, x.size())
{
if (x[j].first > i)
break;
ret[i] += ret[i - x[j].first] * x[j].second * x[j].first;
}
ret[i] *= mint_factorials::mint_inv(i);
}
return ret;
}
static FPS _pow_sparse(const vector<pair<ll, mint>> &x, mint m, ll deg)
{
if (deg == 0)
return {};
FPS ret(deg);
ret[0] = 1;
rep(i, 1, deg)
{
rep(j, 1, x.size())
{
if (x[j].first > i)
break;
ret[i] += ret[i - x[j].first] * x[j].second * ((m + 1) * x[j].first - i);
}
ret[i] *= mint_factorials::mint_inv(i);
}
return ret;
}
static FPS pow_sparse(vector<pair<ll, mint>> x, ll m, ll deg)
{
if (deg == 0)
return {};
if (m == 0)
{
FPS ret(deg);
return ret + 1;
}
if (x.empty())
return FPS(deg);
ll c = x[0].first;
if (c >= (deg + m - 1) / m)
return FPS(deg);
mint a = x[0].second;
mint ai = a.inv();
rep(i, x.size())
{
x[i].first -= c;
x[i].second *= ai;
}
return ((_pow_sparse(x, m, deg - c * m) * a.pow(m)) << (c * m));
}
static FPS _sqrt_sparse(vector<pair<ll, mint>> x, ll deg)
{
return _pow_sparse(x, mint_factorials::mint_inv(2), deg);
}
static FPS sqrt_sparse(vector<pair<ll, mint>> x, ll deg)
{
if (deg == 0)
return {};
if (x.empty())
return FPS(deg);
ll c = x[0].first;
if (c & 1)
return {};
mint a = x[0].second;
mint b = modsqrt(a);
if (b == -1)
return {};
rep(i, x.size())
{
x[i].first -= c;
x[i].second /= a;
}
return ((_sqrt_sparse(x, deg - c / 2) * b) << (c / 2));
}
mint prefix_product(ll k) const
{
ll deg = size() - 1;
if (k <= deg + 1)
{
mint ret = 1;
rep(i, k) ret *= eval(i);
return ret;
}
ll b = 0;
while ((((deg << (b + 1)) + 1) << (b + 1)) < k)
b++;
// mint B = (1ll << b);
vm G(deg + 1);
rep(i, deg + 1)
{
G[i] = eval(i << b);
}
mint twoi = mint(2).inv();
mint t = 1;
rep(i, b)
{
t *= twoi;
vm G1 = shift_sampling(G, t, (deg << i) + 1);
vm G2 = shift_sampling(G, (deg << i) + 1, deg << i);
vm G3 = shift_sampling(G, (deg << i) + 1 + t, deg << i);
rep(j, G.size()) G[j] *= G1[j];
rep(j, G2.size()) G.push_back(G2[j] * G3[j]);
}
mint ret = 1;
G = shift_sampling(G, 0, k >> b);
for (mint a : G)
ret *= a;
rep(i, k & ~((1ll << b) - 1), k) ret *= eval(i);
return ret;
}
static mint Bostan_Mori(FPS F, FPS G, ll k)
{
while (k)
{
if (F.empty())
return 0;
FPS G2(G.size());
rep(i, G.size()) G2[i] = (i & 1 ? -G[i] : G[i]);
FPS F2 = F * G2;
G2 *= G;
rep(i, G.size()) G[i] = G2[i * 2];
if (k & 1)
{
F.resize(F2.size() / 2);
rep(i, F.size()) F[i] = F2[i * 2 + 1];
}
else
{
F.resize((F2.size() + 1) / 2);
rep(i, F.size()) F[i] = F2[i * 2];
}
k >>= 1;
}
return F[0] / G[0];
}
// returns [x^{N-1}]] GF^i for i=0,...,M-1
static FPS PowerProjection(FPS F, FPS G, ll N, ll M)
{
static constexpr int g = internal::primitive_root<mint::mod()>;
static int cnt2 = __builtin_ctz(mint::mod() - 1);
static constexpr int inv2 = (mint::mod() + 1) >> 1;
static mint e2 = mint(g).pow((mint::mod() - 1) >> cnt2), ie2 = e2.inv();
static bool first = true;
static vm es(cnt2 + 1), ies(cnt2 + 1), sum_ie(cnt2 + 1);
if (first)
{
first = false;
es[cnt2] = e2, ies[cnt2] = ie2;
rrep(i, cnt2) es[i] = es[i + 1] * es[i + 1], ies[i] = ies[i + 1] * ies[i + 1];
mint e = 1;
rep(i, cnt2 + 1)
{
sum_ie[i] = ies[i] * e;
e *= es[i];
}
}
ll n = 64 - __builtin_clzll(N - 1);
G = G << ((1ll << n) - N);
N = 1ll << n;
F.resize(N);
G.resize(N);
ll n2 = 0, N2 = 1;
mint invN = mint(N).inv(), invN2 = 1;
vector<FPS> P = {G}, Q = {-F};
while (N > 1)
{
vector<FPS> NP(N2 << 1, FPS(N << 1, 0)), NQ(N2 << 1, FPS(N << 1, 0));
rep(i, N2) rep(j, N) NP[i][j] = P[i][j], NQ[i][j] = Q[i][j];
rep(j, N)
{
vm p(N2), q(N2);
rep(i, N2) p[i] = P[i][j], q[i] = Q[i][j];
internal::butterfly_inv(p), internal::butterfly_inv(q);
mint e = 1;
rep(i, N2)
{
p[i] *= e;
q[i] *= e;
e *= es[n2 + 1];
}
internal::butterfly(p), internal::butterfly(q);
rep(i, N2) NP[i + N2][j] = p[i] * invN2, NQ[i + N2][j] = q[i] * invN2;
}
P = vector<FPS>(N2 << 1, FPS(N)), Q = vector<FPS>(N2 << 1, FPS(N));
rep(i, N2 << 1)
{
NQ[i][0] += (i & N2 ? -1 : 1);
internal::butterfly(NP[i]), internal::butterfly(NQ[i]);
mint e = 1;
rep(j, N)
{
Q[i][j] = (NQ[i][2 * j] * NQ[i][2 * j + 1]) * invN;
P[i][j] = (NP[i][2 * j] * NQ[i][2 * j + 1] - NP[i][2 * j + 1] * NQ[i][2 * j]) * invN * inv2 * e;
e *= -sum_ie[__builtin_ctz(~(unsigned int)(j)) + 2];
}
internal::butterfly_inv(P[i]), internal::butterfly_inv(Q[i]);
P[i].resize(N >> 1);
Q[i].resize(N >> 1);
Q[i][0]--;
}
N >>= 1;
n--;
invN *= 2;
N2 <<= 1;
n2++;
invN2 *= inv2;
}
FPS p(N2), q(N2);
rep(i, N2) p[i] = P[i][0] * invN2, q[i] = Q[i][0] * invN2;
internal::butterfly_inv(p);
internal::butterfly_inv(q);
q.push_back(1);
reverse(p.begin(), p.end());
reverse(q.begin(), q.end());
p.resize(M);
q.resize(M);
return p / q;
}
static FPS compositional_inverse(FPS F, ll m = -1)
{
if (m == -1)
m = F.size();
assert(F.size() > 1);
assert(F[0].val() == 0);
mint c = F[1].inv();
F *= c;
FPS G = PowerProjection(F, {1}, m, m);
rep(i, 1, m) G[i] *= mint_factorials::mint_inv(i) * (m - 1);
reverse(G.begin(), G.end());
G.pop_back();
G = G.pow((-mint_factorials::mint_inv(m - 1)).val()) << 1;
mint e = c;
rep(i, 1, m)
{
G[i] *= e;
e *= c;
}
return G;
}
// return F(G(X))
static FPS composition(FPS F, FPS G, ll m = -1)
{
static constexpr int g = internal::primitive_root<mint::mod()>;
static int cnt2 = __builtin_ctz(mint::mod() - 1);
static constexpr int inv2 = (mint::mod() + 1) >> 1;
static mint e2 = mint(g).pow((mint::mod() - 1) >> cnt2), ie2 = e2.inv();
static bool first = true;
static vm es(cnt2 + 1), ies(cnt2 + 1), sum_ie(cnt2 + 1);
if (first)
{
first = false;
es[cnt2] = e2, ies[cnt2] = ie2;
rrep(i, cnt2) es[i] = es[i + 1] * es[i + 1], ies[i] = ies[i + 1] * ies[i + 1];
mint e = 1;
rep(i, cnt2 + 1)
{
sum_ie[i] = ies[i] * e;
e *= es[i];
}
}
if (m == -1)
m = max(F.size(), G.size());
if (F.empty())
return FPS(1, m);
if (G.empty())
{
FPS ret(1, F[0]);
ret.resize(m);
return ret;
}
if (m == 1)
return {F[0]};
ll n = 64 - __builtin_clzll(m - 1);
ll N = 1ll << n;
F.resize(N);
G.resize(N);
ll n2 = 0, N2 = 1;
mint invN = mint(N).inv(), invN2 = 1;
vector<FPS> Q = {-G};
vector<vector<FPS>> NQ(n);
while (N > 1)
{
NQ[n2] = vector<FPS>(N2 << 1, FPS(N << 1, 0));
rep(i, N2) rep(j, N) NQ[n2][i][j] = Q[i][j];
rep(j, N)
{
vm q(N2);
rep(i, N2) q[i] = Q[i][j];
internal::butterfly_inv(q);
mint e = 1;
rep(i, N2)
{
q[i] *= e;
e *= es[n2 + 1];
}
internal::butterfly(q);
rep(i, N2) NQ[n2][i + N2][j] = q[i] * invN2;
}
Q = vector<FPS>(N2 << 1, FPS(N));
rep(i, N2 << 1)
{
NQ[n2][i][0] += (i & N2 ? -1 : 1);
internal::butterfly(NQ[n2][i]);
rep(j, N) Q[i][j] = (NQ[n2][i][2 * j] * NQ[n2][i][2 * j + 1]) * invN;
internal::butterfly_inv(Q[i]);
Q[i].resize(N >> 1);
Q[i][0]--;
}
N >>= 1;
n--;
invN *= 2;
N2 <<= 1;
n2++;
invN2 *= inv2;
}
FPS q(N2);
rep(i, N2) q[i] = Q[i][0] * invN2;
internal::butterfly_inv(q);
q.push_back(1);
reverse(q.begin(), q.end());
q = q.inv(N2);
reverse(q.begin(), q.end());
F *= q;
F[0] = F.back();
rep(i, 1, N2) F[i] = F[i + N2 - 2];
F.resize(N2);
internal::butterfly(F);
vector<FPS> P(N2, FPS(1));
rep(i, N2) P[i][0] = F[i] * invN2;
while (N2 > 1)
{
N <<= 1;
n++;
invN *= inv2;
N2 >>= 1;
n2--;
invN2 *= 2;
vector<FPS> NP(N2 << 1, FPS(N << 1, 0));
rep(i, N2 << 1)
{
P[i].resize(N);
reverse(P[i].begin() + 1, P[i].end());
internal::butterfly(P[i]);
mint e = 1;
rep(j, N)
{
NP[i][2 * j] = P[i][j] * NQ[n2][i][2 * j + 1] * invN * inv2 * e;
NP[i][2 * j + 1] = -P[i][j] * NQ[n2][i][2 * j] * invN * inv2 * e;
e *= -sum_ie[__builtin_ctz(~(unsigned int)(j)) + 2];
}
internal::butterfly_inv(NP[i]);
reverse(NP[i].begin() + 1, NP[i].end());
}
P = vector<FPS>(N2, FPS(N, 0));
rep(i, N2) rep(j, N) P[i][j] = NP[i][j];
rep(j, N)
{
vector<mint> p(N2);
rep(i, N2) p[i] = NP[i + N2][j] * invN2;
internal::butterfly_inv(p);
mint e = es[n2 + 1];
rep(i, 1, N2)
{
p[N2 - i] *= e;
e *= es[n2 + 1];
}
internal::butterfly(p);
rep(i, N2) P[i][j] += p[i];
}
}
reverse(P[0].begin(), P[0].end());
return P[0].pre(m);
}
};
struct FPS_Mat
{
FPS a, b, c, d;
FPS_Mat() = default;
FPS_Mat(const FPS &a0, const FPS &b0, const FPS &c0, const FPS &d0)
: a(a0), b(b0), c(c0), d(d0) {};
FPS_Mat &operator*=(const FPS_Mat &r)
{
FPS A = a * r.a + b * r.c;
FPS B = a * r.b + b * r.d;
FPS C = c * r.a + d * r.c;
FPS D = c * r.b + d * r.d;
A.shrink();
B.shrink();
C.shrink();
D.shrink();
swap(A, a);
swap(B, b);
swap(C, c);
swap(D, d);
return *this;
}
static FPS_Mat I() { return {{1}, {}, {}, {1}}; }
FPS_Mat operator*(const FPS_Mat &r) const
{
return FPS_Mat(*this) *= r;
}
pair<FPS, FPS> app(const pair<FPS, FPS> &p) const
{
FPS q1 = a * p.first + b * p.second;
FPS q2 = c * p.first + d * p.second;
q1.shrink();
q2.shrink();
return {q1, q2};
}
};
void Naive_GCD(FPS_Mat &M, pair<FPS, FPS> &p, mint &res, int pre)
{
pair<FPS, FPS> qr = p.first.div(p.second);
FPS C = M.a - M.c * qr.first;
FPS D = M.b - M.d * qr.first;
res *= p.second.back().pow(p.first.size() - qr.second.size());
res *= ((~(p.first.size() ^ pre) & 1) && (~(p.second.size() ^ pre) & 1) ? -1 : 1);
C.shrink();
D.shrink();
swap(C, M.c);
swap(D, M.d);
swap(C, M.a);
swap(D, M.b);
p = {p.second, qr.second};
}
FPS_Mat Half_GCD(pair<FPS, FPS> p, mint &res, int pre = 0)
{
int n = p.first.size();
int m = p.second.size();
int k = (n + 1) / 2;
if (m <= k)
return FPS_Mat::I();
auto R1 = Half_GCD({p.first >> k, p.second >> k}, res, pre + k);
p = R1.app(p);
if ((int)p.second.size() < k)
res *= p.first.back().pow(k - p.second.size());
if ((int)p.second.size() <= k)
return R1;
Naive_GCD(R1, p, res, pre);
if ((int)p.second.size() <= k)
return R1;
int l = p.first.size() - 1;
int j = 2 * k - l;
p.first = p.first >> j;
p.second = p.second >> j;
return Half_GCD(p, res, pre + j) * R1;
}
pair<FPS_Mat, mint> GCD_Mat(const FPS &f, const FPS &g)
{
pair<FPS, FPS> p{f, g};
mint res = 1;
p.first.shrink();
p.second.shrink();
int n = p.first.size();
int m = p.second.size();
if (n < m)
{
auto [ret, res0] = GCD_Mat(p.second, p.first);
if ((~n & 1) && (~m & 1))
res0 *= -1;
swap(ret.a, ret.b);
swap(ret.c, ret.d);
return {ret, res0 * res};
}
FPS_Mat ret = FPS_Mat::I();
while (true)
{
auto m1 = Half_GCD(p, res);
p = m1.app(p);
if (p.second.empty())
return {m1 * ret, res};
Naive_GCD(m1, p, res, 0);
if (p.second.empty())
return {m1 * ret, res};
ret = m1 * ret;
}
}
FPS PolyGCD(const FPS &f, const FPS &g)
{
pair<FPS, FPS> p{f, g};
FPS_Mat m = GCD_Mat(f, g).first;
p = m.app(p);
if (!p.first.empty())
p.first *= p.first.back().inv();
return p.first;
}
pair<bool, FPS> PolyInv(const FPS &f, const FPS &g)
{
pair<FPS, FPS> p{f, g};
FPS_Mat m = GCD_Mat(f, g).first;
FPS gcd_ = m.app(p).first;
if ((int)gcd_.size() != 1)
return {false, FPS()};
return {true, m.a * gcd_.back().inv()};
}
mint Resolutant(const FPS &f, const FPS &g)
{
pair<FPS, FPS> p{f, g};
auto [m, res] = GCD_Mat(f, g);
FPS gcd_ = m.app(p).first;
if ((int)gcd_.size() != 1)
return 0;
return {res * gcd_.back().inv()};
}
vector<mint> FindRoot(FPS F)
{
static random_device rnd;
static mt19937 mt(rnd());
static uniform_int_distribution<> randp(0, mint::mod() - 1);
static constexpr long long q = (mint::mod() - 1) >> 1;
F.shrink();
int n = F.size();
if (n <= 1)
return {};
if (n == 2)
return {-F[0] / F[1]};
FPS Frevinv = F.rev().inv();
static function<void(FPS &, FPS &, FPS &)> modf = [&](FPS &G, FPS &f, FPS &finv)
{
if (G.size() >= f.size())
{
FPS Q = (G.rev() * finv).pre(G.size() - f.size() + 1).rev();
G -= f * Q;
G.shrink();
}
};
FPS G = {1};
FPS G2 = {0, 1};
long long k = mint::mod();
while (k)
{
if (k & 1)
{
G *= G2;
modf(G, F, Frevinv);
}
G2 *= G2;
modf(G2, F, Frevinv);
k = k >> 1;
}
G -= {0, 1};
F = PolyGCD(F, G);
n = F.size();
if (n <= 1)
return {};
if (n == 2)
return {-F[0] / F[1]};
Frevinv = F.rev().inv();
vector<mint> ans = {};
vector<FPS> FS = {F};
while (!FS.empty())
{
mint s = randp(mt);
vector<FPS> NFS;
for (FPS &f : FS)
{
if (f.eval(s) == 0)
ans.push_back(s);
FPS frevinv = f.rev().inv();
FPS g = {1};
FPS g2 = {-s, 1};
long long k = q;
while (k)
{
if (k & 1)
{
g *= g2;
modf(g, f, frevinv);
}
g2 *= g2;
modf(g2, f, frevinv);
k = k >> 1;
}
g -= 1;
FPS f1 = PolyGCD(f, g);
g += 2;
FPS f2 = PolyGCD(f, g);
if (f1.size() == 2)
ans.push_back(-f1[0]);
else if (f1.size() > 2)
NFS.push_back(f1);
if (f2.size() == 2)
ans.push_back(-f2[0]);
else if (f2.size() > 2)
NFS.push_back(f2);
}
FS = NFS;
}
return ans;
}
int main()
{
std::cin.tie(nullptr), std::ios_base::sync_with_stdio(false);
readll(N, M, K, L);
if(N == K || L == M)
{
dame(0);
}
mint ans = 0;
// FPS f(M + 1, 1);
{
for(ll i = 0; i*(M + 1) <= M*K && i <= N; i++)
{
ans += ((i&1) == 0 ? 1 : (-1))*mint_factorials::nCk(N, i)*mint_factorials::nCk(M*K-(i*(M+1)) + N - 1, N - 1);
}
}
// debug(ans);
{
// FPS g(L, 1);
// ans -= (g * f.pow(N-1, M*K + 1))[M*K];
REP(k, 0, L-1)
{
for(ll i = 0; i*(M + 1) <= M*K && i <= N - 1; i++) ans -= ((i&1) == 0 ? 1 : (-1))*mint_factorials::nCk(N-1, i)*mint_factorials::nCk(M*K-k -(i*(M+1)) + N - 2, N - 2);
}
}
// debug(ans);
REP(k, L, M)
{
mint tmp = 0;
for(ll i = 0; i*(k+1) <= M*K - k && i <= N - 1; i++)
{
tmp += ((i&1) == 0 ? 1 : (-1))*mint_factorials::nCk(N-1, i)*mint_factorials::nCk(M*K-k-(i*(k+1)) + N - 2, N - 2);
}
// debug(tmp);
ans -= tmp;
}
out(ans);
}