Submission #7571733

---

Source Code Expand

#include <bits/stdc++.h>
using namespace std;

using i8  = int8_t;using i16 = int16_t;using i32 = int32_t;using i64 = int64_t;using u8  = uint8_t;using u16 = uint16_t;using u32 = uint32_t;using u64 = uint64_t;using f32 = float;using f64 = double;
using vi32 = vector<i32>;using vi64 = vector<i64>;using vvi32 = vector<vi32>;using vvi64 = vector<vi64>;using pi32 = pair<i32,i32>;using pi64 = pair<i64,i64>;

#define FOR(i,a,b) for(i64 i=(a), i##_len=(b); i<i##_len; ++i)
#define REP(i,n) FOR(i,0,n)
#define REPS(i,x) for(i64 i=1;i<=(i64)(x);i++)
#define RREPS(i,x) for(i64 i=((i64)(x));i>0;i--)
#define ALL(obj) (obj).begin(),(obj).end()
#define SZ(x) ((i32)(x).size())
#define CLR(ar,val) memset(ar, val, sizeof(ar))
#define cauto const auto&
#define pb push_back
#define mp make_pair
#define UNIQUE(v) v.erase( unique(v.begin(), v.end()), v.end() );
#define perm(c) sort(ALL(c));for(bool c##p=1;c##p;c##p=next_permutation(all(c)))

i32 dx[4]={1,0,-1,0};
i32 dy[4]={0,1,0,-1};
const i32 INF32 = 0x3F3F3F3F;
const i64 INF64 = 0x3F3F3F3F3F3F3F3F;

template<class T> inline bool chmax(T &a, const T &b) { if (a<b) { a=b; return 1; } return 0; }
template<class T> inline bool chmin(T &a, const T &b) { if (b<a) { a=b; return 1; } return 0; }

/// pair vectorに対して直接入出力できるようにする
template<typename A, typename B> istream &operator>>(istream &is, pair<A, B> &p) { is >> p.first >> p.second;return is; }
template<typename A, typename B> ostream &operator<<(ostream &os, const pair<A, B>& p) { os << p.first << ' ' << p.second;return os; }
template<typename T> istream &operator>>(istream &is, vector<T> &v) { for(T& in : v) is >> in; return is; }
template<typename T> ostream &operator<<(ostream &os, const vector<T> &v) { for(i32 i = 0; i < (i32)v.size(); i++) os << v[i] << (i+1 != (i32)v.size() ? " " : ""); return os; }



/// # べき乗
i64 mod_pow(i64 n, i64 k, i64 mod) {
    i64 res = 1;
    for(; k>0; k>>=1) {
        if(k & 1) (res *= n) %= mod;
        (n *= n) %= mod;
    }
    return res;
}


/// # 余りで四則演算する構造体
template<i64 mod>
struct ModInt {
    i64 v;
    i64 mod_pow(i64 x, i64 n) const {
        return (!n) ? 1 : (mod_pow((x*x)%mod,n/2) * ((n&1)?x:1)) % mod;
    }
    ModInt(i64 a = 0) : v(a >= mod ? a % mod : a) {}
    ModInt operator+ ( const ModInt& b ) const {
        return (v + b.v >= mod ? ModInt(v + b.v - mod) : ModInt(v + b.v));
    }
    ModInt operator- () const {
        return ModInt(-v);
    }
    ModInt operator- ( const ModInt& b ) const {
        return (v - b.v < 0 ? ModInt(v - b.v + mod) : ModInt(v - b.v));
    }
    ModInt operator* ( const ModInt& b ) const {return (v * b.v) % mod;}
    ModInt operator/ ( const ModInt& b ) const {return (v * mod_pow(b.v, mod-2)) % mod;}
    
    bool operator== ( const ModInt &b ) const {return v == b.v;}
    ModInt& operator+= ( const ModInt &b ) {
        v += b.v;
        if(v >= mod) v -= mod;
        return *this;
    }
    ModInt& operator-= ( const ModInt &b ) {
        v -= b.v;
        if(v < 0) v += mod;
        return *this;
    }
    ModInt& operator*= ( const ModInt &b ) {
        (v *= b.v) %= mod;
        return *this;
    }
    ModInt& operator/= ( const ModInt &b ) {
        (v *= mod_pow(b, mod-2)) %= mod;
        return *this;
    }
    ModInt pow(i64 x) { return ModInt(mod_pow(v, x)); }
    operator i32() const { return int(v); }
    operator i64() const { return v; }
};
 
