Number Theory

Fast Exponentiation

Computes base raised to the power of exp modulo mod efficiently using binary decomposition.

long long fast_power(long long base, long long exp, long long mod) {
    long long result = 1;
    base %= mod;
    while (exp > 0) {
        if (exp & 1) result = (result * base) % mod;
        base = (base * base) % mod;
        exp >>= 1;
    }
    return result;
}

Fast Multiplication

Computes (a * b) % mod using binary methods to avoid 64-bit overflow before the modulo operation (useful when values can exceed 2^63).

long long fast_mul(long long a, long long b, long long mod) {
    long long result = 0;
    while (b > 0) {
        if (b & 1) result = (result + a) % mod;
        a = (a << 1) % mod;
        b >>= 1;
    }
    return result;
}

Linear Sieve for Primes

Generates all prime numbers up to n in O(n) time complexity.

const int MAXN = 1000005;
int is_composite[MAXN] = {0};
std::vector<int> prime_list;

void linear_sieve(int n) {
    for (int i = 2; i <= n; ++i) {
        if (!is_composite[i]) {
            prime_list.push_back(i);
        }
        for (int p : prime_list) {
            if (i * p > n) break;
            is_composite[i * p] = 1;
            if (i % p == 0) break;
        }
    }
}
</int>

Linear Inverse

Computes modular inverses for all numbers from 1 to n modulo MOD, assuming MOD is prime.

const int MOD = 1e9 + 7;
long long inv[MAXN];

void compute_linear_inverses(int n) {
    inv[1] = 1;
    for (int i = 2; i <= n; ++i) {
        inv[i] = (MOD - MOD / i) * inv[MOD % i] % MOD;
    }
}

Extended Eucldiean Algorithm

Solves the equation ax + by = gcd(a, b) and returns the GCD along with coefficeints x and y.

long long extended_gcd(long long a, long long b, long long &x, long long &y) {
    if (b == 0) {
        x = 1;
        y = 0;
        return a;
    }
    long long d = extended_gcd(b, a % b, y, x);
    y -= a / b * x;
    return d;
}

Chinese Remainder Theorem (CRT)

Solves a system of congruences x ≡ a_i (mod m_i) where m_i are paiwrise coprime.

long long crt(const std::vector<long long="">& remainders, const std::vector<long long="">& moduli) {
    long long product = 1;
    for (long long m : moduli) product *= m;

    long long result = 0;
    int n = remainders.size();
    for (int i = 0; i < n; ++i) {
        long long partial_product = product / moduli[i];
        long long inv, _;
        extended_gcd(partial_product, moduli[i], inv, _);
        inv = (inv % moduli[i] + moduli[i]) % moduli[i];
        result = (result + remainders[i] * partial_product % product * inv % product) % product;
    }
    return (result + product) % product;
}
</long></long>

Discrete Logarithm (BSGS)

Solves a^x ≡ b (mod p) using the Baby-Step Giant-Step algorithm.

#include <unordered_map>
#include <cmath>

long long bsgs(long long a, long long b, long long p) {
    std::unordered_map<long long=""> table;
    b %= p;
    long long t = std::sqrt(p) + 1;
    
    for (long long j = 0; j < t; ++j) {
        long long val = b * fast_power(a, j, p) % p;
        table[val] = j;
    }
    
    long long a_t = fast_power(a, t, p);
    if (a_t == 0) return b == 0 ? 1 : -1;
    
    for (long long i = 0; i <= t; ++i) {
        long long val = fast_power(a_t, i, p);
        if (table.count(val)) {
            long long ans = i * t - table[val];
            if (ans >= 0) return ans;
        }
    }
    return -1;
}
</long></cmath></unordered_map>

Gaussian Elimination

Solves a system of linear equations.

const int MAXN = 105;
double matrix[MAXN][MAXN];
double solution[MAXN];

void gaussian_elimination(int n) {
    for (int i = 1; i <= n; ++i) {
        int pivot = i;
        for (int j = i + 1; j <= n; ++j) {
            if (std::abs(matrix[j][i]) > std::abs(matrix[pivot][i])) pivot = j;
        }
        if (std::abs(matrix[pivot][i]) < 1e-9) {
            // No unique solution
            return; 
        }
        for (int j = 1; j <= n + 1; ++j) std::swap(matrix[i][j], matrix[pivot][j]);
        
        for (int j = i + 1; j <= n; ++j) {
            double factor = matrix[j][i] / matrix[i][i];
            for (int k = i; k <= n + 1; ++k) {
                matrix[j][k] -= matrix[i][k] * factor;
            }
        }
    }
    
    for (int i = n; i >= 1; --i) {
        solution[i] = matrix[i][n + 1];
        for (int j = i + 1; j <= n; ++j) {
            solution[i] -= matrix[i][j] * solution[j];
        }
        solution[i] /= matrix[i][i];
    }
}

Graph Theory

Dijkstra's Algorithm

Computes the shortest paths from a source node to all other nodes in a graph with non-negative edge weights.

#include <queue>
#include <vector>
#include <cstring>

const int INF = 0x3f3f3f3f;
struct Edge { int to, weight; };
std::vector<edge> adj[MAXN];
long long dist[MAXN];
bool visited[MAXN];

