一些常用的或不常用的代码模板。

base.cpp

#pragma region base

#include "bits/stdc++.h"

using namespace std;

using i64 = int64_t;
using u64 = uint64_t;
using i32 = int32_t;
using u32 = uint32_t;
using i16 = int16_t;
using u16 = uint16_t;
using i8 = int8_t;
using u8 = uint8_t;

#define $ auto
#define $$ auto&&
#define $C auto const&
#define C const

$ constexpr inf = 0x3f3f3f3f;

#ifndef ONLINE_JUDGE
// 给单调的黑框框加点色彩
$ constexpr input_style = "\033[34;1m"; // 蓝色加亮
$ constexpr print_style = "\033[m"; // 默认白色字
$ constexpr debug_style = "\033[32;1m"; // 绿色加亮
#define DBG(x) cout << debug_style << #x << ": " << (x) << "\n" << input_style
#else
#define DBG(x)
#endif

$ ccxxxi()
{
#ifdef ONLINE_JUDGE
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
#else
    cout << input_style;
#endif
}

// cout << pair
template <typename K, typename V>
auto& operator<<(ostream& o, C pair<K, V>& p)
{
    return o << p.first << " " << p.second << " ";
}

// cout << vector
template <typename E>
auto& operator<<(ostream& o, C vector<E>& v)
{
    for ($$ i : v)
    {
        o << i << " ";
    }
    return o;
}

// cout << array
template <typename T, u32 N>
auto& operator<<(ostream& o, C array<T, N>& arr)
{
    for ($$ i : arr)
    {
        o << i << " ";
    }
    return o;
}

// py风格print单参数特化
template <typename T>
$ print(C T& a)
{
#ifndef ONLINE_JUDGE
    cout << print_style;
#endif
    cout << a << "\n";
#ifndef ONLINE_JUDGE
    cout << input_style;
#endif
}

// py风格print，多参数，空格分隔，尾随换行
template <typename T, typename... Ts>
$ print(C T& a, C Ts&... args)
{
#ifndef ONLINE_JUDGE
    cout << print_style;
#endif
    cout << a << " ";
    print(args...);
}

#pragma endregion

int main()
{
    ccxxxi();

    DBG(1 + 1);

    $C p = pair{"ans", 42};
    $C v = vector{2, 3, 1};
    $C hello = "Hello";
    $C world = "world!";

    print(p, v);
    print(hello, world);
}

binary_indexed_tree.cpp

#pragma region binary indexed tree

$ constexpr low_bit(C u32& i)
{
    return i & ~i + 1;
}

// 树状数组，点更新、范围查询
template <typename T = int>
class bit
{
    vector<T> tree_; // tree是原始数组a的树状数组

public:
    u32 n; // 有效下标1...n

    explicit bit(C u32& n_in) : n(n_in)
    {
        tree_.resize(n + 2);
    }

    // a[i] += d
    $ add(u32 i, C T& d)
    {
        for (; i <= n; i += low_bit(i))
        {
            tree_[i] += d;
        }
    }

    // a[1] + ... + a[i]
    $ prefix_sum(u32 i)C
    {
        T res = 0;
        for (; i >= 1; i -= low_bit(i))
        {
            res += tree_[i];
        }
        return res;
    }

    // a[i] + ... + a[j]
    $ sum(C u32& i, C u32& j)C
    {
        return prefix_sum(j) - prefix_sum(i - 1);
    }

    // a[i]
    $ operator[](C u32& i)C
    {
        return sum(i, i);
    }
};

// 树状数组，范围更新、点查询
template <typename T = int>
class bit_kai
{
    bit<T> tree_; // tree是a的差分数组b的树状数组

public:
    explicit bit_kai(C u32& n) : tree_(n)
    {
    }

    // a[i...j] += d <=> b[i] += d, b[j + 1] -= d
    $ add(C u32& i, C u32& j, C T& d)
    {
        tree_.add(i, d);
        tree_.add(j + 1, -d);
    }

    // a[i] += d
    $ add(C u32& i, C T& d)
    {
        add(i, i, d);
    }

    // a[i] == b[1] + ... + b[i]
    $ operator[](C u32& i)C
    {
        return tree_.prefix_sum(i);
    }
};

// 树状数组，范围更新、范围查询
template <typename T = int>
class bit_ultimate
{
    vector<T> tree1_; // tree1是b的树状数组
    vector<T> tree2_; // tree2是b*i的树状数组
    u32 n_; // 有效下标1...n

    // a[i...n] += d <=> b[i] += d, bi[i] += d * i
    $ suffix_add(u32 i, C T& d)
    {
        $ di = d * static_cast<T>(i);
        for (; i <= n_; i += low_bit(i))
        {
            tree1_[i] += d;
            tree2_[i] += di;
        }
    }

public:
    explicit bit_ultimate(C u32& n_in) : n_(n_in)
    {
        tree1_.resize(n_ + 2);
        tree2_.resize(n_ + 2);
    }

    // a[i...j] += d <=> a[i...n] += d, a[(j+1)...n] -= d
    $ add(C u32& i, C u32& j, C T& d)
    {
        suffix_add(i, d);
        suffix_add(j + 1, -d);
    }

    // a[i] += d
    $ add(C u32& i, C T& d)
    {
        add(i, i, d);
    }

    // a[1] + ... + a[i]
    $ prefix_sum(u32 i)C
    {
        T res = 0;
        $ x = i + 1;
        for (; i >= 1; i -= low_bit(i))
        {
            res += x * tree1_[i] - tree2_[i];
        }
        return res;
    }