template<i64 mod>
ostream& operator<< (ostream& out, ModInt<mod> a) {return out << a.v;}
template<i64 mod>
istream& operator>> (istream& in, ModInt<mod>& a) {
    in >> a.v;
    return in;
}

template <typename NumType>
struct Combination {
    i32 LIMIT;
    vector<NumType> fact_, finv_;

    Combination() {}
    Combination(i32 LIMIT_) : LIMIT(LIMIT_), fact_(LIMIT+1), finv_(LIMIT+1) {
        fact_[0] = finv_[LIMIT] = NumType(1);
        for(i32 i=1; i<=LIMIT; i++) {
            fact_[i] = fact_[i-1] * NumType(i);
        }
        
        finv_[LIMIT] /= fact_[LIMIT];
        for(i32 i=LIMIT-1; i>=0; i--) {
            finv_[i] = finv_[i+1] * NumType(i+1);
        }
    }

    inline NumType fact(i32 k) const { return fact_[k]; }
    inline NumType finv(i32 k) const { return finv_[k]; }
    NumType P(i32 n, i32 r) const {
        if(r < 0 or n < r) return NumType(0);
        return fact_[n] * finv_[n-r];
    }
    NumType C(i32 n, i32 r) const {
        if(r < 0 or n < r) return NumType(0);
        return fact_[n] * finv_[n-r] * finv_[r];
    }
    // 重複組み合わせ
    NumType H(i32 n, i32 r) const {
        if(n < 0 or r < 0) return NumType(0);
        return r == 0 ? NumType(0) : C(n + r - 1, r);
    }
    // ベル数 (区別できる n 個のボールを区別できない k 個以下の箱に分割)
    // B(n, n) := n 個のボールを任意個のグループに分割する場合の数
    NumType B(i32 n, i32 k) const {
        if(n == 0) return NumType(1);
        k = min(n, k);
        NumType ret(0);
        vector<NumType> pref(k + 1); pref[0] = NumType(1);
        for(i32 i=1; i<=k; i++) {
            if(i & 1) pref[i] = pref[i-1] - finv_[i];
            else pref[i] = pref[i-1] + finv_[i];
        }
        for(i32 i=1; i<=k; i++) {
            ret += NumType(i).pow(n) * finv_[i] * pref[k-i];
        }
        return ret;
    }
    // スターリング数 (区別できる n 個のボールを区別できない k 個の箱に分割)
    NumType S(i32 n, i32 k) const {
        if(n < k) return NumType(0);
        NumType ans(0);
        for(i32 i=0; i<=k; i++) {
            NumType val = C(k, i) * NumType(i).pow(n);
            if((k - i) % 2) ans -= val;
            else ans += val;
        }
        return ans * finv_[k];
    }
};

template <typename NumType>
vector<NumType> get_inv_table(i32 N, i32 P) {
    vector<NumType> res;
    for(i32 i=0; i<N; i++) {
        res.emplace_back(NumType(mod_pow(i, P-2, P)));
    }
    return res;
}

// f(0) ... f(N) の値がわかっているときに
// 多項式補間を O(N) で実現するバージョンのラグランジュ補間
template <typename NumType, i32 LIMIT = 2000010>
struct SpecificLagrangePolynomial {
    Combination<NumType> comb;
    vector<NumType> y, inv_table;
    bool use_inv_table;
    SpecificLagrangePolynomial() : comb(), y() {}
    SpecificLagrangePolynomial(vector<NumType> y_, i32 P=LIMIT)
        : comb(min(LIMIT, P-1)), y(y_), use_inv_table(false) {
        if(P < LIMIT) use_inv_table = true;
        
        if(use_inv_table) {
            inv_table = get_inv_table<NumType>(P, P);
        }
    }

    NumType interpolate(NumType p) {
        NumType pre(1);
        i32 N = y.size(), pv = i32(p);
        if(pv < N) return y[pv];
        for(i32 i=0; i<N; i++) {
            pre *= NumType(pv - i);
        }
        
        NumType res(0);
        for(i32 i=0; i<N; i++) {
            NumType numer(0), denom(1);
            if(use_inv_table)
                numer = pre * inv_table[pv - i];
            else
                numer = pre / NumType(pv - i);

            denom *= comb.finv(N - 1 - i) * comb.finv(i);
            if((N - 1 - i) % 2) res -= numer * denom * y[i];
            else res += numer * denom * y[i];
        }
        return res;
    }
};


