#include"bits/stdc++.h" #include"atcoder/all" using namespace atcoder; #pragma GCC target("avx2,fma") #pragma GCC optimize("Ofast") #pragma GCC optimize("unroll-loops") #pragma GCC target("sse,sse2,sse3,ssse3,sse4,popcnt,abm,mmx,avx,tune=native") using namespace std; using ll = long long; using vl = vector; using vvl = vector; using vvvl = vector; using vvvvl = vector; using vi = vector; using vvi = vector; using vvvi = vector; using vvvvi = vector; typedef pair P; typedef pair Pd; typedef tuple PP; typedef tuple PPP; using vs = vector; using vvs = vector; using vvvs = vector; using vp = vector

; using vvp = vector; using vvvp = vector; using vb = vector; using vvb = vector; using vvvb = vector; #define lb(v, k) (lower_bound(all(v), (k)) - v.begin()) #define ub(v, k) (upper_bound(all(v), (k)) - v.begin()) #define eb emplace_back #define fi first #define se second #define pq(T) priority_queue #define pqr(T) priority_queue, greater> #define pcount(i) __builtin_popcountll(i) #define rep(i, n) for (ll i = 0; i < (ll)(n); i++) #define repi(i, a, b) for (ll i = (ll)(a); i < (ll)(b); i++) #define all(a) (a).begin(), (a).end() #define rll(a) (a).rbegin(), (a).rend() #define double long double #define yesno(i) if(i)cout<<"yes"<> x.at(hfuaig); \ } #define coutvece(x) \ for (ll hfuaig = 0; hfuaig < x.size(); hfuaig++) \ { \ cout << x.at(hfuaig) << endl; \ } #define coutvec(x) \ for (ll hfuaig = 0; hfuaig < x.size(); hfuaig++) \ { \ cout << x.at(hfuaig) << ' '; \ }\ cout << endl; vl dx = {-1,1 ,0 ,0 ,-1, 1, -1, 1}; vl dy = {0 ,0 ,-1,1 ,-1, -1, 1, 1}; pair dxy(pair a, ll i){ return make_pair(a.fi + dx[i], a.se + dy[i]); } template bool chmax(T &a, const T& b){ if(a < b){ a = b; return true; } return false; } template bool chmin(T &a, const T& b){ if(a > b){ a = b; return true; } return false; } ll mod = 998244353; int modd = 998244353; const ll Mod = 1000000007; const ll inf = 999999999999999999LL; const int INF = 999999999; using mint = modint998244353; using mmint = modint1000000007; using intm = static_modint<1000000009>; using vm = vector; using vvm = vector; using vvvm = vector; using vvvvm = vector; using vvvvvm = vector; using vmm = vector; using vvmm = vector; using vvvmm = vector; ll mygcd(ll A, ll B){ if(A == -1)return B; if(B == -1)return A; if(B == 0 || A == 0){ return A ^ B; } if(A % B == 0){ return B; } else{ return mygcd(B, A % B); } } ll mylcm(ll A, ll B){ if(A == -1){ return B; } if(B == -1){ return A; } return A / mygcd(A, B) * B; } typedef struct Point_Coordinates { double x, y; } point; typedef struct Point_Coordinatesll { ll x, y; } pointl; long long modpow(long long a, long long n, long long mo) { long long res = 1; while (n > 0) { if (n & 1) res = res * a % mo; a = a * a % mo; n >>= 1; } return res % mo; } //円の位置関係 bool isd(double a, double b, double c){ if(c > (a + b) * (a + b)){ return false; } if(c == (a + b) * (a + b)){ return true; } if(abs(a - b) * abs(a - b) < c && c < (a + b) * (a + b)){ return true; } if(c == abs(a - b) * abs(a - b)){ return true; } if(c < abs(a - b) * abs(a - b)){ return false; } return true; } const int MAX = 201; int MOD1 = mod; ll fact[MAX], inv_fact[MAX], inv[MAX], perm[MAX]; void init() { MOD1 = mod; // 初期値設定と1はじまりインデックスに直す fact[0] = 1; fact[1] = 1; inv[0] = 1; inv[1] = 1; inv_fact[0] = 1; inv_fact[1] = 1; // メモの計算 repi(i, 2, MAX){ // 階乗 fact[i] = fact[i - 1] * i % MOD1; // 逆元 inv[i] = MOD1 - inv[MOD1%i] * (MOD1 / i) % MOD1; // 逆元の階乗 inv_fact[i] = inv_fact[i - 1] * inv[i] % MOD1; } perm[0] = 1; rep(i, MAX - 1){ perm[i + 1] = ((ll)(i + 1) * perm[i]) % MOD1; } } ll nck(int n, int k) { if (n < k) return 0; // 例外処理 if (n < 0 || k < 0) return 0; // 例外処理 if(k == 0)return 1; ll x = fact[n]; // n!の計算 ll y = inv_fact[n-k]; // (n-k)!の計算 ll z = inv_fact[k]; // k!の計算 return x * ((y * z) % MOD1) % MOD1; //二項係数の計算 } ll npk(ll n, ll k){ return (nck(n, k) * perm[k]) % MOD1; } ll nhk(ll n, ll k){ return nck(n + k - 1, k); } struct graph{ ll N; vvp G; vl dis; vl prev; graph(ll n) : N(n){ G.resize(n); } void push(ll a, ll b){ G[a].push_back({b, 1}); } void push(ll a, ll b, ll c){ G[a].push_back({b, c}); } vl dijkstra(ll i){ dis.assign(N, inf); prev.assign(N, -1); priority_queue, greater