    // a[i] + ... + a[j]
    $ sum(C u32& i, C u32& j)C
    {
        return prefix_sum(j) - prefix_sum(i - 1);
    }

    // a[i]
    $ operator[](C u32& i)C
    {
        return sum(i, i);
    }
};

#pragma endregion

graph.cpp

#pragma region graph

// 朴素的并查集
class dsu
{
protected:
    vector<u32> pa_, size_;

    // 路径压缩
    u32 find(C u32& x)
    {
        return pa_[x] == x ? x : pa_[x] = find(pa_[x]);
    }

    // 仅供继承用
    dsu() = default;

public:
    explicit dsu(C u32& sz)
    {
        pa_.resize(sz);
        iota(pa_.begin(), pa_.end(), 0);
        size_.resize(sz, 1);
    }

    // 启发式合并
    $ unite(u32 x, u32 y)
    {
        x = find(x);
        y = find(y);
        if (x != y)
        {
            if (size_[x] > size_[y])
            {
                swap(x, y);
            }
            pa_[x] = y;
            size_[y] += size_[x];
        }
    }

    $ same(C u32& x, C u32& y)
    {
        return find(x) == find(y);
    }

    $ size_of(C u32& x)
    {
        return size_[find(x)];
    }
};

// 带删除操作的并查集
class dsu_with_erase : public dsu
{
public:
    explicit dsu_with_erase(C u32& sz)
    {
        pa_.resize(sz + sz);
        for (u32 i = 0; i < sz; ++i)
        {
            pa_[i] = pa_[i + sz] = i + sz;
        }
        size_.resize(sz + sz, 1);
    }

    $ erase(C u32& x)
    {
        pa_[x] = x;
    }
};

// 带边权的图，存储为邻接表
template <typename Cost = int, Cost Inf = inf>
class graph_with_cost
{
    class edge
    {
    public:
        u32 to;
        Cost cost;
    };

public:
    vector<vector<edge>> g;
    vector<vector<Cost>> mutable dis;
    vector<vector<u32>> mutable pre;
    u32 sz;

    explicit graph_with_cost(C u32& sz_in) : sz(sz_in)
    {
        g.resize(sz + 1);
        dis.resize(sz + 1);
        pre.resize(sz);

        g[sz].reserve(sz);
        for (u32 i = 0; i != sz; ++i)
        {
            g[sz].push_back(edge{i, 0});
        }
    }

    // 添加u到v，权值为cost的单向边
    $ add(C u32& u, C u32& v, C Cost& cost)
    {
        g[u].push_back(edge{v, cost});
    }

    // 添加uv之间，权值为cost的双向边
    $ add2(C u32& u, C u32& v, C Cost& cost)
    {
        add(u, v, cost);
        add(v, u, cost);
    }

    // 检查是否为二分图，即无奇环
    $ binary_check()C
    {
        vector<i8> color(sz);
        function<bool(u32, i8)> dfs = [&](C u32& v, C i8& c)
        {
            color[v] = c;
            for ($C e : g[v])
            {
                $C u = e.to;
                if (color[u] == c)
                {
                    return false;
                }
                if (color[u] == 0 and !dfs(u, -c))
                {
                    return false;
                }
            }
            return true;
        };
        for (u32 i = 0; i != sz; ++i)
        {
            if (color[i] == 0)
            {
                if (!dfs(i, 1))
                {
                    return false;
                }
            }
        }
        return true;
    }

private:
    // 以u到v的边e更新dis[start][v]
    template <bool Path = false>
    $ relax(C u32& start, C u32& u, C edge& e) C
    {
        if (dis[start][u] == Inf)
        {
            return false;
        }
        $C tmp = dis[start][u] + e.cost;
        if (tmp < dis[start][e.to])
        {
            dis[start][e.to] = tmp;
            if constexpr (Path)
            {
                pre[start][e.to] = u;
            }
            return true;
        }
        return false;
    }

public:
    // 求有负权时的单源最短路，基于Bellman-Ford算法，队列优化，O(nm)
    template <bool Path = false>
    $ dis_bf(C u32& start) C
    {
        dis[start].resize(sz, Inf);
        dis[start][start] = 0;
        $ que = queue<u32>{};
        que.push(start);
        $ in_que = vector<bool>(sz);
        if constexpr (Path)
        {
            pre[start].resize(sz, Inf);
        }
        while (!que.empty())
        {
            $ C u = que.front();
            que.pop();
            in_que[u] = false;
            for ($C e : g[u])
            {
                if (relax<Path>(start, u, e) and !in_que[e.to])
                {
                    que.push(e.to);
                    in_que[e.to] = true;
                }
            }
        }
    }

    // 检查是否存在负环，基于Bellman-Ford算法
    $ nl_check()C
    {
        dis[sz].resize(sz + 1, Inf);
        dis[sz][sz] = 0;
        $ que = queue<u32>{};
        que.push(sz);
        $ in_que = vector<bool>(sz + 1);
        $ upd_tm = vector<u32>(sz, 0);
        while (!que.empty())
        {
            $ C u = que.front();
            que.pop();
            in_que[u] = false;
            for ($C e : g[u])
            {
                if (relax(sz, u, e) and !in_que[e.to])
                {
                    if (++upd_tm[e.to] == sz)
                    {
                        return true;
                    }
                    que.push(e.to);
                    in_que[e.to] = true;
                }
            }
        }
        return false;
    }