/// # 組み合わせ
/// nCr をmodで割った余りを求める
i64 CombinationMod(i32 n, i32 r, i32 mod=1) {
    if (n < 0 || r < 0 || r > n) return -1;
    if (n - r < r) r = n - r;
    if (r == 0) return 1;
    if (r == 1) return n;
    vi64 numerator(r);
    vi64 denominator(r);

    for (i32 k = 0; k < r; k++) {
        numerator[k] = n - r + k + 1;
        denominator[k] = k + 1;
    }

    for (i32 p = 2; p <= r; p++) {
        i32 pivot = denominator[p-1];
        if (pivot > 1) {
            i32 offset = (n - r) % p;
            for (i32 k = p - 1; k < r; k += p) {
                numerator[k-offset] /= pivot;
                denominator[k] /= pivot;
            }
        }
    }
    i64 ret = 1LL;
    if (mod == 1) {
        for (i32 k = 0; k < r; k++) {
            if (numerator[k] > 1) ret *= numerator[k];
        }
    } else {
        for (i32 k = 0; k < r; k++) {
            if (numerator[k] > 1) ret = (ret * numerator[k]) % mod;
        }
    }
    return ret;
}


/// # 拡張ユークリッドの互除法
/// gcd(a, b) = a * x + b * yのxとyを求める
template<typename T>
T extgcd(T a, T b, T &x, T &y) {
    T d = a;
    if(b != 0) {
        d = extgcd(b, a % b, y, x);
        y -= (a / b) * x;
    } else {
        x = 1;
        y = 0;
    }
    return d;
}


/// # オイラーのφ
i64 euler_phi(i64 n) {
    i64 ret = n;
    for (i64 i = 2; i * i <= n; i++) {
        if (n % i == 0) {
            ret -= ret / i;
            while (n % i == 0) n /= i;
        }
    }
    if (n > 1) ret -= ret / n;
    return ret;
}
vi32 euler_phi_table(i32 n) {
    vi32 euler(n + 1);
    for (i32 i = 0; i <= n; i++) {
        euler[i] = i;
    }
    for (i32 i = 2; i <= n; i++) {
        if (euler[i] == i) {
            for (i32 j = i; j <= n; j += i) {
                euler[j] = euler[j] / i * (i - 1);
            }
        }
    }
    return euler;
}

/// # 進数変換
template<typename T>
vector<T> convert_base(T x, T b) {
    vector<T> ret;
    T t = 1, k = abs(b);
    while(x) {
        ret.emplace_back((x * t) % k);
        if(ret.back() < 0) ret.back() += k;
        x -= ret.back() * t;
        x /= k;
        t *= b / k;
    }
    if(ret.empty()) ret.emplace_back(0);
    reverse(ALL(ret));
    return ret;
}


/// # 素数テーブル
vector<bool> prime_table(i32 n) {
    vector<bool> prime(n + 1, true);
    if(n >= 0) prime[0] = false;
    if(n >= 1) prime[1] = false;
    for(i32 i = 2; i * i <= n; i++) {
        if(!prime[i]) continue;
        for(i32 j = i + i; j <= n; j += i) {
            prime[j] = false;
        }
    }
    return prime;
}


/// # 素数判定
bool is_prime(i64 x) {
  for(i64 i = 2; i * i <= x; i++) {
    if(x % i == 0) return false;
  }
  return true;
}

/// # 離散対数
/// mod_log(a, b, p) : a^x === b (mod p) を満たすx(>=0)の最小値 存在しない場合-1
i64 mod_log(i64 a, i64 b, i64 p) {
    i64 ok = p, ng = -1;
    while(ok - ng > 1) {
        auto mid = (ok + ng) / 2;
        if(mid * mid >= p) ok = mid;
        else ng = mid;
    }
    
    unordered_map<i64, i64> baby;
    baby.reserve(ok);
    i64 factor = 1;
    for(i64 i = 0, e = b; i < ok; i++) {
        baby[e] = i;
        (factor *= a) %= p;
        (e *= a) %= p;
    }
    for(i64 i = 1, e = factor; i <= ok; i++) {
        auto it = baby.find(e);
        if(it != end(baby)) return i * ok - it->second;
        (e *= factor) %= p;
    }
    return -1;
}



