**ICPC NOTEBOOK**

**School : Le Hong Phong High School for The Gifted**
**Team : No_Zone_for_Loser**
**Member : Nguyễn Tấn Minh (11CTIN-A), Đặng Thành Nghĩa (11CTIN-A), Lê Lâm (10CTIN)**

[TOC]
# Main template

#include<bits/stdc++.h>

using namespace std;

#define range(v) begin(v), end(v)
#define compact(v) v.erase(unique(range(v)), end(v))
#define len(x) (int)(x).size()
#define BIT(x, i) ((x)>>(i)&1)
#define MASK(x) (1LL<<(x))
#define forw(i,j,z) for(int i=(int)j;i<=(int)z;i++)
#define ford(i,j,z) for(int i=(int)j;i>=(int)z;i--)
#define mp make_pair
//define lucian nghiaaa's cuties
//define zero_op nghiaaa's cuties
//define No_Zone_for_Loser winner team 

typedef long long ll;
typedef long double ld;
typedef unsigned long long ull;

template<class T> bool maximize(T& a, const T& b){ if(a<b) return a=b, true; return false; }
template<class T> bool minimize(T& a, const T& b){ if(a>b) return a=b, true; return false; }

void No_Zone_for_Loser(){

}

signed main(){
  ios_base::sync_with_stdio(0);
  cin.tie(NULL); cout.tie(NULL);

  #define task "a"
  if(fopen(task".inp", "r")){
    freopen(task".inp", "r", stdin);
    freopen(task".out", "w", stdout);
  }

  int TC=1; //cin>>TC;
  while(TC--) No_Zone_for_Loser();

  return 0;
}

# MATH

## Miller Rabin

// From https://github.com/SnapDragon64/ContestLibrary/blob/master/math.h
// which also has specialized versions for 32-bit and 42-bit
//
// Tested:
// - https://oj.vnoi.info/problem/icpc22_national_c (fastest solution)
// - https://www.spoj.com/problems/PON/

// Rabin miller {{{
inline uint64_t mod_mult64(uint64_t a, uint64_t b, uint64_t m) {
    return __int128_t(a) * b % m;
}
uint64_t mod_pow64(uint64_t a, uint64_t b, uint64_t m) {
    uint64_t ret = (m > 1);
    for (;;) {
        if (b & 1) ret = mod_mult64(ret, a, m);
        if (!(b >>= 1)) return ret;
        a = mod_mult64(a, a, m);
    }
}

// Works for all primes p < 2^64
bool is_prime(uint64_t n) {
    if (n <= 3) return (n >= 2);
    static const uint64_t small[] = {
        2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67,
        71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139,
        149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199,
    };
    for (size_t i = 0; i < sizeof(small) / sizeof(uint64_t); ++i) {
        if (n % small[i] == 0) return n == small[i];
    }

    // Makes use of the known bounds for Miller-Rabin pseudoprimes.
    static const uint64_t millerrabin[] = {
        2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37,
    };
    static const uint64_t A014233[] = {  // From OEIS.
        2047LL, 1373653LL, 25326001LL, 3215031751LL, 2152302898747LL,
        3474749660383LL, 341550071728321LL, 341550071728321LL,
        3825123056546413051LL, 3825123056546413051LL, 3825123056546413051LL, 0,
    };
    uint64_t s = n-1, r = 0;
    while (s % 2 == 0) {
        s /= 2;
        r++;
    }
    for (size_t i = 0, j; i < sizeof(millerrabin) / sizeof(uint64_t); i++) {
        uint64_t md = mod_pow64(millerrabin[i], s, n);
        if (md != 1) {
            for (j = 1; j < r; j++) {
                if (md == n-1) break;
                md = mod_mult64(md, md, n);
            }
            if (md != n-1) return false;
        }
        if (n < A014233[i]) return true;
    }
    return true;
}
// }}}

## FFT

using cd=complex<double>;
const double PI=acos(-1);

void fft(vector<cd>& a, bool is_inverse){
    int n=a.size();
    int lg=0; while((1<<lg)<n) ++lg;

    for(int i=1, j=0; i<n; ++i){
        int bit=n>>1;
        for(; j&bit; bit>>=1) j^=bit;
        j^=bit;
        if(i<j) swap(a[i], a[j]);
    }

    for(int len=2; len<=n; len<<=1){
        double ang=2.0*PI/len*(is_inverse?-1:1);
        cd wlen(cos(ang), sin(ang));
        for(int i=0; i<n; i+=len){
            cd w(1);
            for(int j=0; j*2<len; ++j){
                cd u=a[i+j], v=a[i+j+len/2]*w;
                a[i+j]=u+v;
                a[i+j+len/2]=u-v;
                w*=wlen;
            }
        }
    }
    if(is_inverse){
        for(cd &x:a) x/=n;
    }
}

vector<int> multiply(const vector<int>& a, const vector<int>& b){
    vector<cd> fa(all(a)), fb(all(b));
    int n=1; while(a.size()+b.size()>n) n<<=1;

    fa.resize(n), fb.resize(n);
    fft(fa, 0), fft(fb, 0);
    for(int i=0; i<n; ++i) fa[i]*=fb[i];
    fft(fa, 1);

    vector<int> result(n);
    for(int i=0; i<n; ++i) result[i]=round(fa[i].real());
    return result;
}

## NTT

// NTT {{{
//
// Faster than NTT_chemthan.h
//
// Usage:
// auto c = multiply(a, b);
// where a and b are vector<ModInt<ANY_MOD>>
// (If mod is NOT NTT_PRIMES, it does 3 NTT and combine result)

constexpr int NTT_PRIMES[] = {998244353, 167772161, 469762049};