    // 求无负权时的单源最短路，基于Dijkstra算法，O(mlgm)
    template <bool Path = false>
    $ dis_dij(C u32& start) C
    {
        dis[start].resize(sz, Inf);
        dis[start][start] = 0;
        using p_t = pair<Cost, u32>;
        $ que = priority_queue<p_t, vector<p_t>, greater<>>{};
        que.emplace(Cost(), start);
        if constexpr (Path)
        {
            pre[start].resize(sz, Inf);
        }
        while (!que.empty())
        {
            $ C p = que.top();
            que.pop();
            $ C u = p.second;
            if (dis[start][u] < p.first)
            {
                continue;
            }
            for ($$ e : g[u])
            {
                if (relax<Path>(start, u, e))
                {
                    que.emplace(dis[start][e.to], e.to);
                }
            }
        }
    }

    // 求全源最短路，基于Floyd–Warshall算法，O(n^3)
    template <bool Path = false>
    $ dis_fw()C
    {
        for (u32 i = 0; i != sz; ++i)
        {
            dis[i].resize(sz, Inf);
            if constexpr (Path)
            {
                pre[i].resize(sz, Inf);
            }
            for ($C e : g[i])
            {
                dis[i][e.to] = e.cost;
                if constexpr (Path)
                {
                    pre[i][e.to] = i;
                }
            }
            dis[i][i] = 0;
        }
        for (u32 k = 0; k != sz; ++k)
        {
            for (u32 i = 0; i != sz; ++i)
            {
                for (u32 j = 0; j != sz; ++j)
                {
                    $C tmp = dis[i][k] + dis[k][j];
                    if (tmp < dis[i][j])
                    {
                        dis[i][j] = tmp;
                        if constexpr (Path)
                        {
                            pre[i][j] = k;
                        }
                    }
                }
            }
        }
    }

    // 路径还原，需先求dis[u][v]且指定Path=true
    $ get_path(C u32& u, u32 v) C
    {
        $ ret = vector<u32>{};
        for (; v != inf; v = pre[u][v])
        {
            ret.push_back(v);
        }
        reverse(ret.begin(), ret.end());
        return ret;
    }
};

// 带边权的图，直接存边，用于求最小生成树
template <typename Cost = int>
class graph_with_cost_mst
{
    class edge
    {
    public:
        u32 u, v;
        Cost cost;

        $ operator<(C edge& e) C
        {
            return this->cost < e.cost;
        }
    };

    vector<edge> es_;
    u32 sz_;

public:
    // 以顶点数初始化
    explicit graph_with_cost_mst(C u32& sz_in) : sz_(sz_in)
    {
    }

    // 添加uv之间，权值为cost的双向边
    $ add(C u32& u, C u32& v, C Cost& cost)
    {
        es_.push_back(edge{u, v, cost});
    }

    // 最小生成树，基于Kruskal算法
    $ mst()
    {
        sort(es_.begin(), es_.end());
        dsu dsu(sz_);
        Cost ret = 0;
        for ($C e : es_)
        {
            if (!dsu.same(e.u, e.v))
            {
                dsu.unite(e.u, e.v);
                ret += e.cost;
            }
        }
        return ret;
    }
};

#pragma endregion

k-dimensional_tree.cpp

#pragma region k-dimensional tree

// 辅助函数，相当于先push再pop
template <typename E>
$ push_pop(vector<E>& v, C E& e)
{
    if (e < v.front())
    {
        pop_heap(v.begin(), v.end());
        v.back() = e;
        push_heap(v.begin(), v.end());
    }
}

// k维树，维数为K，坐标type为Crd（不可为unsigned），其他信息type为Other
template <u32 K, typename Crd, typename Other>
class kdt
{
public:
    using crd_arr_t = array<Crd, K>;

    class point
    {
    public:
        crd_arr_t crd;
        Other other;

        point(C crd_arr_t& crd_in, C Other& other_in) : crd(crd_in), other(move(other_in))
        {
        }
    };

private:
    u32 root_ = 0;
    vector<u32> axis_, lc_, rc_;
    vector<point>& points_;

    // 以 [first, last) 中的点建树，返回此树的root
    u32 build(C u32& first, C u32& last)
    {
        $C r = choose_axis(first, last);
        $C num = last - first;
        $C mid = first + num / 2;
        if (num == 1)
        {
        }
        else if (num == 2)
        {
            axis_[mid] = r;
            (points_[first].crd[r] <= points_[mid].crd[r] ? lc_[mid] : rc_[mid]) = first;
        }
        else
        {
            $C b = points_.begin();
            nth_element(b + first, b + mid, b + last,
                        [&](C point& x, C point& y)
                        {
                            return x.crd[r] < y.crd[r];
                        }
            );
            axis_[mid] = r;
            lc_[mid] = build(first, mid);
            rc_[mid] = build(mid + 1, last);
        }
        return mid;
    }

    // 选择 [first, last) 中方差最大的维度
    $ choose_axis(C u32& first, C u32& last) C
    {
        u32 ret_axis = 0;
        float var_max = 0;
        for (u32 i = 0; i < K; ++i)
        {
            $ var_i = variance(first, last, i);
            if (var_i > var_max)
            {
                var_max = var_i;
                ret_axis = i;
            }
        }
        return ret_axis;
    }