/// # Nの商
/// Nをi (1 <= i <= N)で割ったとき商をすべて列挙する
/// ```c++
///   auto a = quotient_range(8);
///   for (const auto& e : a) cout << e << endl;
///   /***
///     1 1 8
///     2 2 4
///     3 4 2
///     5 8 1
///   ***/
/// ```
template<typename T>
vector<pair<pair<T, T>, T>> quotient_range(T N) {
    T M;
    vector<pair<pair<T, T>, T>> ret;
    for (M = 1; M * M <= N; M++) {
        ret.emplace_back(make_pair(M, M), N / M);
    }
    for (T i = M; i >= 1; i--) {
        T L = N / (i+1) + 1;
        T R = N / i;
        if (L <= R && ret.back().first.second < L) ret.emplace_back(make_pair(L, R), N / L);
    }
    return ret;
}


/// # 素因数分解
map<i64, i32> prime_factor(i64 n) {
    map<i64, i32> ret;
    for (i64 i = 2; i * i <= n; i++) {
        while(n % i == 0) {
            ret[i]++;
            n /= i;
        }
    }
    if (n != 1) ret[n] = 1;
    return ret;
}


/// # 約数
vi64 divisor(i64 n) {
    vi64 ret;
    for (i64 i = 1; i * i <= n; i++) {
        if (n % i == 0) {
            ret.push_back(i);
            if (i * i != n) ret.push_back(n / i);
        }
    }

    /* sort(ALL(ret)); */
    return ret;
}

/// # SegmentTree
/// モノイドに対する区間への演算がO(logN)でできる
///
/// ## 初期化
/// サイズnで初期化 fは2つの区間の要素をマージする演算 M1はモノイドの単位元
/// ```c++:segment-tree
///   SegmentTree(n,f,M1)
/// ```
///
/// ## 代入
/// k番目の要素にxを代入する
/// ```c++:segment-tree
///   set(k, x)
/// ```
///
/// ## 構築
/// セグメント木を構築する
/// ```c++:segment-tree
///   build()
/// ```
///
/// ## 演算
/// 区間[a,b)に二項演算した結果を返す
/// ```c++:segment-tree
///   query(a, b)
/// ```
///
/// ## 変更
/// k番目の要素をxに変更する
/// ```c++:segment-tree
///   update(k, x)
/// ```
///
/// ## 要素取得
/// k番目の要素を取得する
/// ```c++:segment-tree
///   SegmentTree<i32> seg(30, [](i32 a, i32 b) { return min(a, b); }, INT_MAX);
///   seg[10];
/// ```
template<typename Monoid>
struct SegmentTree {
    using F = function<Monoid(Monoid, Monoid)>;
    i32 sz;
    vector<Monoid> seg;
    const F f;
    const Monoid M1;
    SegmentTree(i32 n, const F f, const Monoid &M1) : f(f), M1(M1) {
        sz = 1;
        while(sz < n) sz <<= 1;
        seg.assign(2 * sz, M1);
    }
    void set(i32 k, const Monoid &x) {
        seg[k+sz] = x;
    }
    void build() {
        for(i32 k = sz - 1; k > 0; k--) {
            seg[k] = f(seg[2 * k + 0], seg[2 * k + 1]);
        }
    }
    void update(i32 k, const Monoid &x) {
        k += sz;
        seg[k] = x;
        while(k >>= 1) {
            seg[k] = f(seg[2 * k + 0], seg[2 * k + 1]);
        }
    }
    Monoid query(i32 a, i32 b) {
        Monoid L = M1, R = M1;
        for(a += sz, b += sz; a < b; a >>= 1, b >>= 1) {
            if (a & 1) L = f(L, seg[a++]);
            if (b & 1) R = f(seg[--b], R);
        }
        return f(L, R);
    }
    Monoid operator[](const i32 &k) const {
        return seg[k+sz];
    }
};

