Number Theory

Min-Max Inclusion-Exclusion

The core identities for min-max inclusion-exclusion over finite sets are: $$\max(S) = \sum_{T \subseteq S} (-1)^{|T|+1} \min(T)$$ $$\min(S) = \sum_{T \subseteq S} (-1)^{|T|+1} \max(T)$$ Expanding the summation shows all non-extremal values cancel out. The identity also holds under expectation, making it useful for problems where one extremum is trivial to calculate while the other is not.

Example: Expected OR Completion Time

Each second, a number in $[0, 2^n -1]$ is selected with probability $p_i$ and ORed with the current value. Calculate the expected time to reach $2^n -1$.

Let $E(\max(S))$ denote the expected time for the last unset bit in set $S$ to be set, and $E(\min(S))$ the expected time for the first bit in $S$ to be set. Applying the expectation form of the identity: $$E(\max(S)) = \sum_{T \subseteq S} (-1)^{|T|+1} E(\min(T))$$ $E(\min(T))$ is the reciprocal of the probability of selecting any number with atleast one bit in $T$ set: $$E(\min(T)) = \frac{1}{1 - \sum_{U \cap T = \emptyset} p(U)} = \frac{1}{\sum_{U \subseteq \overline{T}} p(U)}$$ We compute the required sum using Fast Walsh-Hadamard Transform (FWT) for OR convolution:

#include <iostream>
#include <cstdio>
#include <cmath>
using namespace std;

const int MAX = 1 << 21;
double prob[MAX];
int len;

void fwtor(double *f, int size) {
    for (int step = 1; step < size; step <<= 1) {
        for (int block = 0; block < size; block += step * 2) {
            for (int pos = 0; pos < step; pos++) {
                f[block + step + pos] += f[block + pos];
            }
        }
    }
}

int main() {
    scanf("%d", &len);
    int full = 1 << len;
    for (int i = 0; i < full; i++) scanf("%lf", &prob[i]);
    fwtor(prob, full);
    double ans = 0.0;
    for (int mask = 1; mask < full; mask++) {
        int bit_cnt = __builtin_popcount(mask);
        double complement_sum = prob[(full - 1) ^ mask];
        if (fabs(1.0 - complement_sum) < 1e-8) continue;
        double sign = (bit_cnt % 2 == 1) ? 1.0 : -1.0;
        ans += sign / (1.0 - complement_sum);
    }
    if (ans < 1e-7) puts("INF");
    else printf("%.6lf\n", ans);
    return 0;
}

Primitive Roots

Multiplicative Order

The smallest positive integer $n$ such that $a^n \equiv 1 \pmod{m}$ for $a \perp m$ is called the multiplicative order of $a$ modulo $m$, denoted $\delta_m(a)$.

Primitive Root Definition

An integer $g$ is a primitive root modulo $m$ if $g \perp m$ and $\delta_m(g) = \varphi(m)$. Only numbers of the form $2, 4, p^k, 2p^k$ (where $p$ is an odd prime) have primitive roots, and the count of primitive roots modulo $n$ is $\varphi(\varphi(n))$.

To find primitive roots:

1. Enumerate small values to find the minimal primitive root $g$
2. All other primitive roots are of the form $g^k$ where $k \perp \varphi(n)$
3. To verify a candidate $g$, check that $g^{\varphi(n)/p} \not\equiv 1 \pmod{n}$ for all prime divisors $p$ of $\varphi(n)$

#include <iostream>
#include <vector>
#include <algorithm>
using namespace std;

const int MAX = 1e6 + 5;
int prime[MAX], prime_cnt = 0, phi[MAX];
bool is_comp[MAX];

void sieve_phi() {
    is_comp[1] = true;
    phi[1] = 1;
    for (int i = 2; i < MAX; i++) {
        if (!is_comp[i]) {
            prime[++prime_cnt] = i;
            phi[i] = i - 1;
        }
        for (int j = 1; j <= prime_cnt && i * prime[j] < MAX; j++) {
            is_comp[i * prime[j]] = true;
            if (i % prime[j] == 0) {
                phi[i * prime[j]] = phi[i] * prime[j];
                break;
            }
            phi[i * prime[j]] = phi[i] * (prime[j] - 1);
        }
    }
}

bool has_primitive_root(int x) {
    if (x == 2 || x == 4) return true;
    if (x % 2 == 0) x >>= 1;
    if (x % 2 == 0) return false;
    for (int i = 2; i <= prime_cnt; i++) {
        if (x % prime[i] == 0) {
            while (x % prime[i] == 0) x /= prime[i];
            return x == 1;
        }
    }
    return false;
}

long long qpow(long long a, long long b, long long mod) {
    long long res = 1;
    while (b) {
        if (b & 1) res = res * a % mod;
        a = a * a % mod;
        b >>= 1;
    }
    return res;
}

int gcd(int a, int b) {
    return b == 0 ? a : gcd(b, a % b);
}

int main() {
    sieve_phi();
    int T;
    scanf("%d", &T);
    while (T--) {
        int n, d;
        scanf("%d%d", &n, &d);
        if (!has_primitive_root(n)) {
            puts("0\n");
            continue;
        }
        vector<int> prime_div;
        int phi_n = phi[n];
        int tmp = phi_n;
        for (int i = 2; i * i <= tmp; i++) {
            if (tmp % i == 0) {
                prime_div.push_back(i);
                while (tmp % i == 0) tmp /= i;
            }
        }
        if (tmp > 1) prime_div.push_back(tmp);
        int min_g = 1;
        while (true) {
            if (gcd(min_g, n) != 1) {
                min_g++;
                continue;
            }
            bool valid = true;
            for (int p : prime_div) {
                if (qpow(min_g, phi_n / p, n) == 1) {
                    valid = false;
                    break;
                }
            }
            if (valid) break;
            min_g++;
        }
        vector<int> roots;
        long long cur = 1;
        for (int k = 1; roots.size() < phi[phi_n]; k++) {
            cur = cur * min_g % n;
            if (gcd(k, phi_n) == 1) roots.push_back(cur);
        }
        sort(roots.begin(), roots.end());
        printf("%d\n", (int)roots.size());
        for (int i = d - 1; i < roots.size(); i += d) {
            printf("%d ", roots[i]);
        }
        puts("");
    }
    return 0;
}

