#include #define REP(i,n) for(int i=0,i##_len=int(n);i=0;--i) #define rep(i,a,b) for(int i=int(a);i=int(b);--i) #define All(x) (x).begin(),(x).end() #define rAll(x) (x).rbegin(),(x).rend() #define ITR(i,x) for(auto i=(x).begin();i!=(x).end();++i) using std::cin; using std::cout; using std::cerr; using std::endl; using std::vector; using std::string; typedef long long ll; typedef std::pair P; typedef vector vi; typedef vector vvi; typedef vector vl; typedef vector vvl; constexpr ll mod = 1e9+7; constexpr double eps = 1e-9; const double PI = acos(-1); void solve(); namespace Math{ template bool chmax(T &a,T b){ if(a bool chmin(T &a,T b){ if(b check){ while(llabs(ok-ng)>1){ ll mid=ng-((ng-ok)>>1); if(check(mid)) ok=mid; else ng=mid; } return ok; } ll sqrt(ll n){ll s=1; while(s*s>n||n>=(s+1)*(s+1)){ s=(n/s+s)/2; } return s;} ll roundup(ll n,ll div){ if(n>0) return (n-1)/div+1; else return n/div; } bool square(ll a){ll n=Math::sqrt(a);return a==n*n;} template T pow(T x, ll n){ T ans = T(1); while(n != 0){ if(n&1) ans = ans*x; x = x*x; n = n >> 1; } return ans; } int digitsum(ll N,int a){ if(N==0) return 0; int ret=0; ret+=digitsum(N/a,a)+N%a; return ret; } ll gcd(ll x,ll y){return y ? gcd(y,x%y) : x;}; ll lcm(ll x,ll y){return x/Math::gcd(x,y)*y;} ll manhattan(const P &a,const P &b){return llabs(a.first-b.first)+llabs(a.second-b.second);} } namespace Solution{ using Graph = vector>; ll knapsack(int kinds,int MAX_W,const vl &weight,const vl &cost){ vector> dp(kinds+1,vector(MAX_W+1,0)); REP(i,kinds) REP(j,MAX_W+1){ if(j dp(MAX_W+1); REP(i,kinds) for(int j=weight[i];j<=MAX_W;++j){ dp[j] = std::max(dp[j],dp[j-weight[i]]+cost[i]); } return dp[MAX_W]; } ll huge_knapsack(int kinds,ll MAX_W,const vl &weight,const vl &cost){ int n2=kinds/2; vector