/// # BinaryIndexedTree
/// 数列にある要素を加える操作と区間和を求める操作をO(logN)でできるデータ構造
/// 
/// ## 初期化
/// ```c++:BIT
///   BinaryIndexedTree(100)
/// ```
/// 
/// ## 閉区間[0,k]の合計
/// ```c++:BIT
///   sum(k)
/// ```
/// 
/// ## k番目の数列の要素にxを加える
/// ```c++:BIT
///   add(k, x)
/// ```
template<typename T>
struct BinaryIndexedTree {
    vector<T> dat;
    BinaryIndexedTree(i32 sz) {
        dat.assign(++sz, 0);
    }
    T sum(i32 k) {
        T ret = 0;
        for(++k; k > 0; k -= k & -k) ret += dat[k];
        return ret;
    }
    void add(i32 k, T x) {
        for(++k; k < dat.size(); k += k & -k) dat[k] += x;
    }
};

/// # Trie木 [wiki/トライ木](https://ja.wikipedia.org/wiki/%E3%83%88%E3%83%A9%E3%82%A4%E6%9C%A8)
/// あるノードの配下のすべてのノードには
/// 自身に対応する文字列に共通するプレフィックスがあり
/// ルートには空の文字列が対応している
///
/// Trie<char_size, margin>() : 文字列の種類数char_size, 開始文字がmarginのトライ木を作成する
/// add(S)     : 文字列Sをトライ木に追加する
/// query(S,f) : 文字列Sのプレフィックスに一致する文字列を検索する。一致した文字列ごとに関数fが呼び出される
/// size()     : ノード数を返す
/// count()    : 存在する文字列の個数を返す
template<i32 char_size>
struct TrieNode {
    i32 nxt[char_size];
    i32 exist;
    vi32 accept;
    TrieNode() : exist(0) {
        CLR(nxt,-1);
    }
};
template<i32 char_size, i32 margin>
struct Trie {
    using Node = TrieNode<char_size>;
    vector<Node> nodes;
    i32 root;

    Trie() : root(0) {
        nodes.push_back(Node());
    }

    void update_direct(i32 node, i32 id) {
        nodes[node].accept.push_back(id);
    }

    void update_child(i32 node, i32 child, i32 id) {
        ++nodes[node].exist;
    }

    void add(const string& str, i32 str_index, i32 node_index, i32 id) {
        if (str_index == SZ(str)) {
            update_direct(node_index, id);
        } else {
            const i32 c = str[str_index] - margin;
            if (nodes[node_index].nxt[c] == -1) {
                nodes[node_index].nxt[c] = (i32)nodes.size();
                nodes.push_back(Node());
            }
            add(str, str_index + 1, nodes[node_index].nxt[c], id);
            update_child(node_index, nodes[node_index].nxt[c], id);
        }
    }

    void add(const string &str, i32 id) {
        add(str, 0, 0, id);
    }

    void add(const string &str) {
        add(str, nodes[0].exist);
    }

    void query(const string &str, const function<void(i32)> &f, i32 str_index, i32 node_index) {
        for (auto &idx : nodes[node_index].accept) f(idx);
        if (str_index == SZ(str)) {
            return;
        } else {
            const i32 c = str[str_index] - margin;
            if (nodes[node_index].nxt[c] == -1) return;
            query(str, f, str_index + 1, nodes[node_index].nxt[c]);
        }
    }

    void query(const string &str, const function<void(i32)> &f) {
        query(str, f, 0, 0);
    }

    i32 count() const {
        return (nodes[0].exist);
    }

    i32 size() const {
        return ((i32)nodes.size());
    }
};


/// # Binary-Trie木
template<typename T>
struct BinaryTrieNode {
    using Node = BinaryTrieNode<T>;

    BinaryTrieNode<T> *nxt[2];
    i32 max_index;

    BinaryTrieNode() : max_index(-1) {
        nxt[0] = nxt[1] = nullptr;
    }

    void update_direct(i32 id) {
        max_index = max(max_index, id);
    }

    void update_child(Node *child, i32 id) {
        max_index = max(max_index, id);
    }

    Node* add(const T &bit, i32 bit_index, i32 id, bool need=true) {
        Node* node = need ? new Node(*this) : this;
        if (bit_index == -1) {
            node->update_child(id);
        } else {
            const i32 c = (bit >> bit_index) & 1;
            if (node->nxt[c] == nullptr) node->nxt[c] = new Node(), need = false;
            node->nxt[c] = node->nxt[c]->add(bit, bit_index - 1, id, need);
            node->update_child(node->nxt[c], id);
        }
        return node;
    }