Quadratic Residue

A number $n$ is a quadratic residue modulo an odd prime $p$ if there exists $x$ such that $x^2 \equiv n \pmod{p}$. The solution is unique up to sign modulo $p$.

Euler's criterion: $n$ is a non-residue iff $n^{(p-1)/2} \equiv -1 \pmod{p}$.

To find the root:

1. Randomly select $a$ such that $a^2 - n$ is a non-residue
2. Extend the field with $i^2 = a^2 - n$
3. The solution is $(a + i)^{(p+1)/2} \bmod p$

#include <iostream>
#include <cstdlib>
#include <ctime>
using namespace std;

long long mod, n, i_sq;

struct FieldExt {
    long long real, imag;
};

FieldExt mul(FieldExt a, FieldExt b) {
    FieldExt res;
    res.real = ((a.real * b.real % mod) + (a.imag * b.imag % mod) * i_sq % mod) % mod;
    res.imag = ((a.real * b.imag % mod) + (a.imag * b.real % mod)) % mod;
    if (res.real < 0) res.real += mod;
    if (res.imag < 0) res.imag += mod;
    return res;
}

long long qpow_int(long long a, long long b) {
    long long res = 1;
    while (b) {
        if (b & 1) res = res * a % mod;
        a = a * a % mod;
        b >>= 1;
    }
    return res;
}

FieldExt qpow_ext(FieldExt a, long long b) {
    FieldExt res = {1, 0};
    while (b) {
        if (b & 1) res = mul(res, a);
        a = mul(a, a);
        b >>= 1;
    }
    return res;
}

int main() {
    srand(time(0));
    int T;
    scanf("%d", &T);
    while (T--) {
        scanf("%lld%lld", &n, &mod);
        n %= mod;
        if (qpow_int(n, (mod - 1) >> 1) == mod - 1) {
            puts("Hola!");
            continue;
        }
        long long a;
        while (true) {
            a = rand() % mod;
            long long tmp = (a * a % mod - n + mod) % mod;
            if (qpow_int(tmp, (mod - 1) >> 1) == mod - 1) {
                i_sq = tmp;
                break;
            }
        }
        FieldExt base = {a, 1};
        long long root1 = qpow_ext(base, (mod + 1) >> 1).real % mod;
        long long root2 = (mod - root1) % mod;
        if (root1 > root2) swap(root1, root2);
        if (root1 == root2) printf("%lld\n", root1);
        else printf("%lld %lld\n", root1, root2);
    }
    return 0;
}

Möbius Inversion

Dirichlet convolution of two arithmetic functions is defined as: $$(f * g)(n) = \sum_{d|n} f(d) g(n/d)$$

Key identities: $$(\varphi * 1)(n) = n$$ $$(\mu * 1)(n) = [n = 1]$$

The inversion formula follows directly: $$F = 1 * f \iff f = \mu * F$$

Du Sieve

Du sieve computes prefix sums of multiplicative functions in sublinear time. Let $S(n) = \sum_{i=1}^n f(i)$ and $h = f * g$. From convolution definition: $$\sum_{i=1}^n h(i) = g(1)S(n) + \sum_{i=2}^n g(i)S(\lfloor n/i \rfloor)$$ Rearranged to: $$S(n) = \frac{\sum_{i=1}^n h(i) - \sum_{i=2}^n g(i)S(\lfloor n/i \rfloor)}{g(1)}$$

Example implementation for $\sum \varphi(i)$:

#include <unordered_map>
using namespace std;
typedef long long ll;

const int LIM = 2e6 + 5;
ll phi_pre[LIM];
unordered_map<ll, ll> phi_cache;

ll get_phi_sum(ll n) {
    if (n < LIM) return phi_pre[n];
    if (phi_cache.count(n)) return phi_cache[n];
    ll res = n * (n + 1) / 2;
    for (ll l = 2, r; l <= n; l = r + 1) {
        r = n / (n / l);
        res -= (r - l + 1) * get_phi_sum(n / l);
    }
    return phi_cache[n] = res;
}

Dirichlet Generating Functions

The Dirichlet generating function (DGF) of sequence $f_i$ is: $$\text{DGF}{f}(s) = \sum_{i=1}^\infty \frac{f(i)}{i^s}$$

Key properties:

• The DGF of the constant 1 function is the Riemann zeta function $\zeta(s) = \prod_{p \in \text{Prime}} \frac{1}{1 - p^{-s}}$
• For multiplicative functions, $\text{DGF}{f}(s) = \prod_{p \in \text{Prime}} \sum_{k=0}^\infty \frac{f(p^k)}{p^{ks}}$
• $\text{DGF}{f}(s) \cdot \text{DGF}{g}(s) = \text{DGF}{f * g}(s)$ where $*$ is Dirichlet convolution

Polynomial Operations

Core polynomial operations implemented with NTT modulo 998244353:

#include <iostream>
#include <vector>
using namespace std;

const int MOD = 998244353, G = 3, Gi = 332748118, MAX = 1 << 21;
int bit_rev[MAX];

int qpow(int a, int b) {
    int res = 1;
    while (b) {
        if (b & 1) res = 1LL * res * a % MOD;
        a = 1LL * a * a % MOD;
        b >>= 1;
    }
    return res;
}