// assumptions:
// - |a| is power of 2
// - mint::mod() is a valid NTT primes (2^k * m + 1)
template<typename mint> void ntt(std::vector<mint>& a, bool is_inverse) {
    int n = a.size();
    if (n == 1) return;

    static const int mod = mint::mod();
    static const mint root = mint::get_primitive_root();
    assert(__builtin_popcount(n) == 1 && (mod - 1) % n == 0);

    static std::vector<mint> w{1}, iw{1};
    for (int m = w.size(); m < n / 2; m *= 2) {
        mint dw = root.pow((mod - 1) / (4 * m));
        mint dwinv = dw.inv();
        w.resize(m * 2);
        iw.resize(m * 2);
        for (int i = 0; i < m; ++i) {
            w[m + i] = w[i] * dw;
            iw[m + i] = iw[i] * dwinv;
        }
    }

    if (!is_inverse) {
        for (int m = n; m >>= 1; ) {
            for (int s = 0, k = 0; s < n; s += 2 * m, ++k) {
                for (int i = s; i < s + m; ++i) {
                    mint x = a[i], y = a[i + m] * w[k];
                    a[i] = x + y;
                    a[i + m] = x - y;
                }
            }
        }
    } else {
        for (int m = 1; m < n; m *= 2) {
            for (int s = 0, k = 0; s < n; s += 2 * m, ++k) {
                for (int i = s; i < s + m; ++i) {
                    mint x = a[i], y = a[i + m];
                    a[i] = x + y;
                    a[i + m] = (x - y) * iw[k];
                }
            }
        }
        int n_inv = mint(n).inv().x;
        for (auto& v : a) v *= n_inv;
    }
}

template<typename mint>
std::vector<mint> ntt_multiply(int sz, std::vector<mint> a, std::vector<mint> b) {
    a.resize(sz);
    b.resize(sz);
    if (a == b) { // optimization for squaring polynomial
        ntt(a, false);
        b = a;
    } else {
        ntt(a, false);
        ntt(b, false);
    }
    for (int i = 0; i < sz; ++i) a[i] *= b[i];
    ntt(a, true);
    return a;
}

template<int MOD, typename mint>
std::vector<ModInt<MOD>> convert_mint_and_multiply(
        int sz,
        const std::vector<mint>& a,
        const std::vector<mint>& b) {
    using mint2 = ModInt<MOD>;

    std::vector<mint2> a2(a.size()), b2(b.size());
    for (size_t i = 0; i < a.size(); ++i) {
        a2[i] = mint2(a[i].x);
    }
    for (size_t i = 0; i < b.size(); ++i) {
        b2[i] = mint2(b[i].x);
    }
    return ntt_multiply(sz, a2, b2);
}

long long combine(int r0, int r1, int r2, int mod) {
    using mint2 = ModInt<NTT_PRIMES[2]>;
    static const long long m01 = 1LL * NTT_PRIMES[0] * NTT_PRIMES[1];
    static const long long m0_inv_m1 = ModInt<NTT_PRIMES[1]>(NTT_PRIMES[0]).inv().x;
    static const long long m01_inv_m2 = mint2(m01).inv().x;

    int v1 = (m0_inv_m1 * (r1 + NTT_PRIMES[1] - r0)) % NTT_PRIMES[1];
    auto v2 = (mint2(r2) - r0 - mint2(NTT_PRIMES[0]) * v1) * m01_inv_m2;
    return (r0 + 1LL * NTT_PRIMES[0] * v1 + m01 % mod * v2.x) % mod;
}

template<typename mint>
std::vector<mint> multiply(const std::vector<mint>& a, const std::vector<mint>& b) {
    if (a.empty() || b.empty()) return {};
    int sz = 1, sz_a = a.size(), sz_b = b.size();
    while (sz < sz_a + sz_b) sz <<= 1;
    if (sz <= 16) {
        std::vector<mint> res(sz_a + sz_b - 1);
        for (int i = 0; i < sz_a; ++i) {
            for (int j = 0; j < sz_b; ++j) {
                res[i + j] += a[i] * b[j];
            }
        }
        return res;
    }

    int mod = mint::mod();
    std::vector<mint> res;
    if (std::find(std::begin(NTT_PRIMES), std::end(NTT_PRIMES), mod) != std::end(NTT_PRIMES)) {
        res = ntt_multiply(sz, a, b);
    } else {
        auto c0 = convert_mint_and_multiply<NTT_PRIMES[0], mint> (sz, a, b);
        auto c1 = convert_mint_and_multiply<NTT_PRIMES[1], mint> (sz, a, b);
        auto c2 = convert_mint_and_multiply<NTT_PRIMES[2], mint> (sz, a, b);

        res.resize(sz);
        for (int i = 0; i < sz; ++i) {
            res[i] = combine(c0[i].x, c1[i].x, c2[i].x, mod);
        }
    }

    res.resize(sz_a + sz_b - 1);
    return res;
}
// }}}

## Lucas Theorem

//Calculate C(n, k) modulo M with smal M and big n, k
const int M=1e5+3;

int sub(int a, int b){ return (a<b?a-b+M:a-b); }
int add(int a, int b){ return (a+b<M?a+b:a+b-M); }
int mul(int a, int b){ return 1LL*a*b%M; }
int binpow(int a, int b){ int ans=1; while(b){ if(b&1) ans=mul(ans, a); a=mul(a, a); b>>=1; } return ans; }
int inverse(int a){ return binpow(a, M-2); }

int fact[M], inv[M];

int C(int NN, int KK){
  if(NN<KK) return 0;
  return mul(fact[NN], mul(inv[KK], inv[NN-KK]));
}

int Lucas(ll n, ll m){
  if(n<m) return 0;

  int ans=1;
  for(int i=0; i<4; ++i){
    int mod_n=n%M, mod_m=m%M;
    n/=M; m/=M;
    ans=mul(ans, C(mod_n, mod_m));
  }
  return ans;
}

# DP Trick

## Line Container

struct Line {
  mutable ll m, b, p;
  bool operator<(const Line& o) const { return m < o.m; }
  bool operator<(ll x) const { return p < x; }
};
 