    inline T min_query(T bit, i32 bit_index, i32 bit2, i32 l) {
        if (bit_index == -1) return bit;
        i32 c = (bit2 >> bit_index) & 1;
        if (nxt[c] != nullptr && l <= nxt[c]->max_index) {
            return nxt[c]->min_query(bit, bit_index - 1, bit2, l);
        } else {
            return nxt[1^c]->min_query(bit | (1LL << bit_index), bit_index - 1, bit2, l);
        }
    }
};

/// # 永続Binary-Trie木
template<typename T, i32 MAX_LOG>
struct PersistentBinaryTrie {
    using Node = BinaryTrieNode<T>;
    Node* root;

    PersistentBinaryTrie(Node* root) : root(root) { }
    PersistentBinaryTrie() : root(new Node()) { }
    PersistentBinaryTrie add(const T &bit, i32 id) {
        return PersistentBinaryTrie(root->add(bit, MAX_LOG, id));
    }
    T min_query(i32 bit, i32 l) {
        return root->min_query(0, MAX_LOG, bit, l);
    }
};


/// # 配列のk個ずつの最小値を求める
///
/// ```c++:slide minimum
///   vector<i32> A = {0,1,2,3,2,1,2,3,4,3,2,1,0};
///   auto B = slideMin(A, 4);
///   //B == {0,1,1,1,1,1,2,2,1,0}
/// ```
template<typename T>
vector<T> slideMin(const vector<T>& v, i32 k) {
    deque<i32> deq;
    vector<T> ret;
    for(i32 i=0;i < (i32)v.size(); i++) {
        while(!deq.empty() && v[deq.back()] >= v[i]) {
            deq.pop_back();
        }
        deq.push_back(i);
        if (i - k + 1 >= 0) {
            ret.emplace_back(v[deq.front()]);
            if (deq.front() == i - k + 1) deq.pop_front();
        }
    }
    return ret;
}

/// # 累積和を計算する
///
/// ## 初期化
/// ```c++
///   CumulativeSum<i64> A(sz);
/// ```
///
/// ## インデックスkに値vを加える
/// ```c++
///   A.add(k,v);
/// ```
///
/// ## 累積和を構築
/// ```c++
///   A.build();
/// ```
///
/// ## 配列を表示
/// ```c++
///   A.printf();
/// ```
///
/// ## 閉区間[0,k]の和を計算
/// ```c++
///   A.query(k);
/// ```
template<class T>
struct CumulativeSum {
    vector<T> dat;
    CumulativeSum(i32 sz) : dat(sz, 0) { };

    void add(i32 k, T x) {
        dat[k] += x;
    }

    void build() {
        for (i32 i = 1; i < (i32)dat.size(); i++) {
            dat[i] += dat[i-1];
        }
    }

    void printf() {
        for (i32 i = 0; i < (i32)dat.size(); i++) {
            if (i == 0) {
                cout << "[" << dat[i];
            } else if (i == ((i32)dat.size() - 1)) {
                cout << ", " << dat[i] << "]" << endl;
            } else {
                cout << ", " << dat[i];
            }
        }
    }

    T query(i64 k) {
        if (k < 0) return 0;
        return (dat[min(k, (i64)(dat.size() - 1))]);
    }
};
/// 累積和の2次元版
template<class T>
struct CumulativeSum2D {
    vector<vector<T>> dat;
    CumulativeSum2D(i32 w, i32 h) : dat(w + 1, vector<T>(h + 1, 0)) { }

    void add(i32 x, i32 y, T z) {
        ++x; ++y;
        if (x >= SZ(dat) || y >= (SZ(dat[0]))) return;
        dat[x][y] += z;
    }

    void build() {
        for (i32 i = 1; i < SZ(dat); i++) {
            for (i32 j = 1; j < SZ(dat[i]); j++) {
                dat[i][j] += dat[i][j-1] + dat[i - 1][j] - dat[i - 1][j - 1];
            }
        }
    }

    void printf() {
        for (i32 i = 0; i < SZ(dat); i++) {
            for (i32 j = 0; j < SZ(dat[i]); j++) {
                cout << dat[i][j] << (j == (SZ(dat[i]) - 1) ? "" : " ");
            }
            cout << endl;
        }
    }

    T query(i32 sx, i32 sy, i32 gx, i32 gy) {
        return (dat[gx][gy] - dat[sx][gy] - dat[gx][sy] + dat[sx][sy]);
    }
};