> piq; // 「仮の最短距離, 頂点」が小さい順に並ぶ dis[i] = 0; piq.emplace(dis[i], i); while (!piq.empty()) { P p = piq.top(); piq.pop(); ll v = p.second; if (dis[v] < p.first) { // 最短距離で無ければ無視 continue; } for (auto &e : G[v]) { if (dis[e.fi] > dis[v] + e.se) { // 最短距離候補なら priority_queue に追加 dis[e.fi] = dis[v] + e.se; prev[e.fi] = v; piq.emplace(dis[e.fi], e.fi); } } } return dis; } vl get_path(ll t){ vl path; for (ll cur = t; cur != -1; cur = prev[cur]) { path.push_back(cur); } reverse(path.begin(), path.end()); // 逆順なのでひっくり返す return path; } }; struct BIT { private: vector bit; ll N; public: BIT(ll size) { N = size; bit.resize(N + 1, 0); } // 一点更新です void add(ll a, ll w) { for (int x = a; x <= N; x += x & -x) bit[x] += w; } // 1~Nまでの和を求める。 ll sum(ll a) { ll ret = 0; for (ll x = a; x > 0; x -= x & -x) ret += bit[x]; return ret; } }; //転倒数ライブラリ 1-indexed!! ll numfalls(vl &A){ ll ans = 0; ll N = A.size(); BIT b(N); rep(i, N){ ans += i - b.sum(A.at(i)); b.add(A.at(i), 1); } return ans; } //#define _GLIBCXX_DEBUG #define mint998 modint998244353 #define mint107 modint1000000007 //約数列挙 vector div(long long n) { vector ret; set re; ll N = sqrt(n) + 1; for (long long i = 1; i <= N; i++) { if (n % i == 0) { re.insert(n / i); re.insert(i); } } for(auto value :re){ ret.push_back(value); } return ret; } ll digit_sum(ll X){ ll ans = 0; while(X > 0){ ans += X%10; X/=10; } return ans; } void xypress(vl &A){ vl B = A; sort(all(B)); B.erase(unique(B.begin(), B.end()), B.end()); vector res(A.size()); for (int i = 0; i < A.size(); ++i) { res[i] = lower_bound(B.begin(), B.end(), A[i]) - B.begin(); } A = res; } //素数列挙 std::vector prime( const ll N ) { std::vector is_prime( N + 1 ); for( ll i = 0; i <= N; i++ ) { is_prime[ i ] = true; } std::vector P; for( ll i = 2; i <= N; i++ ) { if( is_prime[ i ] ) { for( ll j = 2 * i; j <= N; j += i ) { is_prime[ j ] = false; } P.push_back( i ); } } return P; } //mod m での逆元 long long modinv(long long a, long long m) { long long b = m, u = 1, v = 0; while (b) { long long t = a / b; a -= t * b; swap(a, b); u -= t * v; swap(u, v); } u %= m; if (u < 0) u += m; return u; } //mod 998244353での割り算 ll inve(ll a, ll b){ a %= mod; return (a * modinv(b, mod) % mod); } vl vecsum(vl x){ vl s = {0}; rep(i, x.size()){ s.push_back(s.back() + x.at(i)); } return s; } const int N_MAX = 10; ll spf[N_MAX]; // smallest prime factors void prepare_factorize() { rep(i, N_MAX) spf[i] = i; for (int p = 2; p * p <= N_MAX; p++) { for (int i = p; i < N_MAX; i += p) { if (spf[i] == i) spf[i] = p; } } } // 素因数分解 // その素因数が何個あるかのmapを返す map allprime(ll n) { map c; while (n != 1) { ll p = spf[n]; while (n % p == 0) { n /= p; c[p]++; } } return c; }/////////////////////////////////// vl prime2(ll n) { vl c(500); while (n != 1) { ll p = spf[n]; while (n % p == 0) { n /= p; c[p]++; } } return c; }/////////////////////////////////// vector primes(ll n){ set c; for(ll i = 2;i * i <= n;i++){ if(n % i == 0)c.insert(i); while(n % i == 0){ n /= i; } } if(n != 1){ c.insert(n); } vl d; for(auto a : c){ d.push_back(a); } return d; } vector SA_IS(vector str, ll var) { if(str.size() == 1) { vector ret(1,0); return ret; } str.push_back(0); ll si = str.size(); vector st(var, 0), en(var, 0); vector SL(si, 0); //s..0, l..1 vector SA(si, -1); vector LMS; vector is_LMS(si, -1); rep(i,str.size()) en[str[i]]++; for(ll i = 1; i < var; i++) en[i] += en[i-1]; for(ll i = 1; i < var; i++) st[i] = en[i-1]; SL[str.size()-1] = 0; for(ll i = str.size()-2; i >= 0; i--) { if(str[i] == str[i+1]) { SL[i] = SL[i+1]; continue; } if(str[i] > str[i+1]) SL[i] = 1; else SL[i] = 0; } for(ll i = 1; i < str.size(); i++) { if(SL[i] == 0 && SL[i-1] == 1) { SA[--en[str[i]]] = i; LMS.push_back(i); is_LMS[i] = 1; } } rep(i,var-1) en[i] = st[i+1]; en[var-1] = str.size(); rep(i,str.size()) if(SA[i] > 0 && SL[SA[i]-1] == 1) { SA[st[str[SA[i]-1]]++] = SA[i]-1; } st[0] = 0; for(ll i = 1; i < var; i++) st[i] = en[i-1]; for(ll i = 1; i < str.size(); i++) if(SA[i] != -1 && SL[SA[i]] == 0) { SA[i] = -1; } for(ll i = str.size()-1; i >= 1; i--) if(SA[i] > 0 && SL[SA[i]-1] == 0) { SA[--en[str[SA[i]-1]]] = SA[i]-1; } rep(i,var-1) en[i] = st[i+1]; en[var-1] = str.size(); ll counter = 0; vector pre_sa, new_sa; rep(i,SA.size()) if(is_LMS[SA[i]] != -1) { is_LMS[SA[i]] = ++counter; new_sa.clear(); for(ll j = SA[i]; j < SA.size(); j++) { new_sa.push_back(str[j]); if(j != SA[i] && is_LMS[j] != -1) { break; } } if(pre_sa == new_sa) { is_LMS[SA[i]] = --counter; } pre_sa = new_sa; } vector new_str; vector rev((ll)LMS.size()+1, 0); counter = 0; rep(i,is_LMS.size()) { if(is_LMS[i] != -1) { new_str.push_back(is_LMS[i]); rev[counter++] = i; } } vector rec = SA_IS(new_str, new_str.size()+1); rep(i,SA.size()) SA[i] = -1; for(ll i = rec.size()-1; i >= 0; i--) { SA[--en[str[rev[rec[i]]]]] = rev[rec[i]]; } rep(i,var-1) en[i] = st[i+1]; en[var-1] = str.size(); rep(i,str.size()) if(SA[i] > 0 && SL[SA[i]-1] == 1) { SA[st[str[SA[i]-1]]++] = SA[i]-1; } for(ll i = 1; i < str.size(); i++) if(SA[i] != -1 && SL[SA[i]] == 0) { SA[i] = -1; } for(ll i = str.size()-1; i >= 1; i--) if(SA[i] > 0 && SL[SA[i]-1] == 0) { SA[--en[str[SA[i]-1]]] = SA[i]-1; } SA.erase(SA.begin()); return SA; } int Judge(point &a, point &b, point &c, point &d) { double s, t; s = (a.x - b.x) * (c.y - a.y) - (a.y - b.y) * (c.x - a.x); t = (a.x - b.x) * (d.y - a.y) - (a.y - b.y) * (d.x - a.x); if (s * t > 0) return false; s = (c.x - d.x) * (a.y - c.y) - (c.y - d.y) * (a.x - c.x); t = (c.x - d.x) * (b.y - c.y) - (c.y - d.y) * (b.x - c.x); if (s * t > 0) return false; return true; } int Judgel(pointl &a, pointl &b, pointl &c, pointl &d) { ll s, t; s = (a.x - b.x) * (c.y - a.y) - (a.y - b.y) * (c.x - a.x); t = (a.x - b.x) * (d.y - a.y) - (a.y - b.y) * (d.x - a.x); if (s * t > 0) return false; s = (c.x - d.x) * (a.y - c.y) - (c.y - d.y) * (a.x - c.x); t = (c.x - d.x) * (b.y - c.y) - (c.y - d.y) * (b.x - c.x); if (s * t > 0) return false; return true; } double kyori(point &a, point &b){ double x = abs(a.x - b.x), y = abs(a.y - b.y); return sqrt(x * x + y * y); } double eps = 1; struct Edge { ll to; }; using Graph = vector>; vector topo_sort(const Graph &G) { // bfs vector ans; ll n = (ll)G.size(); vector ind(n); // ind[i]: 頂点iに入る辺の数(次数) for (ll i = 0; i < n; i++) { // 次数を数えておく for (auto e : G[i]) { ind[e.to]++; } } queue que; for (ll i = 0; i < n; i++) { // 次数が0の点をキューに入れる if (ind[i] == 0) { que.push(i); } } while (!que.empty()) { // 幅優先探索 ll now = que.front(); ans.push_back(now); que.pop(); for (auto e : G[now]) { ind[e.to]--; if (ind[e.to] == 0) { que.push(e.to); } } } return ans; } struct LCA { vector> parent; // parent[k][u]:= u の 2^k 先の親 vector dist; // root からの距離 LCA(const Graph &G, int root = 0) { init(G, root); } // 初期化 void init(const Graph &G, int root = 0) { int V = G.size(); int K = 1; while ((1 << K) < V) K++; parent.assign(K, vector(V, -1)); dist.assign(V, -1); dfs(G, root, -1, 0); for (int k = 0; k + 1 < K; k++) { for (int v = 0; v < V; v++) { if (parent[k][v] < 0) { parent[k + 1][v] = -1; } else { parent[k + 1][v] = parent[k][parent[k][v]]; } } } } // 根からの距離と1つ先の頂点を求める void dfs(const Graph &G, int v, int p, int d) { parent[0][v] = p; dist[v] = d; for (auto e : G[v]) { if (e.to != p) dfs(G, e.to, v, d + 1); } } int query(int u, int v) { if (dist[u] < dist[v]) swap(u, v); // u の方が深いとする int K = parent.size(); // LCA までの距離を同じにする for (int k = 0; k < K; k++) { if ((dist[u] - dist[v]) >> k & 1) { u = parent[k][u]; } } // 二分探索で LCA を求める if (u == v) return u; for (int k = K - 1; k >= 0; k--) { if (parent[k][u] != parent[k][v]) { u = parent[k][u]; v = parent[k][v]; } } return parent[0][u]; } int get_dist(int u, int v) { return dist[u] + dist[v] - 2 * dist[query(u, v)]; } bool is_on_path(int u, int v, int a) { return get_dist(u, a) + get_dist(a, v) == get_dist(u, v); } }; struct matrix{ ll h, w; vvl X; matrix(ll H, ll W) : h(H), w(W){ init(); } void init(){ X.assign(h, vl(w, 0)); } void set(ll H, ll W, ll c){ X.at(H).at(W) = c; } ll get(ll H, ll W){ return X.at(H).at(W); } }; struct rollinghash{ __int128_t m = (1LL << 61) - 1; __int128_t base; string S; int N; vl h, sh, pw; rollinghash(string x) : S(x), N((int)x.size()){ base = rand() % (m - 2) + 2; st(); } void st(){ h.assign(N + 1, 0); sh.assign(N + 1, 0); pw.assign(N + 1, 1); for(int i = 0;i < N;i++){ pw.at(i + 1) = (__int128_t)pw.at(i) * base % m; } for(int i = 0;i < N;i++){ h.at(i + 1) = ((__int128_t)h.at(i) * base + (__int128_t)S.at(i)) % m; } for(int i = N - 1;i >= 0;i--){ sh.at(i) = ((__int128_t)sh.at(i + 1) * base + (__int128_t)S.at(i)) % m; } } ll hash(int l, int r){ if(l == r){ return 0; } else if(l < r){ return (h.at(r) - (__int128_t)h.at(l) * pw.at(r - l) % m + m) % m; } else{ return (sh.at(r) - (__int128_t)sh.at(l) * pw.at(l - r) % m + m) % m; } } ll range_hash(vector> p){ ll now = 0, ret = 0; reverse(all(p)); for(int i = 0;i < p.size();i++){ ret += (__int128_t)hash(p.at(i).fi, p.at(i).se) * pw.at(now) % m; ret %= m; now += abs(p.at(i).fi - p.at(i).se); } return ret; } ll range_hash(vector> p){ ll now = 0, ret = 0; reverse(all(p)); for(int i = 0;i < p.size();i++){ ret += (__int128_t)hash(p.at(i).fi, p.at(i).se) * pw.at(now) % m; ret %= m; now += abs(p.at(i).fi - p.at(i).se); } return ret; } }; class CentroidDecomposition { public: int V; vector > G; vector used; //sz:重心分解後の最大部分木に含まれる頂点の数(自分を含める) //par:重心分解後の親の頂点 vector sz, par; //部分木のサイズを計算 void calcSize(ll u,ll p){ sz[u] = 1; for(ll v : G[u]){ if(!used[v] && v != p){ calcSize(v,u); sz[u] += sz[v]; } } } void cdBuild(ll u,ll p){ calcSize(u,-1); ll tot = sz[u]; bool ok = false; ll pp = -1; //いま見ている部分木での重心を見つける while(!ok){ ok = true; for(ll v : G[u]){ if(!used[v] && v != pp && 2*sz[v] > tot){ pp = u, u = v, ok = false; break; } } } par[u] = p; //何らかの操作 used[u] = true; //深さ優先でたどる for(ll v : G[u]){ if(!used[v]){ cdBuild(v,u); } } } CentroidDecomposition(ll node_size) : V(node_size), G(V), used(V, false) , sz(V, 0), par(V, -1){} void add_edge(ll u,ll v){ G[u].push_back(v), G[v].push_back(u); } void build(){ cdBuild(0,-1); } }; template istream &operator>>(istream &a, pair &b){ a >> b.fi >> b.se; return a; } template ostream& operator<<(ostream& os, const pair& p) { return os << p.first << ' ' << p.second; } template istream &operator>>(istream &a, vector &b){ rep(i, b.size()){ a >> b.at(i); } return a; } template ostream& operator<<(ostream& os, const vector& v) { for (int i = 0; i < (int)v.size(); i++) { if (i) os << ' '; os << v[i]; } return os; } template ostream &operator<<(ostream &a, vector> &b){ rep(i, b.size()){ a << b.at(i); if(i != b.size() - 1){ a< istream &operator>>(istream &a, set &b){ T tmp; a >> tmp; b.insert(tmp); return a; } void scan(){} template void scan(T &head, S &... plus){ cin >> head; scan(plus...); } void out(){} template void out(const T &head, const S &... plus){ cout << head << endl; out(plus...); } //行列の乗算 vvl matrix_mult(vvl X, vvl Y) { vvl Z(X.size(), vl(Y[0].size())); rep(i, X.size()) { rep(k, Y.size()) { rep(j, Y[0].size()) { Z[i][j] = (Z[i][j] + X[i][k] * Y[k][j]) % mod; } } } return Z; } //A^nの計算 vvl matrix_pow(vvl A, ll n) { vvl B(A.size(), vl(A[0].size())); //単位行列でBを初期化 rep(i, B.size()) { B[i][i] = 1; } while (n>0) { if (n & 1) { B = matrix_mult(B, A); } A = matrix_mult(A, A); n = n >> 1; } return B; } template< typename T > struct FormalPowerSeries : vector< T > { using vector< T >::vector; using P = FormalPowerSeries; P pre(int deg) const { return P(begin(*this), begin(*this) + min((int) this->size(), deg)); } P rev(int deg = -1) const { P ret(*this); if(deg != -1) ret.resize(deg, T(0)); reverse(begin(ret), end(ret)); return ret; } void shrink() { while(this->size() && this->back() == T(0)) this->pop_back(); } P operator+(const P &r) const { return P(*this) += r; } P operator+(const T &v) const { return P(*this) += v; } P operator-(const P &r) const { return P(*this) -= r; } P operator-(const T &v) const { return P(*this) -= v; } P operator*(const P &r) const { return P(*this) *= r; } P operator*(const T &v) const { return P(*this) *= v; } P operator/(const P &r) const { return P(*this) /= r; } P operator%(const P &r) const { return P(*this) %= r; } P &operator+=(const P &r) { if(r.size() > this->size()) this->resize(r.size()); for(int i = 0; i < r.size(); i++) (*this)[i] += r[i]; return *this; } P &operator-=(const P &r) { if(r.size() > this->size()) this->resize(r.size()); for(int i = 0; i < r.size(); i++) (*this)[i] -= r[i]; return *this; } P &operator*=(const P &r) { if(this->empty() || r.empty()) { this->clear(); return *this; } auto ret = convolution(*this, r); return *this = {begin(ret), end(ret)}; } P &operator/=(const P &r) { if(this->size() < r.size()) { this->clear(); return *this; } int n = this->size() - r.size() + 1; return *this = (rev().pre(n) * r.rev().inv(n)).pre(n).rev(n); } P &operator%=(const P &r) { return *this -= *this / r * r; } // https://judge.yosupo.jp/problem/division_of_polynomials pair< P, P > div_mod(const P &r) { P q = *this / r; return make_pair(q, *this - q * r); } P operator-() const { P ret(this->size()); for(int i = 0; i < this->size(); i++) ret[i] = -(*this)[i]; return ret; } P &operator+=(const T &r) { if(this->empty()) this->resize(1); (*this)[0] += r; return *this; } P &operator-=(const T &r) { if(this->empty()) this->resize(1); (*this)[0] -= r; return *this; } P &operator*=(const T &v) { for(int i = 0; i < this->size(); i++) (*this)[i] *= v; return *this; } P dot(P r) const { P ret(min(this->size(), r.size())); for(int i = 0; i < ret.size(); i++) ret[i] = (*this)[i] * r[i]; return ret; } P operator>>(int sz) const { if(this->size() <= sz) return {}; P ret(*this); ret.erase(ret.begin(), ret.begin() + sz); return ret; } P operator<<(int sz) const { P ret(*this); ret.insert(ret.begin(), sz, T(0)); return ret; } T operator()(T x) const { T r = 0, w = 1; for(auto &v : *this) { r += w * v; w *= x; } return r; } P diff() const { const int n = (int) this->size(); P ret(max(0, n - 1)); for(int i = 1; i < n; i++) ret[i - 1] = (*this)[i] * T(i); return ret; } P integral() const { const int n = (int) this->size(); P ret(n + 1); ret[0] = T(0); for(int i = 0; i < n; i++) ret[i + 1] = (*this)[i] / T(i + 1); return ret; } // https://judge.yosupo.jp/problem/inv_of_formal_power_series // F(0) must not be 0 P inv(int deg = -1) const { assert(((*this)[0]) != T(0)); const int n = (int) this->size(); if(deg == -1) deg = n; P ret({T(1) / (*this)[0]}); for(int i = 1; i < deg; i <<= 1) { ret = (ret + ret - ret * ret * pre(i << 1)).pre(i << 1); } return ret.pre(deg); } // https://judge.yosupo.jp/problem/log_of_formal_power_series // F(0) must be 1 P log(int deg = -1) const { assert((*this)[0] == T(1)); const int n = (int) this->size(); if(deg == -1) deg = n; return (this->diff() * this->inv(deg)).pre(deg - 1).integral(); } // https://judge.yosupo.jp/problem/sqrt_of_formal_power_series P sqrt(int deg = -1, const function< T(T) > &get_sqrt = [](T) { return T(1); }) const { const int n = (int) this->size(); if(deg == -1) deg = n; if((*this)[0] == T(0)) { for(int i = 1; i < n; i++) { if((*this)[i] != T(0)) { if(i & 1) return {}; if(deg - i / 2 <= 0) break; auto ret = (*this >> i).sqrt(deg - i / 2, get_sqrt); if(ret.empty()) return {}; ret = ret << (i / 2); if(ret.size() < deg) ret.resize(deg, T(0)); return ret; } } return P(deg, 0); } auto sqr = T(get_sqrt((*this)[0])); if(sqr * sqr != (*this)[0]) return {}; P ret{sqr}; T inv2 = T(1) / T(2); for(int i = 1; i < deg; i <<= 1) { ret = (ret + pre(i << 1) * ret.inv(i << 1)) * inv2; } return ret.pre(deg); } P sqrt(const function< T(T) > &get_sqrt, int deg = -1) const { return sqrt(deg, get_sqrt); } // https://judge.yosupo.jp/problem/exp_of_formal_power_series // F(0) must be 0 P exp(int deg = -1) const { if(deg == -1) deg = this->size(); assert((*this)[0] == T(0)); const int n = (int) this->size(); if(deg == -1) deg = n; P ret({T(1)}); for(int i = 1; i < deg; i <<= 1) { ret = (ret * (pre(i << 1) + T(1) - ret.log(i << 1))).pre(i << 1); } return ret.pre(deg); } // https://judge.yosupo.jp/problem/pow_of_formal_power_series P pow(int64_t k, int deg = -1) const { const int n = (int) this->size(); if(deg == -1) deg = n; if(k == 0) { P ret(deg, T(0)); ret[0] = T(1); return ret; } for(int i = 0; i < n; i++) { if(i * k > deg) return P(deg, T(0)); if((*this)[i] != T(0)) { T rev = T(1) / (*this)[i]; P ret = (((*this * rev) >> i).log() * k).exp() * ((*this)[i].pow(k)); ret = (ret << (i * k)).pre(deg); if(ret.size() < deg) ret.resize(deg, T(0)); return ret; } } return *this; } // https://yukicoder.me/problems/no/215 P mod_pow(int64_t k, P g) const { P modinv = g.rev().inv(); auto get_div = [&](P base) { if(base.size() < g.size()) { base.clear(); return base; } int n = base.size() - g.size() + 1; return (base.rev().pre(n) * modinv.pre(n)).pre(n).rev(n); }; P x(*this), ret{1}; while(k > 0) { if(k & 1) { ret *= x; ret -= get_div(ret) * g; ret.shrink(); } x *= x; x -= get_div(x) * g; x.shrink(); k >>= 1; } return ret; } // https://judge.yosupo.jp/problem/polynomial_taylor_shift P taylor_shift(T c) const { int n = (int) this->size(); vector< T > fact(n), rfact(n); fact[0] = rfact[0] = T(1); for(int i = 1; i < n; i++) fact[i] = fact[i - 1] * T(i); rfact[n - 1] = T(1) / fact[n - 1]; for(int i = n - 1; i > 1; i--) rfact[i - 1] = rfact[i] * T(i); P p(*this); for(int i = 0; i < n; i++) p[i] *= fact[i]; p = p.rev(); P bs(n, T(1)); for(int i = 1; i < n; i++) bs[i] = bs[i - 1] * c * rfact[i] * fact[i - 1]; p = (p * bs).pre(n); p = p.rev(); for(int i = 0; i < n; i++) p[i] *= rfact[i]; return p; } }; template< typename Mint > using FPS = FormalPowerSeries< Mint >; // Bostan-Mori // find [x^N] P(x)/Q(x), O(K log K log N) // deg(Q(x)) = K, deg(P(x)) < K, Q[0] = 1 template mint BostanMori(const FPS &P, const FPS &Q, long long N) { assert(!P.empty() && !Q.empty()); if (N == 0) return P[0] / Q[0]; int qdeg = (int)Q.size(); FPS P2{P}, minusQ{Q}; P2.resize(qdeg - 1); for (int i = 1; i < (int)Q.size(); i += 2) minusQ[i] = -minusQ[i]; P2 *= minusQ; FPS Q2 = Q * minusQ; FPS S(qdeg - 1), T(qdeg); for (int i = 0; i < (int)S.size(); ++i) { S[i] = (N % 2 == 0 ? P2[i * 2] : P2[i * 2 + 1]); } for (int i = 0; i < (int)T.size(); ++i) { T[i] = Q2[i * 2]; } return BostanMori(S, T, N >> 1); } // Union-Find, we can undo struct UnionFind { // core member vector par; stack> history; // constructor UnionFind() {} UnionFind(int n) : par(n, -1) { } void init(int n) { par.assign(n, -1); } // core methods int root(int x) { if (par[x] < 0) return x; else return root(par[x]); } bool same(int x, int y) { return root(x) == root(y); } bool merge(int x, int y) { x = root(x), y = root(y); history.emplace(x, par[x]); history.emplace(y, par[y]); if (x == y) return false; if (par[x] > par[y]) swap(x, y); // merge technique par[x] += par[y]; par[y] = x; return true; } int size(int x) { return -par[root(x)]; } // 1-step undo void undo() { for (int iter = 0; iter < 2; ++iter) { par[history.top().first] = history.top().second; history.pop(); } } // erase history void snapshot() { while (!history.empty()) history.pop(); } // all rollback void rollback() { while (!history.empty()) undo(); } // debug friend ostream& operator << (ostream &s, UnionFind uf) { map> groups; for (int i = 0; i < uf.par.size(); ++i) { int r = uf.root(i); groups[r].push_back(i); } for (const auto &it : groups) { s << "group: "; for (auto v : it.second) s << v << " "; s << endl; } return s; } }; istream &operator>>(istream &a, mint &b){ ll tmp; a >> tmp; b = tmp; return a; } ostream& operator<<(ostream& os, const mint& x) { return os << x.val(); } istream &operator>>(istream &a, mmint &b){ ll tmp; a >> tmp; b = tmp; return a; } ostream& operator<<(ostream& os, const mmint& x) { return os << x.val(); } #line 1 "math/number-theory/generalized-floor-sum-pq-le-2.hpp" // 一般化 floor sum (p+q<=2) のうち (0,1),(1,1),(0,2) をまとめて求める。 // ans_01 = Σ floor((a i + b)/m), ans_11 = Σ i*floor((a i + b)/m), ans_02 = Σ floor((a i + b)/m)^2。 // n>=0, m>0 を仮定する。a,b は負でもよい。 // 計算量 O(log m)。 #include #include template struct GeneralizedFloorSumPQLe2Result { T ans_01; T ans_11; T ans_02; }; namespace generalized_floor_sum_pq_le_2_internal { template struct is_integral : std::is_integral {}; #ifdef __SIZEOF_INT128__ template <> struct is_integral<__int128_t> : std::true_type {}; template <> struct is_integral<__uint128_t> : std::true_type {}; #endif template struct is_signed : std::is_signed {}; #ifdef __SIZEOF_INT128__ template <> struct is_signed<__int128_t> : std::true_type {}; template <> struct is_signed<__uint128_t> : std::false_type {}; #endif template T floor_div(T x, T y) { assert(y > 0); if constexpr (is_signed::value) { T q = x / y; T r = x % y; if (r < 0) --q; return q; } else { return x / y; } } template T floor_mod(T x, T y) { assert(y > 0); if constexpr (is_signed::value) { T r = x % y; if (r < 0) r += y; return r; } else { return x % y; } } // Σ_{i=0}^{n-1} i = n(n-1)/2 template T sum_0_to_n_minus_1(T n) { if (n == 0) return 0; if ((n & 1) == 0) return (n / 2) * (n - 1); return n * ((n - 1) / 2); } // Σ_{i=0}^{n-1} i^2 = (n-1)n(2n-1)/6 template T sum_0_to_n_minus_1_sq(T n) { if (n == 0) return 0; T a = n - 1, b = n, c = 2 * n - 1; if ((a % 2) == 0) a /= 2; else if ((b % 2) == 0) b /= 2; else c /= 2; if ((a % 3) == 0) a /= 3; else if ((b % 3) == 0) b /= 3; else c /= 3; return a * b * c; } template T sum_range(T l, T r) { // Σ_{i=l}^{r} i if (l > r) return 0; T cnt = r - l + 1; T s = l + r; if ((s & 1) == 0) s /= 2; else cnt /= 2; return s * cnt; } template struct Result { Int ans_01; Int ans_11; Int ans_02; }; template Result solve(Int n, Int m, Int a, Int b) { if constexpr (is_signed::value) assert(n >= 0); assert(m > 0); if (n == 0) return {0, 0, 0}; const Int qa = floor_div(a, m); a = floor_mod(a, m); const Int qb = floor_div(b, m); b = floor_mod(b, m); if constexpr (is_signed::value) { assert(a >= 0); assert(b >= 0); } assert(a < m); assert(b < m); Result base = {0, 0, 0}; if (a != 0) { const Int y_max = (a * n + b) / m; if (y_max != 0) { const Int x_max = y_max * m - b; const Int t = (x_max + a - 1) / a; // ceil(x_max / a) const Int b2 = (a - (x_max % a)) % a; const auto rec = solve(y_max, a, m, b2); const Int head_01 = rec.ans_01; const Int head_11 = ((2 * t - 1) * rec.ans_01 - rec.ans_02) / 2; const Int head_02 = (2 * y_max - 1) * rec.ans_01 - 2 * rec.ans_11; const Int tail_len = n - t; const Int tail_01 = tail_len * y_max; const Int tail_11 = y_max * sum_range(t, n - 1); const Int tail_02 = tail_len * y_max * y_max; base = {head_01 + tail_01, head_11 + tail_11, head_02 + tail_02}; } } const Int si = sum_0_to_n_minus_1(n); const Int si2 = sum_0_to_n_minus_1_sq(n); Result ans; ans.ans_01 = qa * si + qb * n + base.ans_01; ans.ans_11 = qa * si2 + qb * si + base.ans_11; ans.ans_02 = qa * qa * si2 + 2 * qa * qb * si + 2 * qa * base.ans_11 + qb * qb * n + 2 * qb * base.ans_01 + base.ans_02; return ans; } } // namespace generalized_floor_sum_pq_le_2_internal template GeneralizedFloorSumPQLe2Result generalized_floor_sum_pq_le_2(T n, T m, T a, T b) { static_assert(generalized_floor_sum_pq_le_2_internal::is_integral::value, "T must be integer."); if constexpr (generalized_floor_sum_pq_le_2_internal::is_signed::value) assert(n >= 0); assert(m > 0); #ifdef __SIZEOF_INT128__ using I = __int128_t; const auto res = generalized_floor_sum_pq_le_2_internal::solve( static_cast(n), static_cast(m), static_cast(a), static_cast(b)); return {static_cast(res.ans_01), static_cast(res.ans_11), static_cast(res.ans_02)}; #else const auto res = generalized_floor_sum_pq_le_2_internal::solve(n, m, a, b); return {res.ans_01, res.ans_11, res.ans_02}; #endif } using i128 = __int128_t; vector> matrixmul(vector> X, vector> Y) { vector> Z(4); Z = { X[0] * Y[0] + X[1] * Y[2], X[0] * Y[1] + X[1] * Y[3], X[2] * Y[0] + X[3] * Y[2], X[2] * Y[1] + X[3] * Y[3] }; return Z; } ll op(ll a, ll b){ return a + b; } ll e(){ return 0; } ll mapping(ll a, ll b){ return (a == inf ? b : a); } ll composition(ll a, ll b){ return (a == inf ? b : a); } ll id(){ return inf; } struct ss{ long long value; int size; }; ss op2(ss a, ss b){ return {a.value+b.value, a.size+b.size}; } ss e2(){ return {0, 0}; } ss mapping2(ll a, ss b){ if(a != inf) b.value = a*b.size; return b; } ll composition2(ll a, ll b){ return (a == inf ? b : a); } ll id2(){ return inf; } clock_t st; void beg(ll ep, ll mo, bool nc, bool fac){ st = clock(); mod = mo; ios::sync_with_stdio(false); std::cin.tie(nullptr); cout << fixed << setprecision(20); if(nc)init();//二項係数 階乗(fact[]) srand((unsigned int)time(NULL)); if(fac)prepare_factorize();//素因数分解 rep(i, ep){ eps /= 10; } } ll kp = 0; bool ff(ll x){ if(x >= kp){ return true; } else{ return false; } } std::random_device rd; std::mt19937_64 mt(rd()); void nowtime(){ cout << ((double)(clock()) - (double)(st)) / CLOCKS_PER_SEC << endl; } double nowtimed(){ return ((double)(clock()) - (double)(st)) / CLOCKS_PER_SEC; } ll bits_msb(ll v ) { v = v | (v >> 1); v = v | (v >> 2); v = v | (v >> 4); v = v | (v >> 8); v = v | (v >> 16); v = v | (v >> 32); return v ^ (v >> 1); } ll ssqr(ll a){ ll ok = 0, ng = 2000000000; while(ng - ok > 1){ ll mi = (ok + ng) / 2; if(mi * mi <= a){ ok = mi; } else{ ng = mi; } } return ok; } vs ans(663, string(663, '.')); void pt(int x, int y){ x = 2 * x + 1; y = 2 * y + 1; ans[x][y] = '#'; ans[x + 1][y] = '#'; ans[x - 1][y] = '#'; ans[x][y - 1] = '#'; ans[x][y + 1] = '#'; } void solve(){ rep(i, 165){ pt(164, i); } rep(i, 165){ rep(j, 55){ pt(i, j * 3); } } rep(i, 165){ pt(166, i + 165); } rep(i, 165){ rep(j, 55){ pt(i + 166, j * 3 + 165); } } rep(i, 3){ rep(j, 2){ pt(164 + i, 164 + j); } } cout << 663 << endl; for(auto x : ans)cout << x << endl; ans[330][1] = '.'; for(auto x : ans)cout << x << endl; ans[330][1] = '#'; for(auto x : ans)cout << x << endl; return; } int main() { beg(12, 998244353, false, false);// eps mod nck factorize max変更忘れない! define外せ! ll T; T = 1; //cin >> T; while(T--){ solve(); } }