ps(1<<(kinds/2)); REP(i,1<>j&1){ sw += weight[j]; sv += cost[j]; } } ps[i] = P(sw,sv); } std::sort(ps.begin(),ps.begin() + (1<>j)&1){ sw += weight[n2+j]; sv += cost[n2+j]; } } if(sw<=MAX_W){ ll tv = (lower_bound(ps.begin(),ps.begin()+m,P(MAX_W-sw,LLONG_MAX))-1)->second; Math::chmax(res,sv+tv); } } return res; } ll choose_MonN(int N,int M,const vl &cost){ vvl dp(N+1,vl(M+1,-1LL<<58)); REP(i,N+1) dp[i][0]=0; REP(i,N) REP(j,M){ if(j>i) break; dp[i+1][j+1]=std::max(dp[i][j+1],dp[i][j]+cost[i]); } return dp[N][M]; } ll Partition_Func(int n,int k){ vector> dp(k+1,vector (n+1,0)); dp[0][0]=1; rep(i,1,k+1) REP(j,n+1){ if(j-i>=0) dp[i][j]=(dp[i][j-i]+dp[i-1][j]); else dp[i][j]=dp[i-1][j]; } return dp[k][n]; } int LCS(string s,string t){ int n=s.length(),m=s.length(); vector> dp(n+1,vector(m+1)); REP(i,n) REP(j,m){ if (s[i] == t[j]) dp[i+1][j+1] = dp[i][j] + 1; else dp[i+1][j+1] = std::max(dp[i][j+1], dp[i+1][j]); } return dp[n][m]; } int LIS(const vector a){ int n=a.size(); ll INF=1LL<<50; vector dp(n+1,INF); REP(i,n) *std::lower_bound(All(dp),a[i])=a[i]; int k=std::lower_bound(All(dp),INF)-dp.begin(); return k; } int max_flow(int s,int t,const vector> &g){ struct edge_{int to,cap, rev;}; vector used(g.size(),false); vector> G(g.size()); REP(i,g.size()) REP(j,g[i].size()){ int from = i, to = g[i][j].second; int cap = g[i][j].first; G[from].push_back((edge_){to,cap,(int)G[to].size()}); G[to].push_back((edge_){from,0,(int)G[from].size()-1}); } auto dfs = [&](auto&& f,int v,int t,int fl)->int{ if(v==t) return fl; used[v] = true; REP(i,G[v].size()){ edge_ &e = G[v][i]; if(!used[e.to] && e.cap>0){ int d = f(f, e.to,t,std::min(fl,e.cap)); if(d>0){ e.cap -= d; G[e.to][e.rev].cap += d; return d; } } } return 0; }; int flow=0; while(1){ used.assign(used.size(),false); int f = dfs(dfs,s,t,INT_MAX); if(f==0) return flow; flow += f; } } int bipartite_matching_calculate(const Graph &g){ int V = g.size(); vi match(V,-1); vector used(V,false); auto dfs = [&](auto&& f,int v)->bool{ used[v]=true; REP(i,g[v].size()){ int u = g[v][i], w = match[u]; if(w<0 || (!used[w] && f(f,w))){ match[v]=u; match[u]=v; return true; } } return false; }; int res=0; REP(v,V){ if(match[v] < 0){ used.assign(V,false); if(dfs(dfs,v)) res++; } } return res; } int bipartite_matching(int l,int r,std::function F){ int V = l+r; Graph G(V); REP(i,l) REP(j,r) if(F(i,j)){ G[i].push_back(l+j); G[l+j].push_back(i); } return bipartite_matching_calculate(G); } ll dinic(int s,int t,const vector> &graph){ struct max_flow { struct edge_ { int to;ll cap;int rev; }; int V; vector> G; vector itr, level; max_flow(int V) : V(V) { G.assign(V,vector()); } void add_edge(int from, int to, ll cap) { G[from].push_back((edge_) {to, cap, (int) G[to].size()}); G[to].push_back((edge_) {from, 0, (int) G[from].size()-1}); } void bfs(int s) { level.assign(V,-1); std::queue q; level[s] = 0; q.push(s); while (!q.empty()) { int v = q.front(); q.pop(); for(auto &e: G[v]) if(e.cap > 0 and level[e.to] < 0) { level[e.to] = level[v] + 1; q.push(e.to); } } } ll dfs(int v, int t, ll f) { if(v == t) return f; for(int& i=itr[v]; i < (int) G[v].size(); ++i) { edge_& e = G[v][i]; if (e.cap > 0 and level[v] < level[e.to]) { int d = dfs(e.to, t, std::min(f, e.cap)); if (d > 0) { e.cap -= d; G[e.to][e.rev].cap += d; return d; } } } return 0; } int run(int s, int t) { int ret = 0, f; while(bfs(s), level[t] >= 0) { itr.assign(V,0); while ((f = dfs(s, t, 1LL<<59)) > 0) ret += f; } return ret; } }; max_flow d(graph.size()); REP(i,graph.size()) for(auto e:graph[i]){ d.add_edge(i,e.second,e.first); } return d.run(s,t); } } namespace NTT { using uint = uint_fast32_t; // NTT_PRIMES {{{ constexpr ll NTT_PRIMES[][2] = { {1224736769, 3}, // 2^24 * 73 + 1, {1053818881, 7}, // 2^20 * 3 * 5 * 67 + 1 {1051721729, 6}, // 2^20 * 17 * 59 + 1 {1045430273, 3}, // 2^20 * 997 + 1 {1012924417, 5}, // 2^21 * 3 * 7 * 23 + 1 {1007681537, 3}, // 2^20 * 31^2 + 1 {1004535809, 3}, // 2^21 * 479 + 1 {998244353, 3}, // 2^23 * 7 * 17 + 1 {985661441, 3}, // 2^22 * 5 * 47 + 1 {976224257, 3}, // 2^20 * 7^2 * 19 + 1 {975175681, 17}, // 2^21 * 3 * 5 * 31 + 1 {962592769, 7}, // 2^21 * 3^3 * 17 + 1 {950009857, 7}, // 2^21 * 4 * 151 + 1 {943718401, 7}, // 2^22 * 3^2 * 5^2 + 1 {935329793, 3}, // 2^22 * 223 + 1 {924844033, 5}, // 2^21 * 3^2 * 7^2 + 1 {469762049, 3}, // 2^26 * 7 + 1 {167772161, 3}, // 2^25 * 5 + 1 }; ll extgcd(ll a, ll b, ll &x, ll &y) { ll d; return b == 0 ? (x = a < 0 ? -1 : 1, y = 0, a < 0 ? -a : a) : (d = extgcd(b, a % b, y, x), y -= a / b * x, d); } ll modinv(ll a, ll mod) { ll x, y; extgcd(a, mod, x, y); x %= mod; return x < 0 ? x + mod : x; } ll modpow(ll a, ll b, ll mod) { ll r = 1; a %= mod; while(b) { if(b & 1) r = r * a % mod; a = a * a % mod; b >>= 1; } return r; } // NTT Core {{{ template < int MAX_H > struct Pool { static ll *tmp, *A, *B; }; template < int MAX_H > ll *Pool< MAX_H >::tmp = new ll[1 << MAX_H]; template < int MAX_H > ll *Pool< MAX_H >::A = new ll[1 << MAX_H]; template < int MAX_H > ll *Pool< MAX_H >::B = new ll[1 << MAX_H]; template < int MAX_H, ll mod, ll primitive > class Core { public: static_assert((mod & ((1 << MAX_H) - 1)) == 1, "mod is too small; comment out"); // ord zetaList[i] = 2^(i + 1) ll zetaList[MAX_H], zetaInvList[MAX_H]; // constexpr Core() { zetaList[MAX_H - 1] = modpow(primitive, (mod - 1) / (1 << MAX_H), mod); zetaInvList[MAX_H - 1] = modinv(zetaList[MAX_H - 1], mod); for(int ih = MAX_H - 2; ih >= 0; --ih) { zetaList[ih] = zetaList[ih + 1] * zetaList[ih + 1] % mod; zetaInvList[ih] = zetaInvList[ih + 1] * zetaInvList[ih + 1] % mod; } } void fft(ll *a, uint n, uint nh, bool inverse) const { ll *tmp = Pool< MAX_H >::tmp; uint mask = n - 1; for(uint i = n >> 1, ih = nh - 1; i >= 1; i >>= 1, --ih) { ll zeta = inverse ? zetaInvList[nh - 1 - ih] : zetaList[nh - 1 - ih]; ll powZeta = 1; for(uint j = 0; j < n; j += i) { for(uint k = 0; k < i; ++k) { tmp[j | k] = (a[((j << 1) & mask) | k] + powZeta * a[(((j << 1) | i) & mask) | k]) % mod; } powZeta = powZeta * zeta % mod; } std::swap(a, tmp); } if(nh & 1) { std::swap(a, tmp); for(uint i = 0; i < n; ++i) a[i] = tmp[i]; } if(inverse) { ll invN = modinv(n, mod); for(uint i = 0; i < n; ++i) a[i] = a[i] * invN % mod; } } vector< ll > conv(const vector< ll > &a, const vector< ll > &b) const { uint t = a.size() + b.size() - 1; uint n = 1, nh = 0; while(n < t) n <<= 1, ++nh; return convStrict(a, b, n, nh); } vector< ll > convStrict(const vector< ll > &a, const vector< ll > &b, uint n, uint nh) const { ll *A = Pool< MAX_H >::A, *B = Pool< MAX_H >::B; for(uint i = 0; i < n; ++i) A[i] = B[i] = 0; copy(a.begin(), a.end(), A); copy(b.begin(), b.end(), B); fft(A, n, nh, 0), fft(B, n, nh, 0); for(uint i = 0; i < n; ++i) A[i] = A[i] * B[i] % mod; fft(A, n, nh, 1); return vector< ll >(A, A + n); } }; // Convolution With Garner {{{ template < int MAX_H, int I > class ConvolutionWithGarnerCore { public: static void conv_for(uint n, uint nh, const vector< ll > &a, const vector< ll > &b, vector< ll > &mods, vector< ll > &coeffs, vector< vector< ll > > &constants) { static const Core< MAX_H, NTT_PRIMES[I][0], NTT_PRIMES[I][1] > ntt; auto c = ntt.convStrict(a, b, n, nh); mods[I] = NTT_PRIMES[I][0]; ConvolutionWithGarnerCore< MAX_H, I - 1 >::conv_for( n, nh, a, b, mods, coeffs, constants); // garner for(size_t i = 0; i < c.size(); ++i) { ll v = (c[i] - constants[I][i]) * modinv(coeffs[I], mods[I]) % mods[I]; if(v < 0) v += mods[I]; for(size_t j = I + 1; j < mods.size(); ++j) { constants[j][i] = (constants[j][i] + coeffs[j] * v) % mods[j]; } } for(size_t j = I + 1; j < mods.size(); ++j) { coeffs[j] = (coeffs[j] * mods[I]) % mods[j]; } } }; template < int MAX_H > class ConvolutionWithGarnerCore< MAX_H, -1 > { public: static void conv_for(uint, uint, const vector< ll > &, const vector< ll > &, vector< ll > &, vector< ll > &, vector< vector< ll > > &) {} }; template < int MAX_H > class ConvolutionWithGarner { public: template < int USE > static vector< ll > conv(const vector< ll > &a, const vector< ll > &b, ll mod) { static_assert(USE >= 1, "USE must be positive"); static_assert(USE <= sizeof(NTT_PRIMES) / sizeof(*NTT_PRIMES), "USE is too big"); uint nt = a.size() + b.size() - 1; uint n = 1, nh = 0; while(n < nt) n <<= 1, ++nh; vector< ll > coeffs(USE + 1, 1); vector< vector< ll > > constants(USE + 1, vector< ll >(n)); vector< ll > mods(USE + 1, mod); ConvolutionWithGarnerCore< MAX_H, USE - 1 >::conv_for( n, nh, a, b, mods, coeffs, constants); return constants.back(); } }; } // 1st param is MAX_H NTT::Core< 18, NTT::NTT_PRIMES[0][0], NTT::NTT_PRIMES[0][1] > nttBig; NTT::Core< 18, 998244353, 5 > ntt; using nttconv = NTT::ConvolutionWithGarner< 18 >; // nttconv::conv< USE >(a, b, mod) template std::pair, double> min_ball(iter left, iter right, int seed = 1333) { const int n = right - left; using P=std::complex; using ld=double; assert(n >= 1); if (n == 1) { return {*left, ld(0)}; } std::mt19937 mt(seed); std::shuffle(left, right, mt); // std::random_shuffle(left, right); // simple but deprecated iter ps = left; using circle = std::pair; auto cross=[](const P& a, const P& b) { return a.real()*b.imag() - a.imag()*b.real(); }; auto dot=[](const P& a, const P& b) { return a.real()*b.real() + a.imag()*b.imag(); }; auto make_circle_3 = [&](const P &a, const P &b, const P &c) -> circle { ld A = std::norm(b - c), B = std::norm(c - a), C = std::norm(a - b), S = cross(b - a, c - a); P p = (A * (B + C - A) * a + B * (C + A - B) * b + C * (A + B - C) * c) / (4 * S * S); ld r2 = std::norm(p - a); return {p, r2}; }; auto make_circle_2 = [](const P &a, const P &b) -> circle { P c = (a + b) / (ld)2; ld r2 = std::norm(a - c); return {c, r2}; }; auto in_circle = [](const P &a, const circle &c) -> bool { const double eps=1e-9; return std::norm(a - c.first) <= c.second + eps; }; circle c = make_circle_2(ps[0], ps[1]); // MiniDisc for (int i = 2; i < n; ++i) { if (!in_circle(ps[i], c)) { // MiniDiscWithPoint c = make_circle_2(ps[0], ps[i]); for (int j = 1; j < i; ++j) { if (!in_circle(ps[j], c)) { // MiniDiscWith2Points c = make_circle_2(ps[i], ps[j]); for (int k = 0; k < j; ++k) { if (!in_circle(ps[k], c)) { c = make_circle_3(ps[i], ps[j], ps[k]); } } } } } } return c; } template size_t KMP(const V &Search,const V &Word){ size_t i=2,j=0; std::vector Table(Search.size()); Table[0]=-1;Table[1]=0; while(i0) j=Table[j]; else{ Table[i]=0; ++i; } } i=0;j=0; while(j+i0) i=Table[i]; } } return Search.size(); } template vector Zalgo(V s){ vector A(s.size(),0); A[0]=s.size(); int i = 1, j = 0; while (i < s.size()) { while (i+j < s.size() && s[j] == s[i+j]) ++j; A[i] = j; if (j == 0) { ++i; continue;} int k = 1; while (i+k < s.size() && k+A[k] < j) A[i+k] = A[k], ++k; i += k; j -= k; } return A; } template struct rolling_hash { int n; const long long MOD[2] = {999999937LL, 1000000007LL}, base = 9973; vector hs[2], pw[2]; rolling_hash(){} rolling_hash(const V &s) { n = s.size(); for (int i = 0; i < 2; i++) { hs[i].assign(n+1,0); pw[i].assign(n+1,0); hs[i][0] = 0; pw[i][0] = 1; for (int j = 0; j < n; j++) { pw[i][j+1] = pw[i][j]*base%MOD[i]; hs[i][j+1] = (hs[i][j]*base+s[j])%MOD[i]; } } } long long hash(int l, int r, int i) { return ((hs[i][r]-hs[i][l]*pw[i][r-l])%MOD[i]+MOD[i])%MOD[i]; } bool match(int l1, int r1, int l2, int r2) { bool ret = 1; for (int i = 0; i < 2; i++) ret &= hash(l1,r1,i)==hash(l2,r2,i); return ret; } bool match(int l, int r, long long h[]) { bool ret = 1; for (int i = 0; i < 2; i++) ret &= hash(l,r,i)==h[i]; return ret; } }; template struct RH{ using u64 = std::uint_fast64_t; const u64 MASK30 = (1<<30)-1,MASK31 = (1LL<<31) -1, MOD = (1LL<<61)-1; const u64 Posint = MOD*3; const u64 base = 9973; std::vector Hash,POW; RH(V s){ int n=s.size(); Hash.resize(n+1); POW.resize(n+1,1); for(size_t i=0;i>31; u64 ad = a&MASK31; u64 bu = b>>31; u64 bd = b&MASK31; u64 mid = ad*bu + au*bd; u64 midu = mid>>30; u64 midd = mid&MASK30; return au*bu*2 + midu + (midd<<31) + ad*bd; } u64 CalcMod(u64 val){ val = (val & MOD) + (val>>61); if(val>=MOD) val-=MOD; return val; } }; class mint { private: ll _num,_mod=mod; mint set(ll num){ _num = num ; if(_num<0){ if(_num>=-_mod)_num=_mod+_num; else _num=_mod-(-_num)%_mod; } else if(_num>=_mod) _num%=_mod; return *this; } ll imod()const{ ll n=_mod-2; ll ans = 1,x=_num; while(n != 0){ if(n&1) ans = ans*x%_mod; x = x*x%_mod; n = n >> 1; } return ans; } public: explicit mint(){ _num = 0; } explicit mint(ll num){ _num = num; if(_num<0){ if(_num>=-_mod)_num=_mod+_num; else _num=_mod-(-_num)%_mod; } else if(_num>=_mod) _num%=_mod; } explicit mint(ll num,ll M){ _mod=M; _num=num; if(_num<0){ if(_num>=-_mod)_num=_mod+_num; else _num=_mod-(-_num)%_mod; } else if(_num>=_mod) _num%=_mod; } mint(const mint &cp){_num=cp._num;_mod=cp._mod;} mint operator+ (const mint &x)const{ return mint(_num + x._num , _mod); } mint operator- (const mint &x)const{ return mint(_num - x._num , _mod);} mint operator* (const mint &x)const{ return mint(_num * x._num , _mod); } mint operator/ (const mint &x)const{ return mint(_num * x.imod() , _mod);} mint operator+=(const mint &x){ return set(_num + x._num); } mint operator-=(const mint &x){ return set(_num - x._num); } mint operator*=(const mint &x){ return set(_num * x._num); } mint operator/=(const mint &x){ return set(_num * x.imod());} mint operator= (const ll x){ return set(x); } mint operator+ (const ll x)const{return *this + mint(x,_mod); } mint operator- (const ll x)const{ return *this - mint(x,_mod); } mint operator* (const ll x)const{ return *this * mint(x,_mod); } mint operator/ (const ll x)const{ return *this/mint(x);} mint operator+=(const ll x){ *this = *this + x;return *this; } mint operator-=(const ll x){ *this = *this - x;return *this; } mint operator*=(const ll x){ *this = *this * x;return *this;} mint operator/=(const ll x){ *this = *this / x;return *this;} bool operator==(const mint &x)const{return _num==x._num;} bool operator!=(const mint &x)const{return _num!=x._num;} friend mint operator+(ll x,const mint &m){return mint(m._num + x , m._mod);} friend mint operator-(ll x,const mint &m){return mint( x - m._num , m._mod);} friend mint operator*(ll x,const mint &m){return mint(m._num * (x % m._mod) , m._mod);} friend mint operator/(ll x,const mint &m){return mint(m.imod() * (x % m._mod) , m._mod);} explicit operator ll() { return _num; } explicit operator int() { return (int)_num; } friend std::ostream& operator<<(std::ostream &os, const mint &x){ os << x._num; return os; } friend std::istream& operator>>(std::istream &is, mint &x){ll val; is>>val; x.set(val); return is;} }; ll inv_mod(ll a,ll _mod){return (ll)Math::pow(mint(a,_mod),_mod-2);} class Factorial{ private: vector fac; vector ifac; public: Factorial(ll N){ fac.reserve(N+1); fac.push_back(1); REP(i,N) fac.push_back(fac[i]*(i+1)%mod); ifac.resize(N+1); ifac[N]=inv_mod(fac[N],mod); REP(i,N) ifac[N-1-i]=(ifac[N-i]*(N-i))%mod; } ll fact(ll a){return fac[a];} ll ifact(ll a){return ifac[a];} ll cmb(ll a,ll b){ if(a==0&&b==0) return 1; if(a(const rational &a, const rational &b){return b < a;} friend bool operator>=(const rational &a, const rational &b){return !(a < b);} friend bool operator==(const rational &a, const rational &b){return !(a < b) && !(b < a);} friend bool operator!=(const rational &a, const rational &b){return (a < b) || (b < a);} friend std::ostream& operator<<(std::ostream &os, const rational &x){ printf("%.16f",(double)x.p/(double)x.q); return os; } friend std::istream& operator>>(std::istream &is, rational &x){is>>x.p>>x.q; x.normalize(); return is;} }; template class MAT{ private: int row,col; vector> _A; double eps = 1e-9; MAT set(vector> A){_A = A ; return *this;} public: MAT(){ } MAT(int n,int m=0,T x=T(0)){ if(n<1 || m<0){cout << "err Matrix::Matrix" < a(col,x); _A.push_back(a); } if(m==0) REP(i,n) _A[i][i]=1.0; } MAT(vector> A){row=A.size();col=A[0].size();_A=A;} MAT(const MAT &cp){_A=cp._A;row=cp.row;col=cp.col;} T* operator[] (int i){return _A[i].data();} MAT operator= (vector> x) {return set(x);} MAT operator+ (MAT x) const { if(row!=x.row || col!=x.col){ cerr << "err Matrix::operator+" <>(){return _A;} friend std::ostream &operator<<(std::ostream &os,const MAT &x){ REP(i,x.row) REP(j,x.col) os<>(std::istream &is,MAT &x){REP(i,x.row) REP(j,x.col) is>>x._A[i][j];return is;} size_t size_row()const{return row;} size_t size_col()const{return col;} MAT transpose()const{ MAT r(col,row); REP(i,col) REP(j,row) r[i][j]=_A[j][i]; return r; } MAT inverse()const{ T buf; MAT inv_a(row,0); vector> a=_A; //row reduction REP(i,row){ buf=1/a[i][i]; REP(j,row){ a[i][j]*=buf; inv_a[i][j]*=buf; } REP(j,row){ if(i!=j){ buf=a[j][i]; REP(k,row){ a[j][k]-=a[i][k]*buf; inv_a[j][k]-=inv_a[i][k]*buf; } } } } return inv_a; } MAT Jacobi(MAT b)const{//ヤコビ法によって解を求める size_t sz=row; MAT D(sz,sz),inD(sz,sz),H(sz,sz),N(sz,sz); MAT c(sz,1),x(sz,1),tmp(sz,1); //cout<<"initialized"< DL(row),U(row),inDL(row),H(row),c(row,1),x(row,1),tmp(row,1); for(int i=0;i=j){ DL[i][j] = _A[i][j]; U[i][j] = 0; } else{ DL[i][j] = 0; U[i][j] = _A[i][j]; } } x[i][0] = 1; } inDL=DL.inverse(); c=inDL*b; H=inDL*U; int n=0; while(1){ tmp=x; x=c-H*x; T r = T(0); for(int i=0;i> A=_A; const int n = row, m = col; int r = 0; for(int i = 0; r < n && i < m; ++i) { int pivot = r; for(int j = r+1; j < n; ++j) if(fabs(A[j][i]) > fabs(A[pivot][i])) pivot = j; swap(A[pivot], A[r]); if(fabs(A[r][i]) < eps) continue; for (int k = m-1; k >= i; --k) A[r][k] /= A[r][i]; rep(j,r+1,n) rep(k,i,m) A[j][k] -= A[r][k] * A[j][i]; ++r; } return r; } }; class UnionFind{ private: vector Parent,es; vector diff_weight; public: UnionFind(int N){ es.resize(N,0); Parent.resize(N,-1); diff_weight.resize(N,0LL); } int root(int A){ if(Parent[A]<0) return A; else{ int r = root(Parent[A]); diff_weight[A] += diff_weight[Parent[A]]; // 累積和をとる return Parent[A]=r; } } bool issame(int A,int B){ return root(A)==root(B); } ll weight(int x) { root(x); // 経路圧縮 return diff_weight[x]; } ll diff(int x, int y) { return weight(y) - weight(x); } int size(int A){ return -Parent[root(A)]; } int eize(int A){ return es[root(A)]; } bool connect(int A,int B){ A=root(A); B=root(B); if(A==B) return false; if(size(A)(const edge& y)const{ if(cost!=y.cost) return cost>y.cost; else if(to!=y.to) return to>y.to; else return from>y.from;} bool operator==(const edge& y) const{return !(*thisy);} }; class lca { public: using Graph = vector>; const int n = 0; const int log2_n = 0; std::vector> parent; std::vector depth; lca() {} lca(const Graph &g, int root) : n(g.size()), log2_n(log2(n) + 1), parent(log2_n, std::vector(n)), depth(n) { dfs(g, root, -1, 0); for (int k = 0; k + 1 < log2_n; k++) { for (int v = 0; v < (int)g.size(); v++) { if (parent[k][v] < 0) parent[k + 1][v] = -1; else parent[k + 1][v] = parent[k][parent[k][v]]; } } } void dfs(const Graph &g, int v, int p, int d) { parent[0][v] = p; depth[v] = d; REP(j,g[v].size()) { if (g[v][j] != p) dfs(g, g[v][j], v, d + 1); } } int get(int u, int v) { if (depth[u] > depth[v]) std::swap(u, v); for (int k = 0; k < log2_n; k++) { if ((depth[v] - depth[u]) >> k & 1) { v = parent[k][v]; } } if (u == v) return u; for (int k = log2_n - 1; k >= 0; k--) { if (parent[k][u] != parent[k][v]) { u = parent[k][u]; v = parent[k][v]; } } return parent[0][u]; } }; class SOSU{ private: vector Prime_Number; vector isp; public: SOSU(int N){ isp.resize(N+1,true); isp[0]=isp[1]=false; rep(i,2,N+1) if(isp[i]){ Prime_Number.push_back(i); for(int j=2*i;j<=N;j+=i) isp[j]=false; } } int operator[](int i){return Prime_Number[i];} int size(){return Prime_Number.size();} int back(){return Prime_Number.back();} bool isPrime(int q){return isp[q];} }; class Divisor{//素因数分解をしてから約数列挙、分解結果はP(底,指数)でpfacにまとめている private: vector F; vector

pfactorize; public: Divisor(ll N){ for(ll i = 1; i * i <= N; i++) { if(N % i == 0) { F.push_back(i); if(i * i != N) F.push_back(N / i); } } sort(begin(F), end(F)); SOSU p(Math::sqrt(N)+1); REP(i,p.size()){ pfactorize.push_back(P(p[i],0)); while(N%p[i]==0){ N/=p[i]; pfactorize.back().second++; } if(pfactorize.back().second==0) pfactorize.pop_back(); } if(N>1) pfactorize.push_back(P(N,1)); } int size(){return F.size();} vector

pfac(){return pfactorize;} ll operator[](int k){return F[k];} }; template struct compress{ std::map zip; vector unzip; compress(vector x){ unzip=x; sort(All(x)); x.erase(unique(All(x)),x.end()); REP(i,x.size()){ zip[x[i]]=i; } } size_t operator[](int k){return zip[unzip[k]];} }; template class SegmentTree{ private: typedef std::function F; int n; T d0; vector vertex; F f; F g; public: SegmentTree(F f,F g,T d):d0(d),f(f),g(g){} void init(int _n){ n=1; while(n<_n) n*=2; vertex.resize(2*n-1,d0); } void build(const vector &v){ int n_=v.size(); init(n_); for(int i=0;i=0;i--) vertex[i]=f(vertex[2*i+1],vertex[2*i+2]); } void update(int i,T x){ int k=i+n-1; vertex[k]=g(vertex[k],x); while(k>0){ k=(k-1)/2; vertex[k]=f(vertex[2*k+1],vertex[2*k+2]); } return; } T query(int l,int r){ T vl=d0,vr=d0; l+=n-1; r+=n-1; for(;l<=r;l/=2,r=r/2-1){ if(l%2==0) vl=f(vl,vertex[l]); if(r&1) vr=f(vr,vertex[r]); } return f(vl,vr); } }; template struct LazySegmentTree{ using F = std::function; using G = std::function; using H = std::function; int n,height; F f; G g; H h; T ti; E ei; vector dat; vector laz; LazySegmentTree(F f,G g,H h,T ti,E ei): f(f),g(g),h(h),ti(ti),ei(ei){} void init(int n_){ n=1;height=0; while(n &v){ int n_=v.size(); init(n_); for(int i=0;i>i); } inline void recalc(int k){ while(k>>=1) dat[k]=f(reflect((k<<1)|0),reflect((k<<1)|1)); } void update(int a,int b,E x){ thrust(a+=n); thrust(b+=n-1); for(int l=a,r=b+1;l>=1,r>>=1){ if(l&1) laz[l]=h(laz[l],x),l++; if(r&1) --r,laz[r]=h(laz[r],x); } recalc(a); recalc(b); } void set_val(int a,T x){ thrust(a+=n); dat[a]=x;laz[a]=ei; recalc(a); } T query(int a,int b){ thrust(a+=n); thrust(b+=n-1); T vl=ti,vr=ti; for(int l=a,r=b+1;l>=1,r>>=1) { if(l&1) vl=f(vl,reflect(l++)); if(r&1) vr=f(reflect(--r),vr); } return f(vl,vr); } template int find(int st,C &check,T &acc,int k,int l,int r){ if(l+1==r){ acc=f(acc,reflect(k)); return check(acc)?k-n:-1; } propagate(k); int m=(l+r)>>1; if(m<=st) return find(st,check,acc,(k<<1)|1,m,r); if(st<=l&&!check(f(acc,dat[k]))){ acc=f(acc,dat[k]); return -1; } int vl=find(st,check,acc,(k<<1)|0,l,m); if(~vl) return vl; return find(st,check,acc,(k<<1)|1,m,r); } template int find(int st,C &check){ T acc=ti; return find(st,check,acc,1,0,n); } }; template class FFT{ private: using Complex = std::complex; std::vector C; void DFT(std::vector &F,size_t n,int sig=1){ if(n==1) return; std::vector f0(n/2),f1(n/2); for(size_t i=0;i &f,size_t n){ DFT(f,n,-1); for(size_t i=0;i &A,const std::vector &B){ size_t n=1; while(n<=A.size()+B.size()){ n*=2; } std::vector g(n,Complex(0)); C.resize(n,Complex(0)); size_t i_len=std::max(A.size(),B.size()); for(size_t i=0;i>N>>M>>K>>p>>q; vector b(N); REP(i,N) cin>>b[i]; mint ans(0); mint z=(Math::pow((q-2*p)/q,K)+1)/mint(2); REP(i,N){ if(i