/// # 最長増加部分列 (LIS)
/// a      : 数列aの最長増加部分列の長さを求める
/// strict : trueのとき狭義単調増加列 falseのとき広義単調増加列
template<typename T>
size_t longest_increasing_subsequence(const vector<T> &a, bool strict) {
    vector<T> lis;
    for (auto &p : a){
        typename vector<T>::iterator it;
        if (strict) it = lower_bound(begin(lis), end(lis), p);
        else it = upper_bound(begin(lis), end(lis), p);
        if (end(lis) == it) lis.emplace_back(p);
        else *it = p;
    }
    return lis.size();
}

/// # 個数制限付きナップザック問題 (O(NW))
/// w : i番目の重さの配列
/// m : i番目の品物をいくつまで選べるかの配列
/// v : i番目の品物の価値の配列
/// W : ナップザックの容量
/// NG: 容量として不適当な値 (-1や-INFなど、初期化のため)
template<typename T, typename Compare = greater<T>>
vector<T> knapsack_limitations(const vi32 &w, const vi32 &m, const vector<T> &v, const i32 &W, const T &NG, const Compare &comp = Compare()) {
    const i32 N = (i32) w.size();
    vector<T> dp(W+1, NG), deqv(W+1);
    dp[0] = T();
    vector<i32> deq(W+1);
    for(i32 i = 0; i < N; i++) {
        for (i32 a = 0; a < w[i]; a++) {
            i32 s = 0, t = 0;
            for (i32 j = 0; w[i] * j + a <= W; j++) {
                if (dp[w[i] * j + a] != NG) {
                    auto val = dp[w[i] * j + a] - j * v[i];
                    while(s < t && comp(val, deqv[t - 1])) --t;
                    deq[t] = j;
                    deqv[t++] = val;
                }
                if (s < t) {
                    dp[j*w[i]+a] = deqv[s] + j * v[i];
                    if (deq[s] == j - m[i]) ++s;
                }
            }
        }
    }
    return dp;    
}

/// # 個数制限付きナップザック問題 (O(N^2max(vi)^2))
/// w : i番目の重さの配列
/// m : i番目の品物をいくつまで選べるかの配列
/// v : i番目の品物の価値の配列
/// W : ナップザックの容量
template<typename T>
T knapsack_limitations2(const vector<T> &w, const vector<T> &m, const vector<i32> &v, const T &W) {
    const i32 N = (i32) w.size();
    auto v_max = *max_element(begin(v), end(v));
    if(v_max == 0) return 0;
    vector<i32> ma(N);
    vector<T> mb(N);
    for(i32 i = 0; i < N; i++) {
        ma[i] = min< T >(m[i], v_max - 1);
        mb[i] = m[i] - ma[i];
    }
    T sum = 0;
    for(i32 i = 0; i < N; i++) sum += ma[i] * v[i];
    auto dp = knapsack_limitations(v, ma, w, sum, T(-1), less<>());
    vector<i32> ord(N);
    iota(begin(ord), end(ord), 0);
    sort(begin(ord), end(ord), [&](i32 a, i32 b) { 
            return v[a] * w[b] > v[b] * w[a];
    });
    T ret = T();
    for(i32 i = 0; i < SZ(dp); i++) {
        if(dp[i] > W || dp[i] == -1) continue;
        T rest = W - dp[i], cost = i;
        for(auto &p : ord) {
            auto get = min(mb[p], rest / w[p]);
            if(get == 0) break;
            cost += get * v[p];
            rest -= get * w[p];
        }
        ret = max(ret, cost);
    }
    return ret;
}