struct LineContainer : multiset<Line, less<>> {
  // (for doubles, use inf = 1/.0, div(a,b) = a/b)
  // Caution : This template solve for finding MAX, so if want to find MIN, NEGATE ALL THE LINE FUNCTION
    const ll inf = LLONG_MAX;
  ll div(ll a, ll b) { // floored division
    return a/b-((a^b)<0&&a%b); }
  bool isect(iterator x, iterator y) {
    if (y==end()) { x->p=inf; return false; }
    if (x->m==y->m) x->p=x->b>y->b?inf:-inf;
    else x->p=div(y->b-x->b, x->m-y->m);
    return x->p>=y->p;
  }
  void add(ll m, ll b) {
    auto z=insert({m, b, 0}), y=z++, x=y;
    while(isect(y, z)) z=erase(z);
    if (x!=begin() && isect(--x, y)) isect(x, y = erase(y));
    while((y = x)!=begin() && (--x)->p>=y->p)
    isect(x, erase(y));
  }
  ll query(ll x) { // max value at x
    assert(!empty());
    auto l=*lower_bound(x);
    return l.m*x+l.b;
  }
};

# GRAPH

## Gomory Hu

// Usage:
//   GomoryHu g(n);
//   g.addEdge(u, v, c)  // 1-direction
//   g.calc()
//   --> g.answer[i][j] = min cut i-j
//
// Note:
// - When i -> j is not connected, answer[i][j] = INF
// - Used together with MaxFlowDinic.h
//
// Tested:
// - https://oj.vnoi.info/problem/mcquery
// - https://oj.vnoi.info/problem/icpc21_mn_a

/*
 * Find min cut between every pair of vertices using N max_flow call (instead of N^2)
 * Not tested with directed graph
 * Index start from 0
 */
struct GomoryHu {
    int ok[MN], cap[MN][MN];
    int answer[MN][MN], parent[MN];
    int n;
    MaxFlow flow;

    GomoryHu(int _n) : n(_n), flow(_n) {
        for(int i = 0; i < n; ++i) ok[i] = parent[i] = 0;
        for(int i = 0; i < n; ++i)
            for(int j = 0; j < n; ++j)
                cap[i][j] = 0, answer[i][j] = INF;
    }

    void addEdge(int u, int v, int c) {
        cap[u][v] += c;
    }

    void calc() {
        for(int i = 1; i <= n-1; ++i) {
            flow = MaxFlow(n);
            REP(u,n) REP(v,n)
                if (cap[u][v])
                    flow.addEdge(u, v, cap[u][v]);

            int f = flow.getMaxFlow(i, parent[i]);

            bfs(i);
            for(int j = i+1; j < n; ++j)
                if (ok[j] && parent[j]==parent[i])
                    parent[j]=i;
            
            answer[i][parent[i]] = answer[parent[i]][i] = f;
            for(int j = 0; j < i; ++j)
                answer[i][j]=answer[j][i]=min(f,answer[parent[i]][j]);
        }
    }
    void bfs(int start) {
        memset(ok,0,sizeof ok);
        queue<int> qu;
        qu.push(start);
        while (!qu.empty()) {
            int u=qu.front(); qu.pop();
            for(int xid = 0; xid < (int) flow.g[u].size(); ++xid) {
                int id = flow.g[u][xid];
                int v = flow.e[id].b, fl = flow.e[id].flow, c = flow.e[id].cap;
                if (!ok[v] && fl < cap) {
                    ok[v]=1;
                    qu.push(v);
                }
            }
        }
    }
};

# DATA STRUCTURE

## Manhattan MST

struct Point{
  int x, y, id;
  Point(int x=0, int y=0, int id=0):x(x), y(y), id(id){}
};

struct Edge{
  int u, v, w;
  Edge(int u=0, int v=0, int w=0):u(u), v(v), w(w){}
  bool operator<(const Edge& other) const{
    return w<other.w;
  }
};

int n, lab[N], ID[N], EDGE_SZ;
Point P[N];
Edge E[N*4];

int root(int u){ return (lab[u]<0?u:(lab[u]=root(lab[u]))); }
bool unite(int u, int v){
   u=root(u), v=root(v); if(u==v) return false;
   if(lab[u]>lab[v]) swap(u, v); lab[u]+=lab[v]; lab[v]=u; return true;
}

void get_edge(){
  map<int, int> M;
  iota(ID+1, ID+1+n, 1);
  sort(ID+1, ID+1+n, [&](int u, int v){
    return (P[u].x+P[u].y)<(P[v].x+P[v].y);
  });
  for(int T=1, i; T<=n; ++T){
    i=ID[T];
    for(auto iter=M.lower_bound(-P[i].x); iter!=M.end(); M.erase(iter++)){
      if((P[iter->second].y-P[i].y)>(P[iter->second].x-P[i].x)) break;
      E[EDGE_SZ++]=Edge(iter->second, i, P[i].y-P[iter->second].y+P[i].x-P[iter->second].x);
    }
    M[-P[i].x]=i;
  }
}

void solve(){
  cin>>n;
  for(int i=1; i<=n; ++i){
    cin>>P[i].x>>P[i].y;
    P[i].id=i;
  }
  fill(lab+1, lab+1+n, -1);

  for(int d=0; d<4; ++d){
    get_edge();
    for(int i=1; i<=n; ++i){
      if(d&1) P[i].x*=-1;
      else swap(P[i].x, P[i].y);
    }
  }

  sort(E, E+EDGE_SZ);
  ll ans=0; vector<pair<int, int>> pairs;
  for(int i=0; i<EDGE_SZ; ++i){
    if(unite(E[i].u, E[i].v)){
      ans+=E[i].w;
      pairs.push_back({E[i].u-1, E[i].v-1});
    }
  }
  cout<<ans<<ln;
  for(auto [u, v]:pairs) cout<<u<<" "<<v<<ln;
}

## Wavelet Tree