struct Poly {
    int deg;
    vector<int> coeff;
    void ntt(int flag) {
        for (int i = 0; i < deg; i++) {
            if (i < bit_rev[i]) swap(coeff[i], coeff[bit_rev[i]]);
        }
        for (int mid = 1; mid < deg; mid <<= 1) {
            int wn = qpow(flag == 1 ? G : Gi, (MOD - 1) / (mid << 1));
            for (int block = 0; block < deg; block += mid << 1) {
                int w = 1;
                for (int pos = 0; pos < mid; pos++) {
                    int x = coeff[block + pos], y = 1LL * w * coeff[block + mid + pos] % MOD;
                    coeff[block + pos] = (x + y) % MOD;
                    coeff[block + mid + pos] = (x - y + MOD) % MOD;
                    w = 1LL * w * wn % MOD;
                }
            }
        }
        if (flag == -1) {
            int inv_deg = qpow(deg, MOD - 2);
            for (int i = 0; i < deg; i++) coeff[i] = 1LL * coeff[i] * inv_deg % MOD;
        }
    }
};

Poly poly_mul(Poly a, Poly b) {
    Poly res;
    int tot_deg = a.deg + b.deg, L = 0;
    res.deg = 1;
    while (res.deg <= tot_deg) res.deg <<= 1, L++;
    for (int i = 0; i < res.deg; i++) bit_rev[i] = (bit_rev[i >> 1] >> 1) | ((i & 1) << (L - 1));
    a.coeff.resize(res.deg); b.coeff.resize(res.deg);
    a.ntt(1); b.ntt(1);
    res.coeff.resize(res.deg);
    for (int i = 0; i < res.deg; i++) res.coeff[i] = 1LL * a.coeff[i] * b.coeff[i] % MOD;
    res.ntt(-1);
    res.deg = tot_deg;
    res.coeff.resize(res.deg + 1);
    return res;
}

Poly poly_inv(Poly a) {
    Poly res;
    res.deg = 1;
    res.coeff.resize(1);
    res.coeff[0] = qpow(a.coeff[0], MOD - 2);
    while (res.deg < a.deg + 1) {
        Poly tmp = res;
        int new_deg = res.deg << 1;
        tmp = poly_mul(tmp, tmp);
        a.coeff.resize(new_deg);
        tmp = poly_mul(tmp, a);
        res.coeff.resize(new_deg);
        for (int i = 0; i < res.deg; i++) res.coeff[i] = (2LL * res.coeff[i] % MOD - tmp.coeff[i] + MOD) % MOD;
        res.deg = new_deg;
    }
    res.deg = a.deg;
    res.coeff.resize(res.deg + 1);
    return res;
}

Poly poly_der(Poly a) {
    Poly res;
    res.deg = a.deg - 1;
    res.coeff.resize(res.deg + 1);
    for (int i = 1; i <= a.deg; i++) res.coeff[i - 1] = 1LL * i * a.coeff[i] % MOD;
    return res;
}

Poly poly_int(Poly a) {
    Poly res;
    res.deg = a.deg + 1;
    res.coeff.resize(res.deg + 1);
    res.coeff[0] = 0;
    for (int i = 0; i <= a.deg; i++) res.coeff[i + 1] = 1LL * a.coeff[i] * qpow(i + 1, MOD - 2) % MOD;
    return res;
}

Poly poly_ln(Poly a) {
    Poly der = poly_der(a);
    Poly inv_a = poly_inv(a);
    Poly mul_res = poly_mul(der, inv_a);
    Poly res = poly_int(mul_res);
    res.deg = a.deg;
    res.coeff.resize(res.deg + 1);
    return res;
}

int main() {
    int n;
    scanf("%d", &n);
    Poly a;
    a.deg = n - 1;
    a.coeff.resize(n);
    for (int i = 0; i < n; i++) scanf("%d", &a.coeff[i]);
    Poly ln_a = poly_ln(a);
    for (int i = 0; i <= ln_a.deg; i++) printf("%d ", ln_a.coeff[i]);
    return 0;
}

Data Structures

Lazy Tag Segment Tree (Range Add/Mul + Range Query)

#include <iostream>
using namespace std;
typedef long long ll;

const int MAX = 2e5 + 5;
int arr[MAX], MOD;
ll tree[MAX << 2], add_tag[MAX << 2], mul_tag[MAX << 2];

ll mod_mul(ll a, ll b) {
    return (a * b) % MOD;
}

ll mod_add(ll a, ll b) {
    return (a + b) % MOD;
}

void push_up(int now) {
    tree[now] = mod_add(tree[now << 1], tree[now << 1 | 1]);
}

void build(int now, int l, int r) {
    add_tag[now] = 0;
    mul_tag[now] = 1;
    if (l == r) {
        tree[now] = arr[l] % MOD;
        return;
    }
    int mid = (l + r) >> 1;
    build(now << 1, l, mid);
    build(now << 1 | 1, mid + 1, r);
    push_up(now);
}

void push_down(int now, int l, int r) {
    if (add_tag[now] == 0 && mul_tag[now] == 1) return;
    int mid = (l + r) >> 1;
    int left = now << 1, right = now << 1 | 1;
    mul_tag[left] = mod_mul(mul_tag[left], mul_tag[now]);
    add_tag[left] = mod_add(mod_mul(add_tag[left], mul_tag[now]), add_tag[now]);
    tree[left] = mod_add(mod_mul(tree[left], mul_tag[now]), mod_mul(add_tag[now], mid - l + 1));
    
    mul_tag[right] = mod_mul(mul_tag[right], mul_tag[now]);
    add_tag[right] = mod_add(mod_mul(add_tag[right], mul_tag[now]), add_tag[now]);
    tree[right] = mod_add(mod_mul(tree[right], mul_tag[now]), mod_mul(add_tag[now], r - mid));
    
    add_tag[now] = 0;
    mul_tag[now] = 1;
}