    // 计算 [first, last) 中，维度r的方差
    $ variance(C u32& first, C u32& last, C u32& r) C
    {
        $ sum_x = 0.0f, sum_x2 = 0.0f;
        for ($ i = first; i != last; ++i)
        {
            $C tmp = static_cast<float>(points_[i].crd[r]);
            sum_x += tmp;
            sum_x2 += tmp * tmp;
        }
        return sum_x2 - sum_x * sum_x / static_cast<float>(last - first);
    }

public:
    // 以vector<point>初始化，之后外部不应修改此vector
    explicit kdt(vector<point>& points_in) : points_(points_in)
    {
        $C sz = points_.size();
        axis_.resize(sz);
        lc_.resize(sz, inf);
        rc_.resize(sz, inf);
        root_ = build(0, sz);
    }

private:
    class ret_t
    {
    public:
        double dis;
        Other other;
        $ operator<(C ret_t& a) C
        {
            return tie(this->dis, this->other) < tie(a.dis, a.other);
        }
    };

    ret_t none_{numeric_limits<double>::infinity(), Other()};

    // 返回px的欧氏距离的平方，使用浮点数避免平方后溢出
    $ dis2(C crd_arr_t& p, C u32& x)
    {
        double ret = 0;
        for (u32 i = 0; i != K; ++i)
        {
            $C dis1 = static_cast<double>(p[i]) - static_cast<double>(points_[x].crd[i]);
            ret += dis1 * dis1;
        }
        return sqrt(ret);
    }

public:
    // 返回距离点p最近的k个点，欧氏距离
    $ knn(C crd_arr_t& p, C u32& k)
    {
        vector<ret_t> ret(k, none_);
        function<void(u32)> dfs = [&](C u32& x)
        {
            if (x != inf)
            {
                $C r = axis_[x];
                $C dis_sp = p[r] - points_[x].crd[r];
                $C left = dis_sp <= 0;
                dfs(left ? lc_[x] : rc_[x]);
                $C tmp = ret_t{dis2(p, x), points_[x].other};
                push_pop(ret, tmp);
                if (abs(dis_sp) <= ret.front().dis)
                {
                    dfs(left ? rc_[x] : lc_[x]);
                }
            }
        };
        dfs(root_);
        sort_heap(ret.begin(), ret.end());
        ret.erase(lower_bound(ret.begin(), ret.end(), none_), ret.end());
        return ret;
    }
};

// k维树，维数为2的特化
template <typename Crd, typename Other>
class kdt<2, Crd, Other>
{
public:
    using crd_arr_t = array<Crd, 2>;

    class point
    {
    public:
        crd_arr_t crd;
        Other other;

        point(C crd_arr_t& crd_in, C Other& other_in) : crd(crd_in), other(move(other_in))
        {
        }
    };

private:
    u32 root_ = 0;
    vector<u32> lc_, rc_;
    vector<bool> axis_;
    vector<point>& points_;

    // 以 [first, last) 中的点建树，返回此树的root
    u32 build(C u32& first, C u32& last)
    {
        $C r = choose_axis(first, last);
        $C num = last - first;
        $C mid = first + num / 2;
        if (num == 1)
        {
        }
        else if (num == 2)
        {
            axis_[mid] = r;
            (points_[first].crd[r] <= points_[mid].crd[r] ? lc_[mid] : rc_[mid]) = first;
        }
        else
        {
            $C b = points_.begin();
            nth_element(b + first, b + mid, b + last,
                        [&](C point& x, C point& y)
                        {
                            return x.crd[r] < y.crd[r];
                        }
            );
            axis_[mid] = r;
            lc_[mid] = build(first, mid);
            rc_[mid] = build(mid + 1, last);
        }
        return mid;
    }

    // 选择 [first, last) 中方差最大的维度
    $ choose_axis(C u32& first, C u32& last) C
    {
        return variance(first, last, false) < variance(first, last, true);
    }

    // 计算 [first, last) 中，维度r的方差
    $ variance(C u32& first, C u32& last, C bool& r) C
    {
        $ sum_x = 0.0f, sum_x2 = 0.0f;
        for ($ i = first; i != last; ++i)
        {
            $C tmp = static_cast<float>(points_[i].crd[r]);
            sum_x += tmp;
            sum_x2 += tmp * tmp;
        }
        return sum_x2 - sum_x * sum_x / static_cast<float>(last - first);
    }

public:
    // 以vector<point>初始化，之后外部不应修改此vector
    explicit kdt(vector<point>& points_in) : points_(points_in)
    {
        $C sz = points_.size();
        axis_.resize(sz);
        lc_.resize(sz, inf);
        rc_.resize(sz, inf);
        root_ = build(0, sz);
    }

private:
    class ret_t
    {
    public:
        double dis;
        Other other;
        $ operator<(C ret_t& a) C
        {
            return tie(this->dis, this->other) < tie(a.dis, a.other);
        }
    };

    ret_t none_{numeric_limits<double>::infinity(), Other()};

    // 返回px的欧氏距离的平方，使用浮点数避免平方后溢出
    $ dis2(C crd_arr_t& p, C u32& x)
    {
        $C dis_x = static_cast<double>(p[0]) - static_cast<double>(points_[x].crd[0]);
        $C dis_y = static_cast<double>(p[1]) - static_cast<double>(points_[x].crd[1]);
        return sqrt(dis_x * dis_x + dis_y * dis_y);
    }

public:
    // 返回距离点p最近的k个点，欧氏距离
    $ knn(C crd_arr_t& p, C u32& k)
    {
        vector<ret_t> ret(k, none_);
        function<void(u32)> dfs = [&](C u32& x)
        {
            if (x != inf)
            {
                $C r = axis_[x];
                $C dis_sp = p[r] - points_[x].crd[r];
                $C left = dis_sp <= 0;
                dfs(left ? lc_[x] : rc_[x]);
                $C tmp = ret_t{dis2(p, x), points_[x].other};
                push_pop(ret, tmp);
                if (abs(dis_sp) <= ret.front().dis)
                {
                    dfs(left ? rc_[x] : lc_[x]);
                }
            }
        };
        dfs(root_);
        sort_heap(ret.begin(), ret.end());
        ret.erase(lower_bound(ret.begin(), ret.end(), none_), ret.end());
        return ret;
    }
};