// WaveletMatrix {{{
// Copied from https://github.com/dacin21/dacin21_codebook/blob/master/trees/wavelet_matrix.cpp
//
// Notes:
// - Index from 0
// - k (for k-th query) from 0
// - Need to remove #define int long long
//
// Tested:
// - (kth query) https://judge.yosupo.jp/problem/range_kth_smallest
// - (range_count) https://judge.yosupo.jp/problem/static_range_frequency
class Bit_Presum{
public:
    static constexpr uint32_t omega = CHAR_BIT * sizeof(uint64_t);
    static constexpr uint32_t lg_omega = __lg(omega);
    static_assert(omega == 64u);

    Bit_Presum(vector<uint64_t> mask_) : n(mask_.size()), mask(move(mask_)), presum(n+1){
        build();
    }
    Bit_Presum(uint32_t bits, bool init_val = 0) : n((bits>>lg_omega) + 1), mask(n, init_val ? ~uint64_t{0} : uint64_t{0}), presum(n+1){
        if(init_val) mask.back()<<=((n<<lg_omega) - bits);
        build();
    }
    /// popcount l <= i < r
    uint32_t query(uint32_t l, uint32_t r) const {
        if(__builtin_expect(r < l, false)) return 0;
        return query(r) - query(l);
    }
    /// popcount 0 <= i < x
    uint32_t query(uint32_t x) const {
        uint32_t high = x>>lg_omega, low = x & ((uint64_t{1}<<lg_omega) - 1);
        uint32_t ret = presum_query(high);
        ret+= __builtin_popcountll(mask[high]& ((uint64_t{1} << low)-1));
        return ret;
    }

    void update_pre_build(uint32_t x, bool val){
        uint32_t high = x>>lg_omega, low = x & ((1u<<lg_omega) - 1);
        mask[high] = (mask[high] & ~(uint64_t{1} << low)) | (uint64_t{val}<<low);
    }
    void do_build(){
        build();
    }

    friend ostream& operator<<(ostream&o, Bit_Presum const&b){
        for(auto const&e:b.mask){
            stringstream ss;
            ss << bitset<omega>(e);
            auto s = ss.str();
            reverse(s.begin(), s.end());
            o << s << "|";
        }
        o << " : ";
        for(auto const&e:b.presum) o << e << " ";
        o << "\n";
        return o;
    }

private:
    void presum_build(){
        for(uint32_t x=1;x<=n;++x){
            presum[x]+=presum[x-1];
        }
    }
    // sum 0 <= i < x
    uint32_t presum_query(uint32_t x) const {
        return presum[x];
    }
    void build(){
        for(uint32_t x=0;x<n;++x){
            presum[x+1] = __builtin_popcountll(mask[x]);
        }
        presum_build();
    }

    const uint32_t n;
    vector<uint64_t> mask;
    vector<uint32_t> presum;
};

template<typename T, typename Bit_Ds = Bit_Presum>
class WaveletMatrix{
public:
    static_assert(is_integral<T>::value);
    static constexpr uint32_t height = CHAR_BIT * sizeof(T);

    WaveletMatrix(vector<T> v): n(v.size()), data(height, n){
        build(move(v));
    }
    /// count l <= i < r  s.t.  A <= val[i] < B
    uint32_t range_count(int const&l, int const&r, T const&A, T const&B) const {
        assert(0 <= l && r <= n);
        return count_lower(l, r, B) - count_lower(l, r, A);
    }
    /// count l <= i < r  s.t.  A <= val[i]
    uint32_t range_count_up(int const&l, int const&r, T const&A) const {
        assert(0 <= l && r <= n);
        if(__builtin_expect(l>r, false)) return uint32_t{0};
        return (r-l) - count_lower(l, r, A);
    }
    // k from 0
    T k_th(int const&l, int const&r, int k) const {
        assert(0 <= k && k < n);
        return get_kth(l, r, k);
    }

private:
    void build(vector<T> v){
        m_index.resize(height);
        T const a = numeric_limits<T>::min();
        for(int h = height-1; h>=0;--h){
            T const b = a + (T{1}<<(max(0, h-1))) - !h + (T{1}<<(max(0, h-1)));
            for(int i=0;i<n;++i){
                data[h].update_pre_build(i, v[i]<b);
            }
            data[h].do_build();
            const int m = stable_partition(v.begin(), v.end(), [&b](T const&x){return x < b;}) - v.begin();
            for(int i=m;i<n;++i){
                v[i] = v[i] - (T{1}<<(max(0, h-1))) + !h - (T{1}<<(max(0, h-1)));
            }
            m_index[h] = m;
        }
    }
    /// count l <= i < r  s.t.  val[i] < B
    uint32_t count_lower(int l, int r, T const&B) const {
        assert(0 <= l && r <= n);
        if(__builtin_expect(r<l, false)) return 0;
        uint32_t ret = 0;
        int h = height;
        T a = numeric_limits<T>::min();
        while(h > 0){
            --h;
			bool go_left = B < a + (T{1}<<(max(0, h-1))) - !h + (T{1}<<(max(0, h-1)));
            const int low_l = data[h].query(l), low_r = data[h].query(r);
            if(go_left){
                l = low_l;
                r = low_r;
            } else {
                a = a + (T{1}<<(max(0, h-1))) - !h + (T{1}<<(max(0, h-1)));
                ret+= low_r-low_l;
                l = m_index[h] + l-low_l;
                r = m_index[h] + r-low_r;
            }
        }
        return ret;
    }
    T get_kth(int l, int r, int k) const {
        assert(0 <= l && r <= n);
        assert(0 <= k && k < r-l);
        int h = height;
        T a = numeric_limits<T>::min();
        while(h > 0){
            --h;
            const int low_l = data[h].query(l), low_r = data[h].query(r), low_lr = low_r-low_l;
            bool go_left = k < low_lr;
            if(go_left){
                l = low_l;
                r = low_r;
            } else {
                a+= T{1}<<h;
                k-= low_lr;
                l = m_index[h] + l-low_l;
                r = m_index[h] + r-low_r;
            }
        }
        return a;
    }

    const int n;
    vector<int> m_index;
    vector<Bit_Ds> data;
};
// }}}