void range_mul(int now, int l, int r, int ql, int qr, ll val) {
    if (qr < l || ql > r) return;
    if (ql <= l && r <= qr) {
        mul_tag[now] = mod_mul(mul_tag[now], val);
        add_tag[now] = mod_mul(add_tag[now], val);
        tree[now] = mod_mul(tree[now], val);
        return;
    }
    push_down(now, l, r);
    int mid = (l + r) >> 1;
    range_mul(now << 1, l, mid, ql, qr, val);
    range_mul(now << 1 | 1, mid + 1, r, ql, qr, val);
    push_up(now);
}

void range_add(int now, int l, int r, int ql, int qr, ll val) {
    if (qr < l || ql > r) return;
    if (ql <= l && r <= qr) {
        add_tag[now] = mod_add(add_tag[now], val);
        tree[now] = mod_add(tree[now], mod_mul(val, r - l + 1));
        return;
    }
    push_down(now, l, r);
    int mid = (l + r) >> 1;
    range_add(now << 1, l, mid, ql, qr, val);
    range_add(now << 1 | 1, mid + 1, r, ql, qr, val);
    push_up(now);
}

ll range_query(int now, int l, int r, int ql, int qr) {
    if (qr < l || ql > r) return 0;
    if (ql <= l && r <= qr) return tree[now];
    push_down(now, l, r);
    int mid = (l + r) >> 1;
    return mod_add(range_query(now << 1, l, mid, ql, qr), range_query(now << 1 | 1, mid + 1, r, ql, qr));
}

int main() {
    int n, m;
    scanf("%d%d%d", &n, &m, &MOD);
    for (int i = 1; i <= n; i++) scanf("%d", &arr[i]);
    build(1, 1, n);
    while (m--) {
        int opt, l, r, val;
        scanf("%d%d%d", &opt, &l, &r);
        if (opt == 1) {
            scanf("%d", &val);
            range_mul(1, 1, n, l, r, val);
        } else if (opt == 2) {
            scanf("%d", &val);
            range_add(1, 1, n, l, r, val);
        } else if (opt == 3) {
            printf("%lld\n", range_query(1, 1, n, l, r) % MOD);
        }
    }
    return 0;
}

Segment Tree Merge (Tree Path Counting)

#include <iostream>
#include <vector>
using namespace std;

const int MAX = 1e5 + 5;
vector<int> adj[MAX];
int fa[MAX][20], dep[MAX];

void dfs_lca(int now, int parent) {
    fa[now][0] = parent;
    dep[now] = dep[parent] + 1;
    for (int i = 1; i < 20; i++) fa[now][i] = fa[fa[now][i-1]][i-1];
    for (int v : adj[now]) {
        if (v != parent) dfs_lca(v, now);
    }
}

int lca(int u, int v) {
    if (dep[u] < dep[v]) swap(u, v);
    int diff = dep[u] - dep[v];
    for (int i = 0; i < 20; i++) if (diff & (1 << i)) u = fa[u][i];
    if (u == v) return u;
    for (int i = 19; i >= 0; i--) if (fa[u][i] != fa[v][i]) u = fa[u][i], v = fa[v][i];
    return fa[u][0];
}

struct SegNode {
    int max_cnt, max_val;
    SegNode *left = nullptr, *right = nullptr;
    SegNode() : max_cnt(0), max_val(0) {}
} *root[MAX];

SegNode* push_up(SegNode* now) {
    if (!now->left && !now->right) return now;
    if (!now->left) {
        now->max_cnt = now->right->max_cnt;
        now->max_val = now->right->max_val;
        return now;
    }
    if (!now->right) {
        now->max_cnt = now->left->max_cnt;
        now->max_val = now->left->max_val;
        return now;
    }
    if (now->left->max_cnt >= now->right->max_cnt) {
        now->max_cnt = now->left->max_cnt;
        now->max_val = now->left->max_val;
    } else {
        now->max_cnt = now->right->max_cnt;
        now->max_val = now->right->max_val;
    }
    return now;
}

SegNode* insert(SegNode* now, int l, int r, int pos, int delta) {
    if (!now) now = new SegNode();
    if (l == r) {
        now->max_cnt += delta;
        now->max_val = pos;
        return now;
    }
    int mid = (l + r) >> 1;
    if (pos <= mid) now->left = insert(now->left, l, mid, pos, delta);
    else now->right = insert(now->right, mid + 1, r, pos, delta);
    return push_up(now);
}

SegNode* merge(SegNode* a, SegNode* b, int l, int r) {
    if (!a) return b;
    if (!b) return a;
    SegNode* res = new SegNode();
    if (l == r) {
        res->max_cnt = a->max_cnt + b->max_cnt;
        res->max_val = l;
        return res;
    }
    int mid = (l + r) >> 1;
    res->left = merge(a->left, b->left, l, mid);
    res->right = merge(a->right, b->right, mid + 1, r);
    return push_up(res);
}

int ans[MAX];
void dfs_merge(int now, int parent) {
    for (int v : adj[now]) {
        if (v == parent) continue;
        dfs_merge(v, now);
        root[now] = merge(root[now], root[v], 1, MAX - 5);
    }
    ans[now] = root[now] ? (root[now]->max_cnt > 0 ? root[now]->max_val : 0) : 0;
}

int main() {
    int n, m;
    scanf("%d%d", &n, &m);
    for (int i = 1; i < n; i++) {
        int u, v;
        scanf("%d%d", &u, &v);
        adj[u].push_back(v);
        adj[v].push_back(u);
    }
    dfs_lca(1, 0);
    while (m--) {
        int u, v, z;
        scanf("%d%d%d", &u, &v, &z);
        int l = lca(u, v);
        root[u] = insert(root[u], 1, MAX - 5, z, 1);
        root[v] = insert(root[v], 1, MAX - 5, z, 1);
        root[l] = insert(root[l], 1, MAX - 5, z, -1);
        if (fa[l][0]) root[fa[l][0]] = insert(root[fa[l][0]], 1, MAX - 5, z, -1);
    }
    dfs_merge(1, 0);
    for (int i = 1; i <= n; i++) printf("%d\n", ans[i]);
    return 0;
}