void dijkstra(int start_node) {
    std::priority_queue<:pair int="" long="">, std::vector<:pair int="" long="">>, std::greater<:pair int="" long="">>> pq;
    
    std::memset(dist, 0x3f, sizeof(dist));
    dist[start_node] = 0;
    pq.push({0, start_node});
    
    while (!pq.empty()) {
        int u = pq.top().second;
        long long d = pq.top().first;
        pq.pop();
        
        if (visited[u]) continue;
        visited[u] = true;
        
        for (const auto& edge : adj[u]) {
            int v = edge.to;
            long long w = edge.weight;
            if (dist[v] > dist[u] + w) {
                dist[v] = dist[u] + w;
                pq.push({dist[v], v});
            }
        }
    }
}
</:pair></:pair></:pair></edge></cstring></vector></queue>

Data Structures & Dynamic Programming

Matrix Class

A generic matrix class supporting arithmetic operations, exponentiation, and inversion for solving linear equations or recurrences.

const long long MOD = 1e9 + 7;

class Matrix {
public:
    int rows, cols;
    std::vector<:vector long="">> data;

    Matrix(int r = 0, int c = 0) : rows(r), cols(c) {
        data.assign(rows, std::vector<long long="">(cols, 0));
        if (rows == cols) {
            for (int i = 0; i < rows; ++i) data[i][i] = 1;
        }
    }

    Matrix operator*(const Matrix& other) const {
        Matrix result(rows, other.cols);
        for (int i = 0; i < rows; ++i) {
            for (int k = 0; k < cols; ++k) {
                for (int j = 0; j < other.cols; ++j) {
                    result.data[i][j] = (result.data[i][j] + data[i][k] * other.data[k][j]) % MOD;
                }
            }
        }
        return result;
    }

    Matrix operator^(long long power) const {
        Matrix result(rows, rows);
        Matrix base = *this;
        while (power > 0) {
            if (power & 1) result = result * base;
            base = base * base;
            power >>= 1;
        }
        return result;
    }
};
</long></:vector>

Heavy-Light Decomposition (HLD)

Decomposes a tree into chains to allow efficient queries on paths (e.g., sum, max) using a segment tree.

int parent[MAXN], heavy_child[MAXN], depth[MAXN];
int head[MAXN], position[MAXN], current_pos;
int subtree_size[MAXN];
int base_array[MAXN], arr_copy[MAXN];

struct SegmentTree {
    int tree[MAXN * 4];
    int lazy[MAXN * 4];
    
    void apply(int node, int start, int end, int val) {
        tree[node] = (tree[node] + (end - start + 1) * val) % MOD;
        lazy[node] = (lazy[node] + val) % MOD;
    }

    void push(int node, int start, int end) {
        if (lazy[node]) {
            int mid = (start + end) / 2;
            apply(node * 2, start, mid, lazy[node]);
            apply(node * 2 + 1, mid + 1, end, lazy[node]);
            lazy[node] = 0;
        }
    }

    void update(int node, int start, int end, int l, int r, int val) {
        if (l > end || r < start) return;
        if (l <= start && end <= r) {
            apply(node, start, end, val);
            return;
        }
        push(node, start, end);
        int mid = (start + end) / 2;
        update(node * 2, start, mid, l, r, val);
        update(node * 2 + 1, mid + 1, end, l, r, val);
        tree[node] = (tree[node * 2] + tree[node * 2 + 1]) % MOD;
    }

    int query(int node, int start, int end, int l, int r) {
        if (l > end || r < start) return 0;
        if (l <= start && end <= r) return tree[node];
        push(node, start, end);
        int mid = (start + end) / 2;
        return (query(node * 2, start, mid, l, r) + query(node * 2 + 1, mid + 1, end, l, r)) % MOD;
    }
} seg_tree;

void dfs_decompose(int u, int p) {
    parent[u] = p;
    subtree_size[u] = 1;
    int max_sub_size = 0;
    for (const auto& edge : adj[u]) {
        int v = edge.to;
        if (v == p) continue;
        depth[v] = depth[u] + 1;
        dfs_decompose(v, u);
        subtree_size[u] += subtree_size[v];
        if (subtree_size[v] > max_sub_size) {
            max_sub_size = subtree_size[v];
            heavy_child[u] = v;
        }
    }
}

void dfs_build(int u, int h) {
    head[u] = h;
    position[u] = ++current_pos;
    base_array[current_pos] = arr_copy[u];
    if (heavy_child[u]) {
        dfs_build(heavy_child[u], h);
    }
    for (const auto& edge : adj[u]) {
        int v = edge.to;
        if (v == parent[u] || v == heavy_child[u]) continue;
        dfs_build(v, v);
    }
}

void path_update(int u, int v, int val) {
    while (head[u] != head[v]) {
        if (depth[head[u]] < depth[head[v]]) std::swap(u, v);
        seg_tree.update(1, 1, current_pos, position[head[u]], position[u], val);
        u = parent[head[u]];
    }
    if (depth[u] > depth[v]) std::swap(u, v);
    seg_tree.update(1, 1, current_pos, position[u], position[v], val);
}

int path_query(int u, int v) {
    int result = 0;
    while (head[u] != head[v]) {
        if (depth[head[u]] < depth[head[v]]) std::swap(u, v);
        result = (result + seg_tree.query(1, 1, current_pos, position[head[u]], position[u])) % MOD;
        u = parent[head[u]];
    }
    if (depth[u] > depth[v]) std::swap(u, v);
    result = (result + seg_tree.query(1, 1, current_pos, position[u], position[v])) % MOD;
    return result;
}