## Partially Persistent DSU

struct PersistentDSU{
  vector<int> lab, time_united;
  PersistentDSU(int n):lab(n, -1), time_united(n, INF32){}
 
  int root(int u, int ver){
    if(time_united[u]>ver) return u;
    else return root(lab[u], ver);
  }
 
  void unite(int u, int v, int ver){
    u=root(u, ver), v=root(v, ver);
    if(u!=v){
      if(lab[u]>lab[v]) swap(u, v);
      lab[u]+=lab[v]; lab[v]=u; time_united[v]=ver;
    }
  }
 
  bool same_set(int u, int v, int ver){
    return root(u, ver)==root(v, ver);
  }

# STRING

## Update Hash

struct hash_BIT {
    vector<ll> tr, base;
    ll MOD;

    hash_BIT (int n) : tr(n + 1, 0) {
        MOD = 1e9 + 21;
        base.resize(n + 1);
        base[0] = 1;
        for (int i = 1; i <= n; i++) base[i] = base[i-1] * 29LL % MOD;
        return;
    }

private:
    int p (int k) {
        return k & -k;
    }

    ll pre (int k) {
        ll ans = 0;
        for (int j = k; j > 0; j -= p(j))
            ans = (ans + tr[j] * base[k - j]) % MOD;
        return ans;
    }

public:
    ll get_hash (int l, int r) {
        if (l > r) return 0;
        return (pre(r) - pre(l - 1) * base[r - l + 1] + MOD * MOD) % MOD;
    }

    void update (int k, char c) {
        ll delta = c - 'a' + 1 - get_hash(k, k);
        for (int j = k; j < tr.size(); j += p(j))
            tr[j] = (tr[j] + delta * base[j - k] + MOD * MOD) % MOD;
        return;
    }