Persistent Segment Tree (Static Range K-th Smallest)

#include <iostream>
#include <algorithm>
using namespace std;

const int MAX = 3e5 + 5;
int arr[MAX], sorted[MAX];

struct SegNode {
    int cnt;
    SegNode *left, *right;
    SegNode() : cnt(0), left(nullptr), right(nullptr) {}
} mem[MAX * 20];
int mem_ptr = 0;
SegNode* root[MAX];

SegNode* new_node(SegNode* old) {
    SegNode* res = &mem[mem_ptr++];
    *res = *old;
    return res;
}

void build(SegNode* now, int l, int r) {
    now->cnt = 0;
    if (l == r) return;
    int mid = (l + r) >> 1;
    now->left = &mem[mem_ptr++];
    now->right = &mem[mem_ptr++];
    build(now->left, l, mid);
    build(now->right, mid + 1, r);
}

SegNode* update(SegNode* old, int l, int r, int pos) {
    SegNode* now = new_node(old);
    now->cnt++;
    if (l == r) return now;
    int mid = (l + r) >> 1;
    if (pos <= mid) now->left = update(old->left, l, mid, pos);
    else now->right = update(old->right, mid + 1, r, pos);
    return now;
}

int query(SegNode* a, SegNode* b, int l, int r, int k) {
    if (l == r) return l;
    int mid = (l + r) >> 1;
    int left_cnt = a->left->cnt - b->left->cnt;
    if (left_cnt >= k) return query(a->left, b->left, l, mid, k);
    else return query(a->right, b->right, mid + 1, r, k - left_cnt);
}

int main() {
    int n, m;
    scanf("%d%d", &n, &m);
    for (int i = 1; i <= n; i++) {
        scanf("%d", &arr[i]);
        sorted[i] = arr[i];
    }
    sort(sorted + 1, sorted + n + 1);
    int uniq = unique(sorted + 1, sorted + n + 1) - sorted - 1;
    root[0] = &mem[mem_ptr++];
    build(root[0], 1, uniq);
    for (int i = 1; i <= n; i++) {
        int pos = lower_bound(sorted + 1, sorted + uniq + 1, arr[i]) - sorted;
        root[i] = update(root[i-1], 1, uniq, pos);
    }
    while (m--) {
        int l, r, k;
        scanf("%d%d%d", &l, &r, &k);
        int res_pos = query(root[r], root[l-1], 1, uniq, k);
        printf("%d\n", sorted[res_pos]);
    }
    return 0;
}

Splay Tree (Interval Operations)

Supports interval flip, point update, and maximum subarray sum queries:

#include <iostream>
using namespace std;
typedef long long ll;

const int MAX = 2e5 + 5, INF = 0x3f3f3f3f;
int ch[MAX][2], fa[MAX], siz[MAX], rev_tag[MAX], val[MAX], node_cnt;
ll ls[MAX], rs[MAX], sum[MAX], mx[MAX];

#define ls(x) ch[x][0]
#define rs(x) ch[x][1]
#define is_right(x) (x == ch[fa[x]][1])

ll max2(ll a, ll b) { return a > b ? a : b; }

void push_up(int x) {
    siz[x] = siz[ls(x)] + siz[rs(x)] + 1;
    sum[x] = sum[ls(x)] + sum[rs(x)] + val[x];
    ls[x] = max2(ls[ls(x)], sum[ls(x)] + val[x] + ls[rs(x)]);
    rs[x] = max2(rs[rs(x)], sum[rs(x)] + val[x] + rs[ls(x)]);
    mx[x] = max2(max2(mx[ls(x)], mx[rs(x)]), rs[ls(x)] + val[x] + ls[rs(x)]);
}

void push_down(int x) {
    if (rev_tag[x]) {
        swap(ls[x], rs[x]);
        swap(ls(ls(x)), rs(ls(x)));
        swap(ls(rs(x)), rs(rs(x)));
        swap(ls(x), rs(x));
        rev_tag[ls(x)] ^= 1;
        rev_tag[rs(x)] ^= 1;
        rev_tag[x] = 0;
    }
}

int new_node(ll v) {
    int x = ++node_cnt;
    val[x] = v;
    sum[x] = v;
    ls[x] = rs[x] = mx[x] = max2(v, 0LL);
    siz[x] = 1;
    return x;
}

void rotate(int x) {
    int f = fa[x], ff = fa[f], c = is_right(x);
    if (ff) ch[ff][is_right(f)] = x;
    fa[x] = ff;
    ch[f][c] = ch[x][c ^ 1];
    if (ch[x][c ^ 1]) fa[ch[x][c ^ 1]] = f;
    ch[x][c ^ 1] = f;
    fa[f] = x;
    push_up(f);
    push_up(x);
}

void splay(int x, int top = 0) {
    while (fa[x] != top) {
        int f = fa[x], ff = fa[f];
        if (ff != top) rotate(is_right(x) == is_right(f) ? f : x);
        rotate(x);
    }
    if (!fa[x]) ::ls[0] = ::rs[0] = ::mx[0] = -INF; sum[0] = 0; siz[0] = 0;
}

int kth(int k) {
    int now = 0;
    while (now) {
        push_down(now);
        if (k <= siz[ls(now)]) now = ls(now);
        else if (k == siz[ls(now)] + 1) break;
        else k -= siz[ls(now)] + 1, now = rs(now);
    }
    return now;
}

void reverse_interval(int l, int r) {
    int L = kth(l), R = kth(r + 2);
    splay(L);
    splay(R, L);
    rev_tag[ls(rs(ls(0)))] ^= 1;
}