/// # 数論変換NTT
namespace NTT {
i64 extgcd(i64 a, i64 b, i64 &x, i64 &y) {
    i64 d;
    return b == 0 ? (x = 1, y = 0, a) : (d = extgcd(b, a % b, y, x), y -= a / b * x, d);
}
i64 modinv(i64 a, i64 mod) {
    i64 x, y;
    extgcd(a, mod, x, y);
    x %= mod;
    return x < 0 ? x + mod : x;
}
i64 modpow(i64 a, i64 b, i64 mod) {
    i64 r = 1;
    a %= mod;
    while(b) {
        if(b & 1) r = r * a % mod;
        a = a * a % mod;
        b >>= 1;
    }
    return r;
}
template<i64 mod, i64 primitive>
struct Core {
    Core() {}
    using uint = uint_fast32_t;
    i32 MAX_H;
    vi64 zetaList, zetaInvList;
    Core(i32 MAX_H) : MAX_H(MAX_H), zetaList(MAX_H), zetaInvList(MAX_H) {
        assert((mod & ((1 << MAX_H) - 1)) == 1);

        zetaList.back() = modpow(primitive, (mod - 1) / (1 << MAX_H), mod);
        zetaInvList.back() = modinv(zetaList.back(), mod);
        for(i32 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(vi64 &a, uint n, uint h, bool inv) const {
        static vi64 tmp;
        tmp.resize(n);
        uint mask = n - 1;
        for(uint i = n >> 1, ih = h - 1; i >= 1; i >>= 1, --ih) {
            i64 zeta = inv ? zetaInvList[h - 1 - ih] : zetaList[h - 1 - ih];
            i64 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;
            }
            swap(a, tmp);
        }
        if(inv) {
            i64 invN = modinv(n, mod);
            for(uint i = 0; i < n; ++i) a[i] = a[i] * invN % mod;
        }
    }
    vi64 conv(const vi64 &a, const vi64 &b) const {
        assert(a.size() + b.size() - 1 <= (1 << MAX_H));
 
        if(a.size() == 0) return {};
        if(b.size() == 0) return {};
        return _conv(a, b);
    }
    vi64 _conv(vi64 a, vi64 b) const {
        uint deg = a.size() + b.size() - 1;
        uint h = 0, n = 1;
        while(n < deg) n <<= 1, ++h;
        a.resize(n), b.resize(n);
        return _convStrict(a, b, n, h);
    }
    vi64 _convStrict(vi64 a, vi64 b, uint n, uint h) const {
        fft(a, n, h, 0), fft(b, n, h, 0);
        for(uint i = 0; i < n; ++i) a[i] = a[i] * b[i] % mod;
        fft(a, n, h, 1);
        return a;
    }
};
// NTT_PRIMES {{{
constexpr i64 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
};
// }}}
template <i32 I>
void conv_for(i32 MAX_H, i32 n, i32 h, const vector<i64> &a,
              const vector<i64> &b, vector<i64> &mods,
              vector<i64> &coeffs,
              vector<vector<i64>> &constants) {
    static const Core<NTT_PRIMES[I][0], NTT_PRIMES[I][1]> ntt(MAX_H);
    auto c = ntt._convStrict(a, b, n, h);
    mods[I] = NTT_PRIMES[I][0];
    if(I >= 1) {
        conv_for<I - 1>(MAX_H, n, h, a, b, mods, coeffs, constants);
    }
    // garner
    for(size_t i = 0; i < c.size(); ++i) {
        i64 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 <>
void conv_for< -1 >(i32, i32, i32, const vector<i64> &, const vector<i64> &,
                    vector<i64> &, vector<i64> &,
                    vector<vector<i64>> &) {}
 
template <i32 USE>
vector<i64> conv(i32 MAX_H, vector<i64> a, vector<i64> b,
                  i64 mod) {
    static_assert(USE >= 1, "USE must be positive");
    static_assert(USE <= sizeof(NTT_PRIMES) / sizeof(*NTT_PRIMES), "USE is too big");
    i32 deg = a.size() + b.size() - 1;
    i32 n = 1, h = 0;
    while(n < deg) n <<= 1, ++h;
    a.resize(n), b.resize(n);
    vector<i64> coeffs(USE + 1, 1);
    vector<vector<i64>> constants(USE + 1, vector<i64>(n, 0));
    vector<i64> mods(USE + 1, mod);
    conv_for<USE - 1>(MAX_H, n, h, a, b, mods, coeffs, constants);
    return constants.back();
}
} // namespace NTT
const i32 MAX_H = 18;
NTT::Core<NTT::NTT_PRIMES[0][0], NTT::NTT_PRIMES[0][1]> ntt(MAX_H);


void Main();
i32 main() {cin.tie(nullptr);ios_base::sync_with_stdio(false);cout << fixed << setprecision(15);Main();return 0;}

void Main() {
    //-------------------------------------------------------------------------
    i32 n;
    cin >> n;
    vi64 a(n), b(n);
    for(i32 i = 0; i < n; ++i) cin >> a[i] >> b[i];
    auto c = ntt.conv(a, b);
    cout << 0 << '\n';
    for(i32 i = 0; i < 2 * n - 1; ++i) cout << c[i] << '\n';
    //-------------------------------------------------------------------------
}

Submission Info

Judge Result