    void append (char c) {
        int pos = tr.size();
        tr.push_back((get_hash(pos - p(pos) + 1, pos - 1) * 29LL + c - 'a' + 1) % MOD);
        return;
    }
};

## Manacher

string t; cin >> t;
vector<char> s;
s.push_back('.');
for(char d : t){
    s.push_back(d);
    s.push_back('.');
}
int n = s.size();
vi best(n); // longest palindrome centered at i
int center = 0, lst = 0;
for(int i = 1; i < n; ++i){
    if(lst > i)
    best[i] = min(best[center - (i - center)], lst - i);
    int l = i - best[i], r = i + best[i];
    while(l - 1 >= 0 && r + 1 < n && s[l - 1] == s[r + 1]){
        --l, ++r;
        ++best[i];
    }
    if(i + best[i] > lst){
        lst = i + best[i];
        center = i;
    }
}
ll ans = 0;
for(int i = 0; i < n; ++i){
    best[i] = best[i] * 2 + 1;
    if(s[i] == '.') ans += (best[i] - 1) / 4;
    else ans += (best[i] + 1) / 4;
}
cout << ans;

## Z-function

void z_func(string s, int n){
  vector<int> z(n);
	for(int l = 0, r = 0, i = 1; i < n; ++i){
        if(i <= r) z[i] = min(z[i - l], r - i + 1);
        while(i + z[i] < n && s[z[i]] == s[i + z[i]])
            ++z[i];
        if(i + z[i] - 1 > r){
            l = i;
            r = i + z[i] - 1;
        }
    }
}

## Suffix Array

// Efficient O(N + alphabet_size) time and space suffix array
// For ICPC notebook, it's better to copy a shorter code such as
// https://github.com/kth-competitive-programming/kactl/blob/main/content/strings/SuffixArray.h

// Usage:
// - sa = suffix_array(s, 'a', 'z')
// - lcp = LCP(s, sa)
//   lcp[i] = LCP(sa[i], sa[i+1])
//
// Tested:
// - SA https://judge.yosupo.jp/problem/suffixarray
// - SA https://www.spoj.com/problems/SARRAY/
// - LCP https://judge.yosupo.jp/problem/number_of_substrings
// Suffix Array {{{
// Copied from https://judge.yosupo.jp/submission/52300
// Helper functions {{{
void induced_sort(const std::vector<int>& vec, int val_range,
                  std::vector<int>& SA, const std::vector<bool>& sl,
                  const std::vector<int>& lms_idx) {
    std::vector<int> l(val_range, 0), r(val_range, 0);
    for (int c : vec) {
        if (c + 1 < val_range) ++l[c + 1];
        ++r[c];
    }
    std::partial_sum(l.begin(), l.end(), l.begin());
    std::partial_sum(r.begin(), r.end(), r.begin());
    std::fill(SA.begin(), SA.end(), -1);
    for (int i = (int)lms_idx.size() - 1; i >= 0; --i)
        SA[--r[vec[lms_idx[i]]]] = lms_idx[i];
    for (int i : SA)
        if (i >= 1 && sl[i - 1]) SA[l[vec[i - 1]]++] = i - 1;
    std::fill(r.begin(), r.end(), 0);
    for (int c : vec) ++r[c];
    std::partial_sum(r.begin(), r.end(), r.begin());
    for (int k = (int)SA.size() - 1, i = SA[k]; k >= 1; --k, i = SA[k])
        if (i >= 1 && !sl[i - 1]) {
            SA[--r[vec[i - 1]]] = i - 1;
        }
}

std::vector<int> SA_IS(const std::vector<int>& vec, int val_range) {
    const int n = vec.size();
    std::vector<int> SA(n), lms_idx;
    std::vector<bool> sl(n);
    sl[n - 1] = false;
    for (int i = n - 2; i >= 0; --i) {
        sl[i] = (vec[i] > vec[i + 1] || (vec[i] == vec[i + 1] && sl[i + 1]));
        if (sl[i] && !sl[i + 1]) lms_idx.push_back(i + 1);
    }
    std::reverse(lms_idx.begin(), lms_idx.end());
    induced_sort(vec, val_range, SA, sl, lms_idx);
    std::vector<int> new_lms_idx(lms_idx.size()), lms_vec(lms_idx.size());
    for (int i = 0, k = 0; i < n; ++i)
        if (!sl[SA[i]] && SA[i] >= 1 && sl[SA[i] - 1]) {
            new_lms_idx[k++] = SA[i];
        }
    int cur = 0;
    SA[n - 1] = cur;
    for (size_t k = 1; k < new_lms_idx.size(); ++k) {
        int i = new_lms_idx[k - 1], j = new_lms_idx[k];
        if (vec[i] != vec[j]) {
            SA[j] = ++cur;
            continue;
        }
        bool flag = false;
        for (int a = i + 1, b = j + 1;; ++a, ++b) {
            if (vec[a] != vec[b]) {
                flag = true;
                break;
            }
            if ((!sl[a] && sl[a - 1]) || (!sl[b] && sl[b - 1])) {
                flag = !((!sl[a] && sl[a - 1]) && (!sl[b] && sl[b - 1]));
                break;
            }
        }
        SA[j] = (flag ? ++cur : cur);
    }
    for (size_t i = 0; i < lms_idx.size(); ++i) lms_vec[i] = SA[lms_idx[i]];
    if (cur + 1 < (int)lms_idx.size()) {
        auto lms_SA = SA_IS(lms_vec, cur + 1);
        for (size_t i = 0; i < lms_idx.size(); ++i) {
            new_lms_idx[i] = lms_idx[lms_SA[i]];
        }
    }
    induced_sort(vec, val_range, SA, sl, new_lms_idx);
    return SA;
}
// }}}

template<typename ContainerT = std::string, typename ElemT = unsigned char>
std::vector<int> suffix_array(const ContainerT& s, const ElemT first = 'a',
                         const ElemT last = 'z') {
    std::vector<int> vec(s.size() + 1);
    std::copy(std::begin(s), std::end(s), std::begin(vec));
    for (auto& x : vec) x -= (int)first - 1;
    vec.back() = 0;
    auto ret = SA_IS(vec, (int)last - (int)first + 2);
    ret.erase(ret.begin());
    return ret;
}

std::vector<int> LCP(const std::string& s, const std::vector<int>& sa) {
    int n = s.size(), k = 0;
    std::vector<int> lcp(n), rank(n);
    for (int i = 0; i < n; i++) rank[sa[i]] = i;
    for (int i = 0; i < n; i++, k ? k-- : 0) {
        if (rank[i] == n - 1) {
            k = 0;
            continue;
        }
        int j = sa[rank[i] + 1];
        while (i + k < n && j + k < n && s[i + k] == s[j + k]) k++;
        lcp[rank[i]] = k;
    }
    lcp[n - 1] = 0;
    return lcp;
}

int64_t cnt_distinct_substrings(const std::string& s) {
    auto lcp = LCP(s, suffix_array(s, 0, 255));
    return s.size() * (int64_t) (s.size() + 1) / 2
        - std::accumulate(lcp.begin(), lcp.end(), 0LL);
}

## Prefix Function

vector<int> prefixFunc (string s) {
    int n = s.length();
    vector<int> pi(n);
    for (int i = 1; i < n; i++) {
        int j = pi[i - 1];
        while (j && s[j] != s[i]) j = pi[j - 1];
        if (s[j] == s[i]) pi[i] = j + 1;
    }
    return pi;
}

# GEO

## Basic

template<class T> struct Point{
  T x, y;
  Point(T _x=0, T _y=0):x(_x), y(_y){}

  template<class U> operator Point<U>(){
    return Point<U>(U(x), U(y));
  }

  Point &operator+=(Point p){
    x+=p.x; y+=p.y;
    return *this;
  }

  Point &operator-=(Point p){
    x-=p.x; y-=p.y;
    return *this;
  }

  Point &operator*=(T v){
    x*=v; y*=v;
    return *this;
  }

  Point &operator/=(T v){
    x/=v; y/=v;
    return *this;
  }

  Point operator-() const{
    return Point(-x, -y);
  }

  friend Point operator+(Point a, Point b){
    return a+=b;
  }

  friend Point operator-(Point a, Point b){
    return a-=b;
  }

  friend Point operator*(Point a, T v){
    return a*=v;
  }

  friend Point operator*(T v, Point a){
    return a*=v;
  }

  friend Point operator/(Point a, T v){
    return a/=v;
  }

  friend bool operator==(Point a, Point b){
    return a.x==b.x and a.y==b.y;
  }

  friend istream& operator>>(istream& in, Point& p){
    return in>>p.x>>p.y;
  }

  friend ostream& operator<<(ostream& out, Point p){
    return out<<p.x<<' '<<p.y;
  }
};

template<class T> T dot(Point<T> a, Point<T> b){
  return a.x*b.x+a.y*b.y;
}

template<class T> T cross(Point<T> a, Point<T> b){
  return a.x*b.y-a.y*b.x;
}

template<class T> T dot(Point<T> a, Point<T> b, Point<T> c){
  return dot(b-a, c-a);
}

template<class T> T cross(Point<T> a, Point<T> b, Point<T> c){
  return cross(b-a, c-a);
}

template<class T> T square(Point<T> a){
  return dot(a, a);
}

template<class T> double length(Point<T> p){
  return sqrt(square(p));
}

long double length(Point<long double> p){
  return sqrt(square(p));
}

template<class T> Point<T> normalize(Point<T> p){
  return p/length(p);
}

template<class T> T manhattan(Point<T> a, Point<T> b){
  return abs(a.x-b.x)+abs(a.y-b.y);
}

template<class T> int sign(T x){
  return (x==0?0:(x<0?-1:1));
}

//segment AB, point C
template<class T> bool onSegment(Point<T> a, Point<T> b, Point<T> c){
  return cross(a-c, b-c)==0 and dot(a-c, b-c)<=0;
}

//segment AB, CD
template<class T> bool intersect(Point<T> a, Point<T> b, Point<T> c, Point<T> d){
  return onSegment(a, b, c) or onSegment(a, b, d) or onSegment(c, d, a) or onSegment(c, d, b) or
  (sign(cross(a, b, c))*sign(cross(a, b, d))<0 and sign(cross(c, d, a))*sign(cross(c, d, b))<0);
}

using point=Point<long long>;

template<class T> struct Line{
  T a, b, c; //ax+by+c=0
  Line(Point<T> A, Point<T> B){
    a=B.y-A.y;
    b=A.x-B.x;
    c=-A.x*B.y+A.y*B.x;
  }
  Line(T _a=0, T _b=0, T _c=0){
    a=_a, b=_b, c=_c;
  }
  T func(T x, T y){
    return a*x+b*y+c;
  }
  friend istream& operator>>(istream& in, Line& d){
    return in>>d.a>>d.b>>d.c;
  }
};

//line AB, point C
template<class T> long double LineToPoint(Point<T> a, Point<T> b, Point<T> c){
  return abs(cross(a, b, c))/length(a-b);
}

//line AB, point C
template<class T> long double LineToPoint(Line<T> d, Point<T> p){
  return abs(d.func(p.x, p.y))/sqrt(d.a*d.a+d.b*d.b);
}

using line=Line<long long>;

//line AB, CD
template<class T> bool parallel(Point<T> a, Point<T> b, Point<T> c, Point<T> d){
  Line<T> d1(a, b), d2(c, d);
  return (d1.a*d2.b==d1.b*d2.a);
}

//ray AB, point C
template<class T> long double RayToPoint(Point<T> a, Point<T> b, Point<T> c){
  return (dot(a, b, c)>=0?LineToPoint(a, b, c):min(length(a-c), length(b-c)));
}

## Polygon

// Polygon: convex hull, in polygon, convex diameter ..
// Functions with param vector<P<T>> works with both integers and floating points
// Functions with param Polygon works with floating points only.

typedef vector< Point > Polygon;

// Convex Hull:
// If minimum point --> #define REMOVE_REDUNDANT
// Known issues:
// - Max. point does not work when some points are the same.
// Tested:
// - (min points) https://open.kattis.com/problems/convexhull
// - (max points) https://cses.fi/problemset/task/2195

// #define REMOVE_REDUNDANT
template<typename T>
T area2(P<T> a, P<T> b, P<T> c) { return a%b + b%c + c%a; }

template<typename T>
void ConvexHull(vector<P<T>> &pts) {
    sort(pts.begin(), pts.end());
    pts.erase(unique(pts.begin(), pts.end()), pts.end());
    vector<P<T>> up, dn;
    for (int i = 0; i < (int) pts.size(); i++) {
#ifdef REMOVE_REDUNDANT
        while (up.size() > 1 && area2(up[up.size()-2], up.back(), pts[i]) >= 0) up.pop_back();
        while (dn.size() > 1 && area2(dn[dn.size()-2], dn.back(), pts[i]) <= 0) dn.pop_back();
#else
        while (up.size() > 1 && area2(up[up.size()-2], up.back(), pts[i]) > 0) up.pop_back();
        while (dn.size() > 1 && area2(dn[dn.size()-2], dn.back(), pts[i]) < 0) dn.pop_back();
#endif
        up.push_back(pts[i]);
        dn.push_back(pts[i]);
    }
    pts = dn;
    for (int i = (int) up.size() - 2; i >= 1; i--) pts.push_back(up[i]);
    
#ifdef REMOVE_REDUNDANT
    if (pts.size() <= 2) return;
    dn.clear();
    dn.push_back(pts[0]);
    dn.push_back(pts[1]);
    for (int i = 2; i < (int) pts.size(); i++) {
        if (onSegment(dn[dn.size()-2], pts[i], dn.back())) dn.pop_back();
        dn.push_back(pts[i]);
    }
    if (dn.size() >= 3 && onSegment(dn.back(), dn[1], dn[0])) {
        dn[0] = dn.back();
        dn.pop_back();
    }
    pts = dn;
#endif
}

// Area, perimeter, centroid
template<typename T>
T signed_area2(vector<P<T>> p) {
    T area(0);
    for(int i = 0; i < (int) p.size(); i++) {
        area += p[i] % p[(i + 1) % p.size()];
    }
    return area;
}
double area(const Polygon &p) {
    return std::abs(signed_area2(p) / 2.0);
}
Point centroid(Polygon p) {
    Point c(0,0);
    double scale = 3.0 * signed_area2(p);
    for (int i = 0; i < (int) p.size(); i++){
        int j = (i+1) % p.size();
        c = c + (p[i]+p[j])*(p[i].x*p[j].y - p[j].x*p[i].y);
    }
    return c / scale;
}
double perimeter(Polygon p) {
    double res = 0;
    for(int i = 0; i < (int) p.size(); ++i) {
        int j = (i + 1) % p.size();
        res += (p[i] - p[j]).len();
    }
    return res;
}

// Is convex: checks if polygon is convex.
// Return true for degen points (all vertices are collinear)
template<typename T>
bool is_convex(const vector<P<T>> &P) {
    int sz = (int) P.size();
    int min_ccw = 2, max_ccw = -2;
    for (int i = 0; i < sz; i++) {
        int c = ccw(P[i], P[(i+1) % sz], P[(i+2) % sz]);
        min_ccw = min(min_ccw, c);
        max_ccw = max(max_ccw, c);
    }
    return min_ccw * max_ccw >= 0;
}

// Inside polygon: O(N). Works with any polygon
// Tested:
// - https://open.kattis.com/problems/pointinpolygon
// - https://open.kattis.com/problems/cuttingpolygon
enum PolygonLocation { OUT, ON, IN };
PolygonLocation in_polygon(const Polygon &p, Point q) {
    if ((int)p.size() == 0) return PolygonLocation::OUT;

    // Check if point is on edge.
    int n = p.size();
    REP(i,n) {
        int j = (i + 1) % n;
        Point u = p[i], v = p[j];

        if (u > v) swap(u, v);

        if (ccw(u, v, q) == 0 && u <= q && q <= v) return PolygonLocation::ON;
    }

    // Check if point is strictly inside.
    int c = 0;
    for (int i = 0; i < n; i++) {
        int j = (i + 1) % n;
        if (((p[i].y <= q.y && q.y < p[j].y)
                    || (p[j].y <= q.y && q.y < p[i].y))
                && q.x < p[i].x + (p[j].x - p[i].x) * (q.y - p[i].y) / (double) (p[j].y - p[i].y)) {
            c = !c;
        }
    }
    return c ? PolygonLocation::IN : PolygonLocation::OUT;
}

// Check point in convex polygon, O(logN)
#define Det(a,b,c) ((double)(b.x-a.x)*(double)(c.y-a.y)-(double)(b.y-a.y)*(c.x-a.x))
PolygonLocation in_convex(vector<Point>& l, Point p){
    if (l.empty()) return PolygonLocation::OUT;
    if (l.size() <= 2) {
        return onSegment(l[0], l[1 % l.size()], p) ? PolygonLocation::ON : PolygonLocation::OUT;
    }

    int a = 1, b = l.size()-1, c;
    if (Det(l[0], l[a], l[b]) > 0) swap(a,b);

    if (onSegment(l[0], l[a], p)) return ON;
    if (onSegment(l[0], l[b], p)) return ON;

    if (Det(l[0], l[a], p) > 0 || Det(l[0], l[b], p) < 0) return OUT;
    while(abs(a-b) > 1) {
        c = (a+b)/2;
        if (Det(l[0], l[c], p) > 0) b = c; else a = c;
    }
    int t = cmp(Det(l[a], l[b], p), 0);
    return (t == 0) ? ON : (t < 0) ? IN : OUT;
}


// Cut a polygon with a line. Returns half on left-hand-side.
// To return the other half, reverse the direction of Line l (by negating l.a, l.b)
// The line must be formed using 2 points
Polygon polygon_cut(const Polygon& P, Line l) {
    Polygon Q;
    for(int i = 0; i < (int) P.size(); ++i) {
        Point A = P[i], B = (i == ((int) P.size())-1) ? P[0] : P[i+1];
        if (ccw(l.A, l.B, A) != -1) Q.push_back(A);
        if (ccw(l.A, l.B, A)*ccw(l.A, l.B, B) < 0) {
            Point p; areIntersect(Line(A, B), l, p);
            Q.push_back(p);
        }
    }
    return Q;
}

// Find intersection of 2 convex polygons
// Helper method
bool intersect_1pt(Point a, Point b,
    Point c, Point d, Point &r) {
    double D =  (b - a) % (d - c);
    if (cmp(D, 0) == 0) return false;
    double t =  ((c - a) % (d - c)) / D;
    double s = -((a - c) % (b - a)) / D;
    r = a + (b - a) * t;
    return cmp(t, 0) >= 0 && cmp(t, 1) <= 0 && cmp(s, 0) >= 0 && cmp(s, 1) <= 0;
}
Polygon convex_intersect(Polygon P, Polygon Q) {
    const int n = P.size(), m = Q.size();
    int a = 0, b = 0, aa = 0, ba = 0;
    enum { Pin, Qin, Unknown } in = Unknown;
    Polygon R;
    do {
        int a1 = (a+n-1) % n, b1 = (b+m-1) % m;
        double C = (P[a] - P[a1]) % (Q[b] - Q[b1]);
        double A = (P[a1] - Q[b]) % (P[a] - Q[b]);
        double B = (Q[b1] - P[a]) % (Q[b] - P[a]);
        Point r;
        if (intersect_1pt(P[a1], P[a], Q[b1], Q[b], r)) {
            if (in == Unknown) aa = ba = 0;
            R.push_back( r );
            in = B > 0 ? Pin : A > 0 ? Qin : in;
        }
        if (C == 0 && B == 0 && A == 0) {
            if (in == Pin) { b = (b + 1) % m; ++ba; }
            else            { a = (a + 1) % m; ++aa; }
        } else if (C >= 0) {
            if (A > 0) { if (in == Pin) R.push_back(P[a]); a = (a+1)%n; ++aa; }
            else       { if (in == Qin) R.push_back(Q[b]); b = (b+1)%m; ++ba; }
        } else {
            if (B > 0) { if (in == Qin) R.push_back(Q[b]); b = (b+1)%m; ++ba; }
            else       { if (in == Pin) R.push_back(P[a]); a = (a+1)%n; ++aa; }
        }
    } while ( (aa < n || ba < m) && aa < 2*n && ba < 2*m );
    if (in == Unknown) {
        if (in_convex(Q, P[0])) return P;
        if (in_convex(P, Q[0])) return Q;
    }
    return R;
}

// Find the diameter of polygon.
// Returns:
// <diameter, <ids of 2 points>>
pair<double, pair<int,int>> convex_diameter(const Polygon &p) {
    const int n = p.size();
    int is = 0, js = 0;
    for (int i = 1; i < n; ++i) {
        if (p[i].y > p[is].y) is = i;
        if (p[i].y < p[js].y) js = i;
    }
    double maxd = (p[is]-p[js]).len();
    int i, maxi, j, maxj;
    i = maxi = is;
    j = maxj = js;
    do {
        int ii = (i+1) % n, jj = (j+1) % n;
        if ((p[ii] - p[i]) % (p[jj] - p[j]) >= 0) j = jj;
        else i = ii;
        if ((p[i] - p[j]).len() > maxd) {
            maxd = (p[i] - p[j]).len();
            maxi = i;
            maxj = j;
        }
    } while (i != is || j != js);
    return {maxd, std::minmax(maxi, maxj)}; /* farthest pair is (maxi, maxj). */
}

// Pick theorem
// Given non-intersecting polygon.
// S = area
// I = number of integer points strictly Inside
// B = number of points on sides of polygon
// S = I + B/2 - 1