#pragma endregion

math.cpp

#pragma region math

#pragma region qpow and prime and phi

// 快速幂，不取模
template <typename Int1, typename Int2>
$ constexpr qpow(Int1 base, Int2 e)
{
    Int1 ret = 1;
    for (; e; e >>= 1, base *= base)
    {
        if (e & 1)
        {
            ret *= base;
        }
    }
    return ret;
}

// 快速乘，避免64位整数溢出
template <$ Mod, typename Int1, typename Int2>
$ constexpr mul(Int1 a, Int2 b)
{
    $ t = static_cast<i64>(
        static_cast<u64>(a) * b
        - static_cast<u64>(static_cast<long double>(a) * b / Mod) * Mod);
    return static_cast<Int1>(t < 0 ? t + Mod : t);
}

// 快速幂，不依赖Mod素性
template <$ Mod, typename Int1, typename Int2>
$ constexpr qpow_not_prime(Int1 base, Int2 e)
{
    static_assert(Mod > 1);
    Int1 ret = 1;
    for (base %= Mod; e; e >>= 1, base = mul<Mod>(base, base))
    {
        if (e & 1)
        {
            ret = mul<Mod>(ret, base);
        }
    }
    return ret;
}

// 快速幂，Mod为素数，基于费马小定理优化
template <$ Mod, typename Int1, typename Int2>
$ constexpr qpow_prime(Int1 base, Int2 e)
{
    return qpow_not_prime<Mod>(base, e % (Mod - 1));
}

// 素数直接检测
class prime
{
    // 此class仅用作封装，不应实例化
    prime() = delete;

public:
    // 基于费马小定理，进行Miller–Rabin素性测试
    $ static constexpr check(C i64& n)
    {
        $ checker = array{2, 3, 5, 7, 11, 13, 17, 19, 23};

        for ($$ a : checker)
        {
            if (n % a == 0)
            {
                return n == a;
            }
        }
        if (n < checker.back())
        {
            return false;
        }

        // n == 2^r * d + 1
        i64 r = 0;
        $ d = n - 1;
        while ((d & 1) == 0)
        {
            d >>= 1;
            ++r;
        }

        // 为了constexpr，mul不能直接调用，让我再次感受到了static_cast的语法设计有多糟糕
        for ($$ a : checker)
        {
            // x = a^d % n
            i64 x = 1;
            $ base = a % n, e = d;
            for (; e != 0; e >>= 1)
            {
                if ((e & 1) != 0)
                {
                    x = static_cast<i64>(
                        static_cast<u64>(x) * base
                        - static_cast<u64>(static_cast<long double>(x) * base / n) * n);
                    if (x < 0)
                    {
                        x += n;
                    }
                }
                base = static_cast<i64>(
                    static_cast<u64>(base) * base
                    - static_cast<u64>(static_cast<long double>(base) * base / n) * n);
                if (base < 0)
                {
                    base += n;
                }
            }

            if (x == 1)
            {
                continue;
            }
            $ ok = false;
            for ($ i = 0; i < r && !ok; ++i)
            {
                if (x == n - 1)
                {
                    ok = true;
                }
                x = static_cast<i64>(
                    static_cast<u64>(x) * x
                    - static_cast<u64>(static_cast<long double>(x) * x / n) * n);
                if (x < 0)
                {
                    x += n;
                }
            }
            if (!ok)
            {
                return false;
            }
        }
        return true;
    }
};

// 素数表
template <u32 N>
class prime_arr
{
public:
    array<u32, N> pri{};
    u32 cnt = 0;

    // 欧拉筛法，O(n)，编译期计算
    constexpr prime_arr()
    {
        array<bool, N> vis{};
        for (u32 i = 2; i < N; ++i)
        {
            if (!vis[i])
            {
                pri[cnt++] = i;
            }
            for (u32 j = 0; j < cnt; ++j)
            {
                if (static_cast<u64>(i) * pri[j] >= N)
                {
                    break;
                }
                vis[i * pri[j]] = true;
                if (!(i % pri[j]))
                {
                    break;
                }
            }
        }
    }
};

// 欧拉函数表
template <u32 N>
class phi_arr
{
public:
    array<u32, N> phi{};

    constexpr phi_arr()
    {
        phi[1] = 1;
        for (u32 i = 2; i < N; ++i)
        {
            if (!phi[i])
            {
                for ($ j = i; j < N; j += i)
                {
                    if (!phi[j])
                    {
                        phi[j] = j;
                    }
                    phi[j] = phi[j] / i * (i - 1);
                }
            }
        }
    }
};

// 素数表 + 欧拉函数表
template <u32 N>
class prime_phi_arr
{
public:
    array<u32, N> pri{};
    array<u32, N> phi{};
    u32 cnt = 0;