ll query_max_sub(int l, int r) {
    int L = kth(l), R = kth(r + 2);
    splay(L);
    splay(R, L);
    return mx[ls(rs(ls(0)))];
}

int main() {
    int n, m;
    scanf("%d%d", &n, &m);
    new_node(-INF);
    new_node(-INF);
    for (int i = 1; i <= n; i++) {
        ll v;
        scanf("%lld", &v);
        new_node(v);
    }
    while (m--) {
        char opt;
        int l, r;
        scanf(" %c%d%d", &opt, &l, &r);
        if (opt == 'R') reverse_interval(l, r);
        else if (opt == 'Q') printf("%lld\n", query_max_sub(l, r));
        else if (opt == 'M') {
            int pos = kth(l + 1);
            val[pos] = r;
            push_up(pos);
            splay(pos);
        }
    }
    return 0;
}

CDQ Divide and Conquer (3D Partial Order)

#include <iostream>
#include <algorithm>
using namespace std;

const int MAX = 1e5 + 5;
struct Point {
    int x, y, z, cnt, *ans;
    bool operator == (const Point& other) const {
        return x == other.x && y == other.y && z == other.z;
    }
} a[MAX], tmp1[MAX], tmp2[MAX];

int ans[MAX], res[MAX], n, k;

bool cmp(const Point& a, const Point& b) {
    if (a.x != b.x) return a.x < b.x;
    if (a.y != b.y) return a.y < b.y;
    return a.z < b.z;
}

void cdq2(int l, int r) {
    if (l >= r) return;
    int mid = (l + r) >> 1;
    cdq2(l, mid);
    cdq2(mid + 1, r);
    int i = l, j = mid + 1, ptr = l, add = 0;
    while (ptr <= r) {
        if (i <= mid && (j > r || tmp1[i].z <= tmp1[j].z)) {
            tmp2[ptr] = tmp1[i++];
            add += tmp2[ptr].cnt;
        } else {
            tmp2[ptr] = tmp1[j++];
            if (!tmp2[ptr].cnt) *tmp2[ptr].ans += add;
        }
        ptr++;
    }
    for (int p = l; p <= r; p++) tmp1[p] = tmp2[p];
}

void cdq(int l, int r) {
    if (l >= r) return;
    int mid = (l + r) >> 1;
    cdq(l, mid);
    cdq(mid + 1, r);
    int i = l, j = mid + 1, ptr = l;
    while (ptr <= r) {
        if (i <= mid && (j > r || a[i].y <= a[j].y)) {
            tmp1[ptr] = a[i++];
            tmp1[ptr].cnt = 1;
        } else {
            tmp1[ptr] = a[j++];
            tmp1[ptr].cnt = 0;
        }
        ptr++;
    }
    for (int p = l; p <= r; p++) a[p] = tmp1[p];
    cdq2(l, r);
}

int main() {
    scanf("%d%d", &n, &k);
    for (int i = 1; i <= n; i++) {
        scanf("%d%d%d", &a[i].x, &a[i].y, &a[i].z);
        a[i].ans = &ans[i];
        ans[i] = 0;
    }
    sort(a + 1, a + n + 1, cmp);
    for (int i = n - 1; i >= 1; i--) {
        if (a[i] == a[i + 1]) *a[i].ans = *a[i + 1].ans + 1;
    }
    cdq(1, n);
    for (int i = 1; i <= n; i++) res[ans[i]]++;
    for (int i = 0; i < n; i++) printf("%d\n", res[i]);
    return 0;
}

Graph Theory

Dinic's Max Flow

#include <iostream>
#include <vector>
#include <queue>
#include <cstring>
using namespace std;
typedef long long ll;

const ll INF = 1e18;
struct Edge {
    int to, rev;
    ll cap;
    Edge(int t, int r, ll c) : to(t), rev(r), cap(c) {}
};

template <int MAXN>
class MaxFlow {
private:
    vector<Edge> adj[MAXN];
    int dep[MAXN], ptr[MAXN], s, t;

    bool bfs() {
        memset(dep, -1, sizeof dep);
        queue<int> q;
        dep[s] = 0;
        q.push(s);
        while (!q.empty()) {
            int u = q.front();
            q.pop();
            for (Edge& e : adj[u]) {
                if (e.cap > 0 && dep[e.to] == -1) {
                    dep[e.to] = dep[u] + 1;
                    q.push(e.to);
                    if (e.to == t) return true;
                }
            }
        }
        return false;
    }

    ll dfs(int u, ll flow) {
        if (u == t) return flow;
        for (int& i = ptr[u]; i < adj[u].size(); i++) {
            Edge& e = adj[u][i];
            if (e.cap > 0 && dep[e.to] == dep[u] + 1) {
                ll pushed = dfs(e.to, min(flow, e.cap));
                if (pushed > 0) {
                    e.cap -= pushed;
                    adj[e.to][e.rev].cap += pushed;
                    return pushed;
                }
            }
        }
        return 0;
    }

public:
    void add_edge(int u, int v, ll cap) {
        adj[u].emplace_back(v, adj[v].size(), cap);
        adj[v].emplace_back(u, adj[u].size() - 1, 0);
    }

    ll max_flow(int s_, int t_) {
        s = s_, t = t_;
        ll total = 0;
        while (bfs()) {
            memset(ptr, 0, sizeof ptr);
            while (ll pushed = dfs(s, INF)) total += pushed;
        }
        return total;
    }
};

MaxFlow<405> mf;

int main() {
    int n, m, s, t;
    scanf("%d%d%d%d", &n, &m, &s, &t);
    while (m--) {
        int u, v;
        ll cap;
        scanf("%d%d%lld", &u, &v, &cap);
        mf.add_edge(u, v, cap);
    }
    printf("%lld\n", mf.max_flow(s, t));
    return 0;
}

