1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
|
#include<bits/stdc++.h>
using namespace std;
namespace exfunc{
#ifdef _WIN32
#include<windows.h>
#include<psapi.h>
#else
#include<sys/resource.h>
#endif
#define fs first
#define sc second
#define pb push_back
typedef unsigned int ui;
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int,int> pii;
template<typename T>void chkmin(T& a,T b){if(b<a)a=b;}
template<typename T>void chkmax(T& a,T b){if(b>a)a=b;}
struct dsu{
public:
dsu()=delete;
dsu(int _n):n(_n),f(_n+1,-1){}
int merge(int a,int b){int x=find(a),y=find(b);if(x==y)return x;if(-f[x]<-f[y])swap(x,y);return f[x]+=f[y],f[y]=x;}
bool same(int x,int y){return find(x)==find(y);}
int find(int x){if(f[x]<0)return x;return f[x]=find(f[x]);}
int size(int x){return -f[find(x)];}
private:
int n;
vector<int>f;
};
template<class T>struct BIT{
public:
BIT()=delete;
BIT(int _n):n(_n),f(_n+1){};
BIT(vector<T>a):n(a.size()-1),f(a.size(),0){f=a;for(int i=1;i<=n;i++){int fa=i+(i&-i);if(fa<=n)f[fa]+=f[i];}}
void add(int p,T x){while(p<=n){f[p]+=x;p+=p&-p;}}
T sum(int p){T s=0;while(p>0){s+=f[p];p-=p&-p;}return s;}
T sum(int l,int r){return sum(r)-sum(l-1);}
private:
int n;
vector<T>f;
};
namespace math{
constexpr ll smod(ll x,ll m){x%=m;if(x<0)x+=m;return x;}
constexpr ll qmi(ll x,ll n,int m){if(m==1)return 0;ui _m=static_cast<ui>(m);ull r=1,y=smod(x,m);while(n){if(n&1)r=(r*y)%_m;y=(y*y)%_m;n>>=1;}return static_cast<ll>(r);}
constexpr ll qmi(ll x,ll n,ll m){if(m<=INT_MAX)return qmi(x,n,(int)m);if(m==1)return 0;ull _m=static_cast<ull>(m);__int128 r=1,y=smod(x,m);while(n){if(n&1)r=(r*y)%_m;y=(y*y)%_m;n>>=1;}return static_cast<ll>(r);}
constexpr bool isprime(int n){if(n<=1)return 0;if(n==2||n==7||n==61)return 1;if(n%2==0)return 0;ll d=n-1;while(d%2==0)d/=2;constexpr ll b[3]={2,7,61};for(ll a:b){ll t=d;ll y=qmi(a,t,n);while(t!=n-1&&y!=1&&y!=n-1){y=(y*y)%n;t<<=1;}if(y!=n-1&&t%2==0)return 0;}return 1;}
constexpr bool isprime(ll n){if(n<=INT_MAX)return isprime((int)n);if(n<=1)return 0;if(n==2)return 1;if(n%2==0)return 0;ll d=n-1;while(d%2==0)d/=2;constexpr ll b[7]={2,325,9375,28178,450775,9780504,1795265022};for(ll a:b){ll t=d;ll y=qmi(a,t,n);while(t!=n-1&&y!=1&&y!=n-1){y=(y*y)%n;t<<=1;}if(y!=n-1&&t%2==0)return 0;}return 1;}
constexpr pair<ll,ll>invgcd(ll a,ll b){a=smod(a,b);if(a==0)return {b,0};ll s=b,t=a,m0=0,m1=1;while(t){ll u=s/t;s-=t*u,m0-=m1*u;swap(s,t),swap(m0,m1);}if(m0<0)m0+=b/s;return {s,m0};}
constexpr ll invmod(ll x,ll m){return invgcd(x,m).sc;}
constexpr ll gcd(ll x,ll y){if(!x)return y;if(!y)return x;ll t=__builtin_ctzll(x|y);x>>=__builtin_ctzll(x);while(y){y>>=__builtin_ctzll(y);if(x>y)swap(x,y);y-=x;}return x<<t;}
int ceilpow2(int n){if(n<=1)return 0;return 32-__builtin_clz(n-1);}
int bsf(ui n){return __builtin_ctz(n);}
vector<int>sieve(int n){vector<int>res;vector<bool>isp(n+1,1);isp[0]=isp[1]=false;for(int j=4;j<=n;j+=2)isp[j]=false;for(int i=3;i*i<=n;i+=2)if(isp[i])for(int j=i*i;j<=n;j+=i*2)isp[j]=false;if(n>=2)res.pb(2);for(int j=3;j<=n;j+=2)if(isp[j])res.pb(j);return res;}
pair<ll,ll>crt(const vector<ll>&r,const vector<ll>&m) {int n=int(r.size());ll r0=0,m0=1;for(int i=0;i<n;i++){ll r1=smod(r[i],m[i]),m1=m[i];if(m0<m1){swap(r0,r1);swap(m0,m1);}if(m0%m1==0){if(r0%m1!=r1)return {0,0};continue;}ll g,im;tie(g,im)=invgcd(m0, m1);ll u1=(m1/g);if((r1-r0)%g)return {0,0};ll x=(r1-r0)/g%u1*im%u1;r0+=x*m0;m0*=u1;if(r0<0)r0+=m0;}return {r0,m0};}
};
struct barrett{
ui _m;ull im;
barrett(ui m):_m(m),im((ull)(-1)/m+1){}
ui umod()const{return _m;}
ui mul(ui a,ui b)const{ull z=a;z*=b;ull x=(ull)(((unsigned __int128)(z)*im)>>64);ui v=(ui)(z-x*_m);if(_m<=v)v+=_m;return v;}
};
template<int m>struct static_modint{
using mi=static_modint;
static constexpr int mod(){return m;}
static constexpr ui umod(){return static_cast<ui>(m);}
static constexpr bool prime=math::isprime(m);
static mi raw(int v){mi x;x._v=static_cast<ui>(v);return x;}
static_modint():_v(0){}
template<class T,enable_if_t<is_signed<T>::value,int> =0>
static_modint(T v){ll x=static_cast<ll>(v%static_cast<ll>(umod()));if(x<0)x+=umod();_v=static_cast<ui>(x);}
template<class T,enable_if_t<is_unsigned<T>::value,int> =0>
static_modint(T v){_v=static_cast<ui>(v%umod());}
static_modint(bool v){_v=static_cast<ui>(v%umod());}
ui val()const {return _v;}
mi&operator++(){++_v;if(_v==umod())_v=0;return *this;}
mi&operator--(){if(_v==0)_v=umod();--_v;return *this;}
mi operator++(int){mi tmp=*this;++*this;return tmp;}
mi operator--(int){mi tmp=*this;--*this;return tmp;}
mi&operator+=(const mi&rhs){_v+=rhs._v;if(_v>=umod())_v-=umod();return *this;}
mi&operator-=(const mi&rhs){if(_v<rhs._v)_v+=umod();_v-=rhs._v;return *this;}
mi&operator*=(const mi&rhs){ull z=_v;z*=rhs._v;_v=static_cast<ui>(z%umod());return *this;}
mi&operator/=(const mi&rhs){return *this=*this*rhs.inv();}
mi operator+()const {return *this;}
mi operator-()const {return mi()-*this;}
mi pow(ll n)const{mi x=*this,r=1;while(n){if(n&1)r*=x;x*=x;n>>=1;}return r;}
mi inv()const{if constexpr(prime)return pow(umod()-2);else{auto eg=math::invgcd(_v,m);return raw(static_cast<int>(eg.second));}}
friend mi operator+(const mi& lhs,const mi& rhs){return mi(lhs)+=rhs;}
friend mi operator-(const mi& lhs,const mi& rhs){return mi(lhs)-= rhs;}
friend mi operator*(const mi& lhs,const mi& rhs){return mi(lhs)*=rhs;}
friend mi operator/(const mi& lhs,const mi& rhs){return mi(lhs)/=rhs;}
friend bool operator==(const mi& lhs,const mi& rhs){return lhs._v==rhs._v;}
friend bool operator!=(const mi& lhs,const mi& rhs){return lhs._v!=rhs._v;}
private:
ui _v;
};
template<int id>struct dynamic_modint{
using mi=dynamic_modint;
static int mod(){return static_cast<int>(bt.umod());}
static void set_mod(int m){bt=barrett(m);}
static mi raw(int v){mi x;x._v=static_cast<ui>(v);return x;}
dynamic_modint():_v(0){}
template<class T,enable_if_t<is_integral<T>::value,int> =0>
dynamic_modint(T v){if constexpr(is_signed<T>::value){ll x=static_cast<ll>(v%static_cast<ll>(mod()));if(x<0)x+=mod();_v=static_cast<ui>(x);}else _v=static_cast<ui>(v%mod());}
ui val()const {return _v;}
mi&operator++(){if(++_v==umod())_v=0;return *this;}
mi&operator--(){if(_v==0)_v=umod();--_v;return *this;}
mi operator++(int){mi t=*this;++*this;return t;}
mi operator--(int){mi t=*this;--*this;return t;}
mi&operator+=(const mi& rhs){_v+=rhs._v;if(_v>=umod())_v-=umod();return *this;}
mi&operator-=(const mi& rhs) {if(_v<rhs._v)_v+=umod();_v-=rhs._v;return *this;}
mi&operator*=(const mi& rhs){_v=bt.mul(_v,rhs._v);return *this;}
mi&operator/=(const mi& rhs){return *this=*this*rhs.inv();}
mi operator+()const {return *this;}
mi operator-()const {return mi()-*this;}
mi pow(ll n)const{mi x=*this,r=1;while(n){if(n&1)r*=x;x*=x;n>>=1;}return r;}
mi inv()const{auto eg=math::invgcd(_v,mod());return raw(static_cast<int>(eg.second));}
friend mi operator+(const mi& lhs,const mi& rhs){return mi(lhs)+=rhs;}
friend mi operator-(const mi& lhs,const mi& rhs){return mi(lhs)-=rhs;}
friend mi operator*(const mi& lhs,const mi& rhs){return mi(lhs)*=rhs;}
friend mi operator/(const mi& lhs,const mi& rhs){return mi(lhs)/=rhs;}
friend bool operator==(const mi& lhs,const mi& rhs){return lhs._v==rhs._v;}
friend bool operator!=(const mi& lhs,const mi& rhs){return lhs._v!=rhs._v;}
private:
ui _v;
static barrett bt;
static ui umod(){return bt.umod();}
};
template<int id>barrett dynamic_modint<id>::bt=998244353;
template<class S,S(*op)(S,S),S(*e)()>struct segtree{
segtree()=delete;
explicit segtree(int n):_n(n){log=math::ceilpow2(_n);size=1<<log;d=vector<S>(2*size,e());}
explicit segtree(const vector<S>&v):_n(static_cast<int>(v.size())-1){log=math::ceilpow2(_n);size=1<<log;d.resize(2*size);if(_n>0){copy(v.begin()+1,v.end(),d.begin()+size);fill(d.begin()+size+_n,d.end(),e());for(int i=size-1;i>=1;i--)update(i);}else fill(d.begin(),d.end(),e());}
void set(int p,const S& x){p=p+size-1;d[p]=x;for(int i=1;i<=log;i++)update(p>>i);}
const S& get(int p)const{return d[p+size-1];}
S prod(int l,int r)const{if(l>r)return e();S sml=e(),smr=e();for(int L=l+size-1,R=r+size;L<R;L>>=1,R>>=1){if(L&1)sml=op(sml,d[L++]);if(R&1)smr=op(d[--R],smr);}return op(sml,smr);}
const S& all_prod()const{return d[1];}
template<bool(*f)(const S&)>int max_right(int l)const{return max_right(l,[](const S& x){return f(x);});}
template<class F>int max_right(int l,F f)const{if(l>_n)return _n;int cur=l+size-1;S sm=e();do{cur>>=math::bsf(static_cast<ui>(cur));if(!f(op(sm,d[cur]))){while(cur<size){cur=(2*cur);if(f(op(sm,d[cur]))){sm=op(sm,d[cur]);cur++;}}return cur-size;}sm=op(sm,d[cur]);cur++;}while((cur&(cur-1)));return _n;}
template<bool(*f)(const S&)>int min_left(int r)const{return min_left(r,[](const S& x){return f(x);});}
template<class F>int min_left(int r,F f)const{if(r<1)return 1;int cur=r+size-1;S sm=e();do{while(cur>1&&(cur&1))cur>>=1;if(!f(op(d[cur],sm))){while(cur<size){cur=(2*cur+1);if(f(op(d[cur],sm))){sm=op(d[cur],sm);cur--;}}return cur-size+2;}sm=op(d[cur],sm);cur--;}while((cur&(cur-1)));return 1;}
private:
int _n,size,log;
vector<S> d;
void update(int k){d[k]=op(d[2*k],d[2*k+1]);}
};
template<class S,S(*op)(S,S),S(*e)(),class F,S(*mapping)(F,S),F(*composition)(F,F),F(*id)()>struct lazy_segtree{
lazy_segtree()=delete;
explicit lazy_segtree(int n):_n(n){log=math::ceilpow2(_n);size=1<<log;d=vector<S>(2*size,e());lz=vector<F>(size,id());}
explicit lazy_segtree(const vector<S>& v):_n(static_cast<int>(v.size())-1){log=math::ceilpow2(_n);size=1<<log;d.resize(2*size);lz=vector<F>(size,id());if(_n>0){copy(v.begin()+1,v.end(),d.begin()+size);fill(d.begin()+size+_n,d.end(),e());for(int i=size-1;i>=1;i--)update(i);}else fill(d.begin(), d.end(), e());}
void set(int p,const S& x){p=p+size-1;for(int i=log;i>=1;i--)push(p>>i);d[p]=x;for(int i=1;i<=log;i++)update(p>>i);}
S get(int p){p=p+size-1;for(int i=log;i>=1;i--)push(p>>i);return d[p];}
S prod(int l,int r){if(l>r)return e();int L=l+size-1,R=r+size;for(int i=log;i>=1;i--){if(((L>>i)<<i)!=L)push(L>>i);if(((R>>i)<<i)!=R)push((R-1)>>i);}S sml=e(),smr=e();for(;L<R;L>>=1,R>>=1){if(L&1)sml=op(sml,d[L++]);if(R&1)smr=op(d[--R],smr);}return op(sml,smr);}
const S& all_prod()const{return d[1];}
void apply(int p,const F& f){p=p+size-1;for(int i=log;i>=1;i--)push(p>>i);d[p]=mapping(f,d[p]);for(int i=1;i<=log;i++)update(p>>i);}
void apply(int l, int r, const F& f){if(l>r)return;int L=l+size-1,R=r+size;for(int i=log;i>=1;i--){if(((L>>i)<<i)!=L)push(L>>i);if(((R>>i)<<i)!=R)push((R-1)>>i);}{int l2=L,r2=R;for(;L<R;L>>= 1,R>>=1){if(L&1)all_apply(L++,f);if(R&1)all_apply(--R,f);}L=l2;R=r2;}for(int i=1;i<=log;i++){if(((L>>i)<<i)!=L)update(L>>i);if(((R>>i)<<i)!=R)update((R-1)>>i);}}
template<bool(*f)(const S&)>int max_right(int l){return max_right(l,[](const S& x){return f(x);});}
template<class G>int max_right(int l,G g){if(l>_n)return _n;int cur=l+size-1;for(int i=log;i>=1;i--)push(cur>>i);S sm=e();do{while((cur&1)==0)cur>>=1;if(!g(op(sm,d[cur]))){while(cur<size){push(cur);cur=(2*cur);if(g(op(sm,d[cur]))){sm=op(sm,d[cur]);cur++;}}return cur-size;}sm=op(sm,d[cur]);cur++;}while(cur&(cur-1));return _n;}
template<bool(*f)(const S&)>int min_left(int r){return min_left(r,[](const S& x){return f(x);});}
template<class G> int min_left(int r,G g){if(r<1)return 1;int cur=r+size-1;for(int i=log;i>=1;i--)push(cur>>i);S sm=e();do{while(cur>1&&(cur&1))cur>>=1;if(!g(op(d[cur],sm))){while(cur<size){push(cur);cur=(2*cur+1);if(g(op(d[cur],sm))){sm=op(d[cur],sm);cur--;}}return cur-size+2;}sm=op(d[cur],sm);if((cur&-cur)==cur)break;cur--;}while(1);return 1;}
private:
int _n,size,log;
vector<S> d;
vector<F> lz;
void update(int k){d[k]=op(d[2*k],d[2*k+1]);}
void all_apply(int k,const F& f){d[k]=mapping(f,d[k]);if(k<size)lz[k]=composition(f,lz[k]);}
void push(int k){all_apply(2*k,lz[k]);all_apply(2*k+1,lz[k]);lz[k]=id();}
};
namespace IO{
void fastIO(){ios::sync_with_stdio(0);cin.tie(0);}
void openf(string fileName){freopen((fileName+".in").c_str(),"r",stdin);freopen((fileName+".out").c_str(),"w",stdout);}
#define gc() getchar()
int read(){int k=0,f=1,c=gc();for(;!isdigit(c);c=gc())if(c=='-')f=-1;for(;isdigit(c);c=gc())k=k*10+(c^48);return k*f;}
template<typename T>bool read(T& x){T k=0;int f=1,c=gc();for(;!isdigit(c);c=gc())if(c=='-')f=-1;for(;isdigit(c);c=gc())k=k*10+(c^48);x=k*f;return 1;}
#undef gc
uint8_t buf[1<<20],*p1,*p2;
#define gc() (p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<20,stdin),p1==p2)?0:*p1++)
int readf(){int k=0,f=1,c=gc();for(;!isdigit(c);c=gc())if(c=='-')f=-1;for(;isdigit(c);c=gc())k=k*10+(c^48);return k*f;}
template<typename T>bool readf(T& x){T k=0;int f=1,c=gc();if(!c)return 0;for(;!isdigit(c);c=gc())if(c=='-')f=-1;for(;isdigit(c);c=gc()){k=k*10+(c^48);if(!c)return 0;}x=k*f;return 1;}
#undef gc
template<typename T>void write(T x){if(x<0)putchar('-'),x=-x;if(x>9)write(x/10);putchar(x%10+'0');}
}
namespace rd{
mt19937 rand32(chrono::high_resolution_clock::now().time_since_epoch().count());
int randint(int l,int r){uniform_int_distribution<>di(l,r);return di(rand32);}
double random(){uniform_real_distribution<>di;return di(rand32);}
}
template<const int M>struct Comb_static{
using mi=static_modint<M>;
vector<mi>fact,invfact;
int N=0,unex=0;
Comb_static()=delete;
Comb_static(int n,int _unex):N(n),unex(_unex){fact.resize(N+1),invfact.resize(N+1);fact[0]=fact[1]=invfact[0]=1;for(int i=2;i<=N;i++)fact[i]=fact[i-1]*i;invfact[N]=math::invmod(fact[N].val(),M);for(int i=N-1;i>=1;i--)invfact[i]=invfact[i+1]*(i+1);}
ll comb(ll m,ll n){if(n<0||n>m)return unex;if(n==0||n==m)return 1;return (fact[m]*invfact[n]*invfact[m-n]).val();}
};
template<int id>struct Comb_dynamic{
using mi=dynamic_modint<id>;
vector<mi>fact,invfact;
int N=0,M=0,unex=0;
Comb_dynamic()=delete;
Comb_dynamic(int n,int mod,int _unex):N(n),M(mod),unex(_unex){fact.resize(N+1),invfact.resize(N+1);mi::set_mod(M);fact[0]=fact[1]=invfact[0]=1;for(int i=2;i<=N;i++)fact[i]=fact[i-1]*i;invfact[N]=math::invmod(fact[N].val(),M);for(int i=N-1;i>=1;i--)invfact[i]=invfact[i+1]*(i+1);}
ll comb(ll m,ll n){if(n<0||n>m)return unex;if(n==0||n==m)return 1;return (fact[m]*invfact[n]*invfact[m-n]).val();}
};
template<class T>struct Matrix{
using Ma=Matrix;
ui n,m;
vector<T>data;
Matrix()=delete;
Matrix(ui _n,ui _m):n(_n),m(_m),data(_n*_m,T()){}
static Ma I(ui s){Ma r(s,s);for(ui i=1;i<=s;++i)r[i][i]=T(1);return r;}
void set(ui i,ui j,const T& k){data[(i-1)*m+j-1]=k;}
struct P{T* p;T& operator[](ui c){return p[c-1];}};
struct CP{const T* p;const T& operator[](ui c)const{return p[c-1];}};
P operator[](ui r){return{data.data()+(r-1)*m};}
CP operator[](ui r)const{return{data.data()+(r-1)*m};}
Ma operator+(const Ma& o)const{Ma r(n,m);for(ui i=0;i<n*m;++i)r.data[i]=data[i]+o.data[i];return r;}
Ma operator-(const Ma& o)const{Ma r(n,m);for(ui i=0;i<n*m;++i)r.data[i]=data[i]-o.data[i];return r;}
Ma operator*(const Ma& o)const{Ma r(n,o.m);for(ui i=1;i<=n;++i)for(ui k=1;k<=m;++k){T v=(*this)[i][k];if(!z(v))for(ui j=1;j<=o.m;++j)r[i][j]+=v*o[k][j];}return r;}
Ma operator^(ull p)const{Ma r=I(n),b=*this;for(;p;p>>=1,b=b*b)if(p&1)r=r*b;return r;}
static bool z(T x){if constexpr(is_floating_point<T>::value)return abs(x)<1e-9;return x==T();}
static bool cmp(T a,T b){if constexpr(is_floating_point<T>::value)return abs(a)>abs(b);return !z(a)&&z(b);}
int rref(ui lim=0){if(!lim)lim=m;ui r=1;for(ui c=1;c<=lim&&r<=n;++c){ui p=r;for(ui i=r+1;i<=n;++i)if(cmp((*this)[i][c],(*this)[p][c]))p=i;if(z((*this)[p][c]))continue;if(p!=r)for(ui j=1;j<=m;++j)swap((*this)[r][j],(*this)[p][j]);T inv=T(1)/(*this)[r][c];for(ui j=c;j<=m;++j)(*this)[r][j]*=inv;for(ui i=1;i<=n;++i)if(i!=r&&!z((*this)[i][c])){T f=(*this)[i][c];for(ui j=c;j<=m;++j)(*this)[i][j]-=f*(*this)[r][j];}++r;}return r-1;}
T det()const{Ma a=*this;T res=1;for(ui i=1;i<=n;++i){ui p=i;for(ui j=i+1;j<=n;++j)if(cmp(a[j][i],a[p][i]))p=j;if(z(a[p][i]))return T();if(p!=i){for(ui j=1;j<=n;++j)swap(a[i][j],a[p][j]);res=-res;}res*=a[i][i];T inv=T(1)/a[i][i];for(ui j=i+1;j<=n;++j)if(!z(a[j][i])){T f=a[j][i]*inv;for(ui k=i;k<=n;++k)a[j][k]-=f*a[i][k];}}return res;}
pair<bool,Ma>inv()const{Ma a(n,2*n);for(ui i=1;i<=n;++i){for(ui j=1;j<=n;++j)a[i][j]=(*this)[i][j];a[i][n+i]=1;}if(a.rref(n)<n)return{0,Ma(0,0)};Ma r(n,n);for(ui i=1;i<=n;++i)for(ui j=1;j<=n;++j)r[i][j]=a[i][n+j];return{1,r};}
static int solve(Ma A,vector<T>b,vector<T>&ans){Ma a(A.n,A.m+1);for(ui i=1;i<=A.n;++i){for(ui j=1;j<=A.m;++j)a[i][j]=A[i][j];a[i][A.m+1]=b[i];}ui r=1;ans.assign(A.m+1,T());for(ui c=1;c<=A.m&&r<=A.n;++c){ui p=r;for(ui i=r+1;i<=A.n;++i)if(cmp(a[i][c],a[p][c]))p=i;if(z(a[p][c]))continue;if(p!=r)for(ui j=1;j<=A.m+1;++j)swap(a[r][j],a[p][j]);T inv=T(1)/a[r][c];for(ui j=c;j<=A.m+1;++j)a[r][j]*=inv;for(ui i=1;i<=A.n;++i)if(i!=r&&!z(a[i][c])){T f=a[i][c];for(ui j=c;j<=A.m+1;++j)a[i][j]-=f*a[r][j];}ans[c]=a[r][A.m+1];++r;}for(ui i=r;i<=A.n;++i)if(!z(a[i][A.m+1]))return -1;return r-1<A.m;}
};
template<class S,S(*op)(S,S),S(*e)()>struct ST{
ST()=delete;
int n,W;
vector<S>d;
ST(vector<S>a){if((n=a.size()-1)<=0)return;W=n+1;int K=__lg(n);d.resize((K+1)*W);copy(a.begin()+1,a.end(),d.begin()+1);for(int i=1;i<=K;++i){const S*p=d.data()+(i-1)*W;S*c=d.data()+i*W;for(int j=1,h=1<<(i-1),lim=n-(1<<i)+1;j<=lim;++j)c[j]=op(p[j],p[j+h]);}}
inline S prod(int l,int r)const{if(__builtin_expect(l>r,0))return e();int k=__lg(r-l+1);const S*p=d.data()+k*W;return op(p[l],p[r-(1<<k)+1]);}
};
class Perf{
using C=chrono::steady_clock;
C::time_point s,e;
public:void start(){s=C::now();}
void end(){e=C::now();}
void output(ostream& o=cerr){size_t m=0;
#ifdef _WIN32
PROCESS_MEMORY_COUNTERS p;GetProcessMemoryInfo(GetCurrentProcess(),&p,sizeof(p));m=p.PeakWorkingSetSize;
#else
rusage r;getrusage(0,&r);m=r.ru_maxrss*1024;
#endif
o<<chrono::duration_cast<chrono::microseconds>(e-s).count()<<"us "<<chrono::duration_cast<chrono::microseconds>(e-s).count()/1000.0<<"ms "<<chrono::duration_cast<chrono::microseconds>(e-s).count()/1000000.0<<"s\n"<<m<<"B "<<m/1024.0<<"KiB "<<m/1048576.0<<"MiB\n";}
};
template<class Cap>struct mf_graph{
public:
mf_graph():_n(0){}
mf_graph(int n):_n(n),g(n){}
int add_edge(int fr,int to,Cap cap){int m=pos.size(),fri=g[fr].size(),toi=g[to].size();pos.pb({fr,g[fr].size()});if(fr==to)toi++;g[fr].pb({to,toi,cap});g[to].pb({fr,fri,0});return m;}
struct ed{int from,to;Cap cap,flow;};
ed get_edge(int i){int m=pos.size();auto e=g[pos[i].fs][pos[i].sc];auto re=g[e.to][e.rev];return ed{pos[i].fs,e.to,e.cap+re.cap,re.cap};}
vector<ed>edges(){int m=pos.size();vector<ed>res;for(int i=0;i<m;i++)res.pb(get_edge(i));return res;}
void change_edge(int i,Cap nc,Cap nf){int m=pos.size();auto& e=g[pos[i].fs][pos[i].sc];auto& re=g[e.to][e.rev];e.cap=nc-nf;re.cap=nf;}
Cap flow(int s,int t){return flow(s,t,numeric_limits<Cap>::max());}
Cap flow(int s,int t,Cap fl){vector<int>lev(_n),it(_n);auto bfs=[&](){fill(lev.begin(),lev.end(),-1);lev[s]=0;queue<int>q;q.push(s);while(!q.empty()){int v=q.front();q.pop();for(auto e:g[v]){if(e.cap==0||lev[e.to]>=0)continue;lev[e.to]=lev[v]+1;if(e.to==t)return;q.push(e.to);}}};auto dfs=[&](auto sf,int v,Cap up){if(v==s)return up;Cap res=0;int levv=lev[v];for(int& i=it[v];i<g[v].size();i++){edg& e=g[v][i];if(levv<=lev[e.to]||g[e.to][e.rev].cap==0)continue;Cap d=sf(sf,e.to,min(up-res,g[e.to][e.rev].cap));if(d<=0)continue;g[v][i].cap+=d;g[e.to][e.rev].cap-=d;res+=d;if(res==up)break;}return res;};Cap fw=0;while(fw<fl){bfs();if(lev[t]==-1)break;fill(it.begin(),it.end(),0);while(fw<fl){Cap f=dfs(dfs,t,fl-fw);if(!f)break;fw+=f;}}return fw;}
vector<bool>min_cut(int s){vector<bool>vis(_n);queue<int>q;q.push(s);while(!q.empty()){int p=q.front();q.pop();vis[p]=true;for(auto e:g[p])if(e.cap&&!vis[e.to]){vis[e.to]=true;q.push(e.to);}}return vis;}
private:
int _n;
struct edg{int to,rev;Cap cap;};vector<pii>pos;vector<vector<edg>>g;
};
template<class Cap,class Cost>struct mcf_graph{
mcf_graph(){}
mcf_graph(int n):_n(n),g(n){}
int add_edge(int fr,int to,Cap cap,Cost cst){int m=pos.size();pos.pb({fr,g[fr].size()});int fri=g[fr].size(),toi=g[to].size();if(fr==to)toi++;g[fr].pb({to,toi,cap,cst});g[to].pb({fr,fri,0,-cst});return m;}
struct ed{int from,to;Cap cap,flow;Cost cost;};
ed get_edge(int i){int m=pos.size();auto e=g[pos[i].fs][pos[i].sc];auto re=g[e.to][e.rev];return {pos[i].fs,e.to,e.cap+re.cap,re.cap,e.cost};}
vector<ed>edges(){int m=pos.size();vector<ed>res(m);for(int i=0;i<m;i++)res[i]=get_edge(i);return res;}
pair<Cap,Cost>flow(int s,int t){return flow(s,t,numeric_limits<Cap>::max());}
pair<Cap,Cost>flow(int s,int t,Cap fl) {return slope(s,t,fl).back();}
vector<pair<Cap,Cost>>slope(int s,int t){return slope(s,t,numeric_limits<Cap>::max());}
vector<pair<Cap,Cost>>slope(int s,int t,Cap fl){vector<Cost>du(_n,0),di(_n);vector<int>pv(_n),pe(_n);vector<bool>vis(_n);auto du_ref=[&](){fill(di.begin(),di.end(),numeric_limits<Cost>::max());fill(pv.begin(),pv.end(),-1);fill(pe.begin(),pe.end(),-1);fill(vis.begin(),vis.end(),false);struct Q{Cost key;int to;bool operator<(Q r)const{return key>r.key;}};priority_queue<Q>qu;di[s]=0;qu.push({0,s});while(!qu.empty()){int v=qu.top().to;qu.pop();if(vis[v])continue;vis[v]=true;if(v==t)break;for(int i=0;i<g[v].size();i++){auto e=g[v][i];if(vis[e.to]||!e.cap)continue;Cost cost=e.cost-du[e.to]+du[v];if(di[e.to]-di[v]>cost){di[e.to]=di[v]+cost;pv[e.to]=v;pe[e.to]=i;qu.push(Q{di[e.to],e.to});}}}if(!vis[t])return false;for(int v=0;v<_n;v++){if(!vis[v])continue;du[v]-=di[t]-di[v];}return true;};Cap flow=0;Cost cost=0,pcpf=-1;vector<pair<Cap,Cost>>res;res.pb({flow,cost});while(flow<fl){if(!du_ref())break;Cap c=fl-flow;for(int v=t;v!=s;v=pv[v])chkmin(c,g[pv[v]][pe[v]].cap);for(int v=t;v!=s;v=pv[v]){auto& e=g[pv[v]][pe[v]];e.cap-=c;g[v][e.rev].cap+=c;}Cost d=-du[s];flow+=c;cost+=c*d;if(pcpf==d)res.pop_back();res.pb({flow,cost});pcpf=d;}return res;}
private:
int _n;
struct _ed{int to,rev;Cap cap;Cost cost;};
vector<pii>pos;
vector<vector<_ed>>g;
};
};
using namespace exfunc;
int main(){
}
/*
*/
|