    // 欧拉筛法，O(n)，编译期计算
    constexpr prime_phi_arr()
    {
        array<bool, N> vis{};
        phi[1] = 1;
        for (u32 i = 2; i < N; ++i)
        {
            if (!vis[i])
            {
                pri[cnt++] = i;
                phi[i] = i - 1;
            }
            for (u32 j = 0; j < cnt; ++j)
            {
                if (static_cast<u64>(i) * pri[j] >= N)
                {
                    break;
                }
                $ idx = i * pri[j];
                vis[idx] = true;
                if (i % pri[j])
                {
                    phi[idx] = phi[i] * (pri[j] - 1);
                }
                else
                {
                    phi[idx] = phi[i] * pri[j];
                    break;
                }
            }
        }
    }
};

// 快速幂，自动判断Mod是否为素数
template <$ Mod, typename Int1, typename Int2>
$ constexpr qpow(Int1 base, Int2 e)
{
    if constexpr (prime::check(Mod))
    {
        return qpow_prime<Mod>(base, e);
    }
    else
    {
        return qpow_not_prime<Mod>(base, e);
    }
}

#pragma endregion

#pragma region inv and fac and comb

// 求ax+by=gcd(a,b)的解，返回(x,y,d)
$ constexpr ex_gcd(C i64& a, C i64& b)
{
    if (b == 0)
    {
        return array<i64, 3>{1, 0, a};
    }
    $$[x, y, d] = ex_gcd(b, a % b);
    return array{y, x - a / b * y, d};
}

// 逆元
template <$ Mod>
class inv_helper
{
    // 此class仅用作封装，不应实例化
    inv_helper() = delete;

public:
    // 基于费马小定理
    $ static constexpr inv_prime(C i64& x)
    {
        return qpow<Mod>(x, Mod - 2);
    }

    // 基于扩展欧几里得算法
    $ static constexpr inv_not_prime(C i64& x)
    {
        $ ret = ex_gcd(x, Mod)[0];
        if (ret < 0)
        {
            ret += Mod;
        }
        return ret;
    }
};

// 求x的逆元，自动选择合适的算法
template <$ Mod>
$ constexpr inv(C i64& x)
{
    if constexpr (prime::check(Mod))
    {
        return inv_helper<Mod>::inv_prime(x);
    }
    return inv_helper<Mod>::inv_not_prime(x);
}

// 阶乘及其逆元，Start不为0时推广为前缀积
template <$ Mod, u32 N, i64 Start = 0>
class fac_arr
{
    // 否则fac为0逆元无意义
    static_assert(Start + N <= Mod or Start > Mod);
public:
    array<i64, N> fac{};
    array<i64, N> fac_inv{};

    constexpr fac_arr()
    {
        $C start = Start % Mod;
        fac[0] = start != 0 ? start : 1;
        for (u32 i = 1; i < N; ++i)
        {
            fac[i] = fac[i - 1] * ((start + i) % Mod) % Mod;
        }
        fac_inv[N - 1] = inv<Mod>(fac[N - 1]);
        for ($ i = N - 1; i > 0; --i)
        {
            fac_inv[i - 1] = fac_inv[i] * ((start + i) % Mod) % Mod;
        }
    }
};

// 阶乘及其逆元，Start不为0时推广为前缀积，非constexpr版本
template <$ Mod>
class fac_arr_not_constexpr
{
public:
    vector<i64> fac, fac_inv;

    explicit fac_arr_not_constexpr(C u32& n, i64 start)
    {
        assert(start + n <= Mod or start > Mod);
        fac.resize(n);
        fac_inv.resize(n);
        fac[0] = start != 0 ? (start %= Mod) : 1;
        for (u32 i = 1; i < n; ++i)
        {
            fac[i] = fac[i - 1] * ((start + i) % Mod) % Mod;
        }
        fac_inv[n - 1] = inv<Mod>(fac[n - 1]);
        for ($ i = n - 1; i > 0; --i)
        {
            fac_inv[i - 1] = fac_inv[i] * ((start + i) % Mod) % Mod;
        }
    }
};

// 逆元表，[Start, Start + N)
template <$ Mod, u32 N, i64 Start = 0>
class inv_arr
{
public:
    array<i64, N> arr{};

    // 懒癌写法，能用就行
    constexpr inv_arr()
    {
        constexpr fac_arr<Mod, N, Start> tmp;
        for (u32 i = 1; i < N; ++i)
        {
            arr[i] = tmp.fac[i - 1] * tmp.fac_inv[i] % Mod;
        }
        if constexpr (Start != 0)
        {
            arr[0] = inv<Mod>(Start);
        }
    }
};

// 逆元表，[Start, Start + N)，非constexpr版本
template <$ Mod>
class inv_arr_not_constexpr
{
public:
    vector<i64> arr;

    // 懒癌写法，能用就行
    explicit inv_arr_not_constexpr(u32 n, i64 start)
    {
        arr.resize(n);
        fac_arr_not_constexpr<Mod> tmp(n, start);
        for (u32 i = 1; i < n; ++i)
        {
            arr[i] = tmp.fac[i - 1] * tmp.fac_inv[i] % Mod;
        }
        if (start != 0)
        {
            arr[0] = inv<Mod>(start);
        }
    }
};

// 组合数，应满足 n < N or Mod == N
template <$ Mod, u32 N>
class comb
{
    C fac_arr<Mod, N> fac_{};

    // 直接利用组合数公式，需要 n < N
    $ constexpr comb1(C u32& n, C u32& m) C
    {
        return n < m ? 0 : fac_.fac[n] * fac_.fac_inv[m] % Mod * fac_.fac_inv[n - m] % Mod;
    }