Kruskal Reconstruction Tree

Used for queries of the form: find all nodes reachable from a starting node using edges with weight <= x.

For two nodes u and v, the weight of their LCA on the reconstruction tree equals the minimum possible maximum edge weight along any path between u and v.

#include <iostream>
#include <vector>
#include <algorithm>
using namespace std;

const int MAX = 2e5 + 5;
struct Edge {
    int u, v, w;
    bool operator < (const Edge& other) const {
        return w < other.w;
    }
} edges[MAX * 3];

int dsu[MAX], node_cnt, n, m, q;
int find(int x) { return dsu[x] == x ? x : dsu[x] = find(dsu[x]); }

vector<int> adj[MAX];
int val[MAX], fa[MAX][20], siz[MAX], in[MAX], dfn_cnt, height[MAX];
int a[MAX], sorted[MAX];

struct PartitionTree {
    int num[20][MAX], val[20][MAX];
    void build(int l, int r, int level) {
        if (l == r) return;
        int mid = (l + r) >> 1, same = mid - l + 1;
        for (int i = l; i <= r; i++) if (val[level][i] < sorted[mid]) same--;
        int lp = l, rp = mid + 1;
        for (int i = l; i <= r; i++) {
            num[level][i] = (i == l) ? 0 : num[level][i-1];
            if (val[level][i] < sorted[mid] || (val[level][i] == sorted[mid] && same > 0)) {
                val[level+1][lp++] = val[level][i];
                num[level][i]++;
                if (val[level][i] == sorted[mid]) same--;
            } else val[level+1][rp++] = val[level][i];
        }
        build(l, mid, level + 1);
        build(mid + 1, r, level + 1);
    }

    int query(int l, int r, int ql, int qr, int k, int level) {
        if (l == r) return val[level][l];
        int mid = (l + r) >> 1;
        int left = (ql == l) ? 0 : num[level][ql - 1];
        int cnt = num[level][qr] - left;
        if (cnt >= k) {
            return query(l, mid, l + left, l + num[level][qr] - 1, k, level + 1);
        } else {
            int nl = mid + 1 + (ql - l - left);
            int nr = nl + (qr - ql + 1 - cnt) - 1;
            return query(mid + 1, r, nl, nr, k - cnt, level + 1);
        }
    }
} pt;

void dfs(int u, int p) {
    fa[u][0] = p;
    height[u] = height[p] + 1;
    for (int i = 1; i < 20; i++) fa[u][i] = fa[fa[u][i-1]][i-1];
    if (u <= n) {
        in[u] = ++dfn_cnt;
        pt.val[0][dfn_cnt] = a[u];
        siz[u] = 1;
        return;
    }
    siz[u] = 0;
    for (int v : adj[u]) {
        if (v == p) continue;
        dfs(v, u);
        siz[u] += siz[v];
        in[u] = min(in[u], in[v]);
    }
}

int jump(int u, int limit) {
    for (int i = 19; i >= 0; i--) {
        if (val[fa[u][i]] <= limit) u = fa[u][i];
    }
    return u;
}

int main() {
    scanf("%d%d%d", &n, &m, &q);
    node_cnt = n;
    for (int i = 1; i <= n; i++) {
        scanf("%d", &a[i]);
        sorted[i] = a[i];
        dsu[i] = i;
    }
    sort(sorted + 1, sorted + n + 1);
    for (int i = 1; i <= m; i++) scanf("%d%d%d", &edges[i].u, &edges[i].v, &edges[i].w);
    sort(edges + 1, edges + m + 1);
    val[0] = 1e9;
    for (int i = 1; i <= m; i++) {
        int u = find(edges[i].u), v = find(edges[i].v);
        if (u == v) continue;
        int new_node = ++node_cnt;
        dsu[new_node] = dsu[u] = dsu[v] = new_node;
        val[new_node] = edges[i].w;
        adj[new_node].push_back(u);
        adj[new_node].push_back(v);
    }
    dfs(node_cnt, 0);
    pt.build(1, n, 0);
    while (q--) {
        int u, x, k;
        scanf("%d%d%d", &u, &x, &k);
        int root = jump(u, x);
        if (siz[root] < k) {
            puts("-1");
            continue;
        }
        int res = pt.query(1, n, in[root], in[root] + siz[root] - 1, siz[root] - k + 1, 0);
        printf("%d\n", res);
    }
    return 0;
}

Strict Second Minimum Spanning Tree

The strict second minimum spanning tree differs from the MST by exactly one edge. We enumerate non-MST edges and replace the maximum edge on the path between its endpoints to find the minimal increase.

#include <iostream>
#include <vector>
#include <algorithm>
using namespace std;
typedef long long ll;

const int MAX = 1e5 + 5, INF = 1e9;
struct Edge {
    int u, v, w;
    bool in_mst;
    bool operator < (const Edge& other) const {
        return w < other.w;
    }
} edges[MAX * 3];

vector<pair<int, int>> adj[MAX];
int dsu[MAX], n, m;
ll mst_sum = 0;

int find(int x) { return dsu[x] == x ? x : dsu[x] = find(dsu[x]); }

void kruskal() {
    sort(edges + 1, edges + m + 1);
    for (int i = 1; i <= n; i++) dsu[i] = i;
    for (int i = 1; i <= m; i++) {
        int u = edges[i].u, v = edges[i].v, w = edges[i].w;
        if (find(u) != find(v)) {
            edges[i].in_mst = true;
            dsu[find(u)] = find(v);
            mst_sum += w;
            adj[u].emplace_back(v, w);
            adj[v].emplace_back(u, w);
        }
    }
}

int fa[MAX][20], max1[MAX][20], max2[MAX][20], dep[MAX];