    // 基于Lucas定理，需要 Mod == N
    constexpr i64 comb2(C i64& n, C i64& m) C
    {
        return m != 0 ? comb1(static_cast<u32>(n % Mod), static_cast<u32>(m % Mod)) * comb2(n / Mod, m / Mod) % Mod : 1;
    }

public:
    // 求n取m的组合数，自动选择合适的算法
    $ constexpr operator()(C i64& n, C i64& m) C
    {
        if (n < N)
        {
            return comb1(static_cast<u32>(n), static_cast<u32>(m));
        }
        if constexpr (Mod == N)
        {
            return comb2(n, m);
        }
        else
        {
            throw invalid_argument("表太小不够用");
        }
    }
};

#pragma endregion

// 整数平方根，牛顿法，比 std::sqrt 之后再取整大概会快一点点
$ constexpr sqrt_int(C u32& n)
{
    u32 x = 1;
    $ decreased = false;
    while (true)
    {
        $C nx = (x + n / x) / 2;
        if (x == nx or nx > x and decreased)
        {
            break;
        }
        decreased = nx < x;
        x = nx;
    }
    return x;
}

#pragma endregion

int main()
{
    ccxxxi();

    static_assert(qpow(2, 10) == 1024);
    static_assert(qpow<100>(2, 10) == 24);
    static_assert(qpow<1021>(2, 10) == 3);

    static_assert(prime::check(998'244'353));

    $ constexpr pri = prime_arr<200>{};
    static_assert(pri.cnt == 46);
    static_assert(pri.pri[pri.cnt - 1] == 199);

    static_assert(inv<10>(3) == 7);
    static_assert(inv<11>(3) == 4);

    $ constexpr c = comb<7, 7>();
    static_assert(c(5, 2) == 3);

    static_assert(sqrt_int(42) == 6);
}

segment_tree.cpp

#pragma region segment tree

// 线段树，点修改、区间查询
template <typename T = int>
class sgt
{
    u32 n_ = 1;
    vector<T> tree_;
    T init_val_;
    function<T(T, T)> func_;

    $ build()
    {
        for ($ i = n_ - 1; i != 0;)
        {
            --i;
            tree_[i] = func_(tree_[i * 2 + 1], tree_[i * 2 + 2]);
        }
    }

public:
    sgt(C u32& n, C T& init_val, function<T(T, T)> func)
        : init_val_(init_val), func_(move(func))
    {
        while (n_ < n)
        {
            n_ *= 2;
        }
        tree_.resize(n_ * 2 - 1, init_val_);
    }

    sgt(vector<T>&& a, C T& init_val, C function<T(T, T)>& func)
        : sgt(a.size(), init_val, func)
    {
        move(a.begin(), a.end(), tree_.begin() + n_ - 1);
        build();
    }

    sgt(C vector<T>& a, C T& init_val, C function<T(T, T)>& func)
        : sgt(a.size(), init_val, func)
    {
        copy(a.begin(), a.end(), tree_.begin() + n_ - 1);
        build();
    }

    // a[i] = v
    $ set(u32 i, C T& v)
    {
        i += n_ - 1;
        tree_[i] = v;
        while (i != 0)
        {
            i = (i - 1) / 2;
            tree_[i] = func_(tree_[i * 2 + 1], tree_[i * 2 + 2]);
        }
    }

    // 求[i, j)在func意义下的和
    $ query(C u32& i, C u32& j)C
    {
        T res = init_val_;
        function<void(u32, u32, u32)> dfs = [&](C u32& o, C u32& l, C u32& r)
        {
            if (r <= i or j <= l)
            {
                return;
            }
            if (i <= l and r <= j)
            {
                res = func_(res, tree_[o]);
                return;
            }
            $ C m = (l + r) / 2;
            dfs(o * 2 + 1, l, m);
            dfs(o * 2 + 2, m, r);
        };
        dfs(0, 0, n_);
        return res;
    }
};

// 线段树，区间set&add、sum&min&max
template <u32 N>
class sgt_large
{
    u32 n_;

    using v_t = vector<i64>;
    C i64 no_ = INT64_MAX;
    v_t set_ = v_t(N * 4, no_),
        add_ = v_t(N * 4, 0),
        sum_ = v_t(N * 4, 0),
        min_ = v_t(N * 4, 0),
        max_ = v_t(N * 4, 0);

    // 以子节点更新o
    $ maintain(C u32& o, C u32& l, C u32& r)
    {
        $ C num = r - l;
        if (num == 1) // 叶子节点
        {
            sum_[o] = min_[o] = max_[o] = 0;
        }
        else
        {
            $ C lc = o * 2 + 1, rc = o * 2 + 2;
            sum_[o] = sum_[lc] + sum_[rc];
            min_[o] = min(min_[lc], min_[rc]);
            max_[o] = max(max_[lc], max_[rc]);
        }
        if (set_[o] != no_)
        {
            sum_[o] = set_[o] * num;
            min_[o] = max_[o] = set_[o];
        }
        if (add_[o] != 0)
        {
            sum_[o] += add_[o] * num;
            min_[o] += add_[o];
            max_[o] += add_[o];
        }
    }

    // 下传标记
    $ down(C u32& o)
    {
        $ C lc = o * 2 + 1, rc = o * 2 + 2;
        if (set_[o] != no_)
        {
            set_[lc] = set_[rc] = set_[o];
            add_[lc] = add_[rc] = 0;
            set_[o] = no_;
        }
        if (add_[o] != 0)
        {
            add_[lc] += add_[o];
            add_[rc] += add_[o];
            add_[o] = 0;
        }
    }

    template <bool Set>
    $ set_or_add(C u32& i, C u32& j, C i64& v)
    {
        function<void(u32, u32, u32)> dfs = [&](C u32& o, C u32& l, C u32& r)
        {
            if (r <= i or j <= l)
            {
                // pass
            }
            else if (i <= l and r <= j)
            {
                if constexpr (Set)
                {
                    set_[o] = v;
                    add_[o] = 0;
                }
                else
                {
                    add_[o] += v;
                }
            }
            else
            {
                down(o);
                $ C m = (l + r) / 2;
                dfs(o * 2 + 1, l, m);
                dfs(o * 2 + 2, m, r);
            }
            maintain(o, l, r);
        };
        dfs(0, 0, n_);
    }

public:
    explicit sgt_large(C u32& n) : n_(n)
    {
    }

    explicit sgt_large(C v_t& a) : n_(a.size())
    {
        function<void(u32, u32, u32)> dfs = [&](C u32& o, C u32& l, C u32& r)
        {
            if (r - l == 1)
            {
                add_[o] = a[l];
            }
            else
            {
                $ C m = (l + r) / 2;
                dfs(o * 2 + 1, l, m);
                dfs(o * 2 + 2, m, r);
            }
            maintain(o, l, r);
        };
        dfs(0, 0, n_);
    }

    // a[i...(j-1)] = v
    $ set(C u32& i, C u32& j, C i64& v)
    {
        set_or_add<true>(i, j, v);
    }

    // a[i...(j-1)] += v
    $ add(C u32& i, C u32& j, C i64& v)
    {
        set_or_add<false>(i, j, v);
    }

private:
    struct smm
    {
        i64 sum, min, max;
    };

public:
    // 求[i, j)的sum&min&max
    $ query(C u32& i, C u32& j)C
    {
        $ res = smm{0,INT64_MAX,INT64_MIN};
        function<void(u32, u32, u32, i64)> dfs = [&](C u32& o, C u32& l, C u32& r, C i64& add_val)
        {
            if (r <= i or j <= l)
            {
                return;
            }
            if (set_[o] != no_)
            {
                $ C tmp = set_[o] + add_val + add_[o];
                $ C num = min(r, j) - max(l, i);
                res.sum += tmp * num;
                res.min = min(res.min, tmp);
                res.max = max(res.max, tmp);
            }
            else if (i <= l and r <= j)
            {
                res.sum += sum_[o] + add_val * (r - l);
                res.min = min(res.min, min_[o] + add_val);
                res.max = max(res.max, max_[o] + add_val);
            }
            else
            {
                $ C m = (l + r) / 2;
                dfs(o * 2 + 1, l, m, add_val + add_[o]);
                dfs(o * 2 + 2, m, r, add_val + add_[o]);
            }
        };
        dfs(0, 0, n_, 0);
        return res;
    }
};

#pragma endregion

stdc++.h

#pragma once

#include <algorithm>
#include <any>
#include <array>
#include <atomic>
#include <bit>
#include <bitset>
#include <cassert>
#include <cctype>
#include <cerrno>
#include <cfenv>
#include <cfloat>
#include <charconv>
#include <chrono>
#include <cinttypes>
#include <ciso646>
#include <climits>
#include <clocale>
#include <cmath>
#include <codecvt>
#include <compare>
#include <complex>
#include <concepts>
#include <condition_variable>
#include <csetjmp>
#include <csignal>
#include <cstdarg>
#include <cstddef>
#include <cstdint>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <ctime>
#include <cuchar>
#include <cwchar>
#include <cwctype>
#include <deque>
#include <exception>
#include <execution>
#include <filesystem>
#include <forward_list>
#include <fstream>
#include <functional>
#include <future>
#include <initializer_list>
#include <iomanip>
#include <ios>
#include <iosfwd>
#include <iostream>
#include <istream>
#include <iterator>
#include <limits>
#include <list>
#include <locale>
#include <map>
#include <memory>
#include <memory_resource>
#include <mutex>
#include <new>
#include <numbers>
#include <numeric>
#include <optional>
#include <ostream>
#include <queue>
#include <random>
#include <ranges>
#include <ratio>
#include <regex>
#include <scoped_allocator>
#include <set>
#include <shared_mutex>
#include <span>
#include <sstream>
#include <stack>
#include <stdexcept>
#include <streambuf>
#include <string>
#include <string_view>
#include <system_error>
#include <thread>
#include <tuple>
#include <typeindex>
#include <typeinfo>
#include <type_traits>
#include <unordered_map>
#include <unordered_set>
#include <utility>
#include <valarray>
#include <variant>
#include <vector>
#include <version>

string.cpp

#pragma region string

// Z函数
$ z_func(C string& s)
{
    $C n = s.length();
    $ z = vector<u32>(n);
    for (u32 i = 1, l = 0, r = 0; i != n; ++i)
    {
        if (i <= r)
        {
            z[i] = min(r - i + 1, z[i - l]);
        }
        while (i + z[i] < n and s[z[i]] == s[i + z[i]])
        {
            ++z[i];
        }
        if (i + z[i] - 1 > r)
        {
            l = i;
            r = i + z[i] - 1;
        }
    }
    return z;
}

#pragma endregion