void dfs(int u, int p, int w) {
    fa[u][0] = p;
    max1[u][0] = w;
    max2[u][0] = -1;
    dep[u] = dep[p] + 1;
    for (int i = 1; i < 20; i++) {
        fa[u][i] = fa[fa[u][i-1]][i-1];
        int prev1 = max1[u][i-1], prev2 = max1[fa[u][i-1]][i-1];
        int prev2_1 = max2[u][i-1], prev2_2 = max2[fa[u][i-1]][i-1];
        max1[u][i] = max(prev1, prev2);
        if (prev1 != prev2) max2[u][i] = max(min(prev1, prev2), max(prev2_1, prev2_2));
        else max2[u][i] = max(prev2_1, prev2_2);
    }
    for (auto& e : adj[u]) {
        int v = e.first, ew = e.second;
        if (v != p) dfs(v, u, ew);
    }
}

pair<int, int> get_max(int u, int v) {
    int m1 = 0, m2 = 0;
    if (dep[u] < dep[v]) swap(u, v);
    int diff = dep[u] - dep[v];
    for (int i = 0; i < 20; i++) {
        if (diff & (1 << i)) {
            int tmp1 = max1[u][i], tmp2 = max2[u][i];
            if (tmp1 > m1) {
                m2 = max(m1, tmp2);
                m1 = tmp1;
            } else if (tmp1 < m1) {
                m2 = max(m2, tmp1);
            } else {
                m2 = max(m2, tmp2);
            }
            u = fa[u][i];
        }
    }
    if (u == v) return {m1, m2};
    for (int i = 19; i >= 0; i--) {
        if (fa[u][i] != fa[v][i]) {
            int tmp1_u = max1[u][i], tmp2_u = max2[u][i];
            int tmp1_v = max1[v][i], tmp2_v = max2[v][i];
            int t1 = max(tmp1_u, tmp1_v);
            int t2 = max(min(tmp1_u, tmp1_v), max(tmp2_u, tmp2_v));
            if (t1 > m1) {
                m2 = max(m1, t2);
                m1 = t1;
            } else if (t1 < m1) {
                m2 = max(m2, t1);
            } else {
                m2 = max(m2, t2);
            }
            u = fa[u][i], v = fa[v][i];
        }
    }
    int tmp1_u = max1[u][0], tmp2_u = max2[u][0];
    int tmp1_v = max1[v][0], tmp2_v = max2[v][0];
    int t1 = max(tmp1_u, tmp1_v);
    int t2 = max(min(tmp1_u, tmp1_v), max(tmp2_u, tmp2_v));
    if (t1 > m1) {
        m2 = max(m1, t2);
        m1 = t1;
    } else if (t1 < m1) {
        m2 = max(m2, t1);
    } else {
        m2 = max(m2, t2);
    }
    return {m1, m2};
}

int main() {
    scanf("%d%d", &n, &m);
    for (int i = 1; i <= m; i++) {
        scanf("%d%d%d", &edges[i].u, &edges[i].v, &edges[i].w);
        edges[i].in_mst = false;
    }
    kruskal();
    dfs(1, 0, 0);
    ll min_diff = 1e18;
    for (int i = 1; i <= m; i++) {
        if (edges[i].in_mst) continue;
        int u = edges[i].u, v = edges[i].v, w = edges[i].w;
        auto [m1, m2] = get_max(u, v);
        if (w > m1) min_diff = min(min_diff, (ll)w - m1);
        else if (w > m2) min_diff = min(min_diff, (ll)w - m2);
    }
    printf("%lld\n", mst_sum + min_diff);
    return 0;
}

Search

IDA* for K-th Shortest Path

#include <iostream>
#include <vector>
#include <queue>
#include <cstring>
using namespace std;

const int MAX = 2e4 + 5;
const double INF = 1e7;
vector<pair<int, double>> adj[MAX], rev_adj[MAX];
double h[MAX];
bool vis[MAX];
int n, m;

dijkstra_heuristic() {
    for (int i = 1; i <= n; i++) h[i] = INF;
    h[n] = 0;
    memset(vis, 0, sizeof vis);
    for (int i = 1; i <= n; i++) {
        int u = 0;
        for (int j = 1; j <= n; j++) {
            if (!vis[j] && h[j] < h[u]) u = j;
        }
        if (u == 0) break;
        vis[u] = true;
        for (auto& e : rev_adj[u]) {
            int v = e.first;
            double w = e.second;
            if (h[v] > h[u] + w) h[v] = h[u] + w;
        }
    }
}

struct Node {
    int u;
    double g;
    Node(int u_, double g_) : u(u_), g(g_) {}
    bool operator > (const Node& other) const {
        return g + h[u] > other.g + h[other.u];
    }
};

priority_queue<Node, vector<Node>, greater<Node>> pq;
int cnt[MAX], ans = 0;
double max_cost;

void a_star(int limit) {
    pq.emplace(1, 0.0);
    while (!pq.empty() && pq.size() < 4e6) {
        auto [u, g] = pq.top();
        pq.pop();
        cnt[u]++;
        if (g > max_cost) return;
        if (u == n) {
            ans++;
            max_cost -= g;
            continue;
        }
        if (cnt[u] > limit) continue;
        for (auto& e : adj[u]) {
            int v = e.first;
            double w = e.second;
            pq.emplace(v, g + w);
        }
    }
}

int main() {
    scanf("%d%d%lf", &n, &m, &max_cost);
    for (int i = 1; i <= m; i++) {
        int u, v;
        double w;
        scanf("%d%d%lf", &u, &v, &w);
        adj[u].emplace_back(v, w);
        rev_adj[v].emplace_back(u, w);
    }
    dijkstra_heuristic();
    if (h[1] > max_cost) {
        puts("0");
        return 0;
    }
    int limit = (int)(max_cost / h[1]) + 1;
    a_star(limit);
    printf("%d\n", ans);
    return 0;
}