Mathematical Algorithms

Fast Exponentiation

Reduces the time complexity of computing powers to logarithmic time by leveraging binary decomposition of the exponent.

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

Lucas Theorem

Efficiently computes combinations $\binom{n}{k} \pmod p$ when $p$ is a small prime and $n, k$ are large.

const int MAXN = 1e5 + 5;
const int MOD = 1e9 + 7; // Must be prime
long long fact[MAXN], invFact[MAXN];

void precompute_combinatorics() {
    fact[0] = 1;
    invFact[0] = 1;
    for (int i = 1; i < MAXN; ++i) {
        fact[i] = (fact[i - 1] * i) % MOD;
        invFact[i] = power(fact[i], MOD - 2);
    }
}

long long comb(int n, int m) {
    if (m < 0 || m > n) return 0;
    return fact[n] * invFact[n - m] % MOD * invFact[m] % MOD;
}

long long lucas(long long n, long long m) {
    if (m == 0) return 1;
    return comb(n % MOD, m % MOD) * lucas(n / MOD, m / MOD) % MOD;
}

Extended Euclidean Algorithm

Finds integers $x$ and $y$ such that $ax + by = \gcd(a, b)$. Critical for solving linear Diophantine equations and modular inverses.

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 g = extended_gcd(b, a % b, y, x);
    y -= (a / b) * x;
    return g;
}

Modular Inverse Methods

Two standard approaches: one relies on the Extended Euclidean Algorithm (works for any modulus), the other uses Fermat's Little Theorem (requires prime modulus).

// Using Extended GCD
long long mod_inverse_extgcd(long long a, long long m) {
    long long x, y;
    long long g = extended_gcd(a, m, x, y);
    if (g != 1) return -1; // No inverse exists
    return (x % m + m) % m;
}

// Using Fermat's Little Theorem (m must be prime)
long long mod_inverse_fermat(long long a, long long p) {
    return binary_pow(a, p - 2, p);
}

Linear Inverse Precomputation

Calculates modular inverses for all integers up to $N$ in linear time using the recurrence: $inv[i] = (MOD - \lfloor MOD/i \rfloor) \times inv[MOD \% i] \pmod{MOD}$.

const int MAXN = 1e6 + 5;
long long inv[MAXN];

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

Chinese Remainder Theorem (CRT)

Solves a system of congruences with pairwise coprime moduli:

$x \equiv a_i \pmod{m_i}$

long long crt_solver(const vector<long long="">& remainders, const vector<long long="">& moduli) {
    long long lcm_all = 1;
    for (long long m : moduli) lcm_all *= m;
    
    long long total_solution = 0;
    for (size_t i = 0; i < moduli.size(); ++i) {
        long long mi = moduli[i];
        long long Yi = lcm_all / mi;
        long long yi_inv = mod_inverse_extgcd(Yi, mi);
        total_solution = (total_solution + remainder[i] * Yi % lcm_all * yi_inv % lcm_all) % lcm_all;
    }
    return (total_solution + lcm_all) % lcm_all;
}</long></long>

Extended CRT (Non-Coprime Moduli)

Merges congruence equations iteratively to handle cases where moduli are not necessarily coprime.

long long excrt(const vector<long long="">& a, const vector<long long="">& m) {
    long long ans = a[0], M = m[0];
    for (size_t i = 1; i < a.size(); ++i) {
        long long x, y;
        long long g = extended_gcd(M, m[i], x, y);
        long long diff = ((a[i] - ans) % m[i] + m[i]) % m[i];
        if (diff % g != 0) return -1; // No solution
        
        long long t = (diff / g) % (m[i] / g);
        x = (x * t) % (m[i] / g);
        
        ans += M * x;
        M *= (m[i] / g);
        ans = (ans % M + M) % M;
    }
    return ans;
}</long></long>

Euler's Totient Sieve

Computes $\phi(n)$ for all $n \leq N$ in $O(N)$ time using a linear sieve.

const int MAXN = 1e6 + 5;
int phi[MAXN];
bool is_composite[MAXN];
vector<int> primes;

void euler_sieve() {
    phi[1] = 1;
    for (int i = 2; i < MAXN; ++i) {
        if (!is_composite[i]) {
            primes.push_back(i);
            phi[i] = i - 1;
        }
        for (int p : primes) {
            if (1LL * i * p >= MAXN) break;
            is_composite[i * p] = true;
            if (i % p == 0) {
                phi[i * p] = phi[i] * p;
                break;
            } else {
                phi[i * p] = phi[i] * (p - 1);
            }
        }
    }
}</int>

Miller-Rabin Primality Test

A probabilistic algorithm optimized for competitive programming. Uses a deterministic set of bases for inputs within standard integer ranges.

bool miller_rabin(long long n) {
    if (n < 2) return false;
    if (n == 2 || n == 3) return true;
    if (n % 2 == 0) return false;

    long long d = n - 1, r = 0;
    while (d % 2 == 0) { d /= 2; r++; }

    // Deterministic bases for n < 3,317,044,064,279,371
    static const vector<long long=""> bases = {2, 3, 5, 7, 11, 13, 17, 19, 23};
    
    for (long long a : bases) {
        if (a >= n) break;
        long long x = binary_pow(a, d, n);
        if (x == 1 || x == n - 1) continue;
        bool composite = true;
        for (int i = 0; i < r - 1; ++i) {
            x = (x * x) % n;
            if (x == n - 1) { composite = false; break; }
        }
        if (composite) return false;
    }
    return true;
}</long>

Single Integer Euler Function

Computes $\phi(n)$ for a single number via trial division up to $\sqrt{n}$.

long long get_phi_single(long long n) {
    long long result = n;
    for (long long i = 2; i * i <= n; ++i) {
        if (n % i == 0) {
            while (n % i == 0) n /= i;
            result -= result / i;
        }
    }
    if (n > 1) result -= result / n;
    return result;
}

Number Theory Block (Division Block)

Optimizes summation loops where terms like $\lfloor n/i \rfloor$ remain constant over contiguous ranges.

for (long long l = 1, r; l <= n; l = r + 1) {
    r = n / (n / l);
    // Process range [l, r] where floor(n/k) equals floor(n/l)
    long long term_value = n / l;
    long long range_count = r - l + 1;
    // Accumulate results...
}

Data Structure Algorithms

Monotonic Stack

Maintains elements in increasing or decreasing order. Efficiently solves next-greater-element and histogram problems.

vector<int> monotonic_stack(const vector<int>& arr, bool ascending = true) {
    vector<int> res;
    stack<int> s;
    for (int val : arr) {
        while (!s.empty() && (ascending ? s.top() > val : s.top() < val)) {
            res.push_back(s.top());
            s.pop();
        }
        s.push(val);
    }
    while (!s.empty()) {
        res.push_back(s.top());
        s.pop();
    }
    return res;
}</int></int></int></int>

Monotonic Queue

Implements a sliding window minimum/maximum using a deque. Time complexity remains linear per element.

deque<int> mono_queue_max(const vector<int>& nums, int k) {
    deque<int> dq;
    vector<int> result;
    for (int i = 0; i < nums.size(); ++i) {
        while (!dq.empty() && nums[dq.back()] <= nums[i]) dq.pop_back();
        while (!dq.empty() && dq.front() <= i - k) dq.pop_front();
        dq.push_back(i);
        if (i >= k - 1) result.push_back(nums[dq.front()]);
    }
    return result;
}</int></int></int></int>

Binary Trie (0-1 Trie)

Stores binary representations of numbers. Enables fast retrieval of values that maximize/minimize XOR with a query.

struct TrieNode {
    TrieNode* children[2] = {};
};

class BinaryTrie {
public:
    void insert(int num) {
        TrieNode* curr = root;
        for (int i = 30; i >= 0; --i) {
            int bit = (num >> i) & 1;
            if (!curr->children[bit]) curr->children[bit] = new TrieNode();
            curr = curr->children[bit];
        }
    }

    int max_xor(int num) {
        TrieNode* curr = root;
        int result = 0;
        for (int i = 30; i >= 0; --i) {
            int bit = (num >> i) & 1;
            int desired = 1 - bit;
            if (curr->children[desired]) {
                result |= (1 << i);
                curr = curr->children[desired];
            } else {
                curr = curr->children[bit];
            }
        }
        return result;
    }
private:
    TrieNode* root = new TrieNode();
};

Sparse Table (RMQ)

Provides $O(1)$ range maximum/minimum queries after $O(N \log N)$ preprocessing. Suitable for static arrays.

const int LOG = 20;
int st[MAXN][LOG];
int logs[MAXN];

void build_sparse_table(int n) {
    logs[1] = 0;
    for (int i = 2; i <= n; ++i) logs[i] = logs[i / 2] + 1;
    
    for (int i = 1; i <= n; ++i) st[i][0] = arr[i];
    
    for (int j = 1; j <= logs[n]; ++j) {
        for (int i = 1; i + (1 << j) - 1 <= n; ++i) {
            st[i][j] = max(st[i][j - 1], st[i + (1 << (j - 1))][j - 1]);
        }
    }
}

int query_range_max(int L, int R) {
    int j = logs[R - L + 1];
    return max(st[L][j], st[R - (1 << j) + 1][j]);
}

Fenwick Tree (Binary Indexed Tree)

Supports efficient point updates and prefix sum calculations in $O(\log N)$ time.

class FenwickTree {
private:
    vector<int> tree;
    int size;

public:
    FenwickTree(int n) : size(n), tree(n + 1, 0) {}

    void update(int idx, int delta) {
        for (; idx <= size; idx += idx & (-idx))
            tree[idx] += delta;
    }

    int query(int idx) const {
        int sum = 0;
        for (; idx > 0; idx -= idx & (-idx))
            sum += tree[idx];
        return sum;
    }

    int query_range(int l, int r) const {
        return query(r) - query(l - 1);
    }
};</int>

Segment Tree with Lazy Propagation

Handles range updates and range queries dynamically. Essential for problems requiring interval modifications.

struct SegTree {
    int n;
    vector<long long=""> tree, lazy;

    SegTree(int size, const vector<long long="">& initial_arr) : n(size) {
        tree.resize(4 * n);
        lazy.resize(4 * n, 0);
        build(initial_arr, 1, 0, n - 1);
    }

    void push(int node, int l, int r) {
        if (lazy[node] != 0) {
            tree[node] += lazy[node];
            if (l != r) {
                lazy[2 * node] += lazy[node];
                lazy[2 * node + 1] += lazy[node];
            }
            lazy[node] = 0;
        }
    }

    void build(const vector<long long="">& arr, int node, int l, int r) {
        if (l == r) { tree[node] = arr[l]; return; }
        int mid = (l + r) / 2;
        build(arr, 2 * node, l, mid);
        build(arr, 2 * node + 1, mid + 1, r);
        tree[node] = tree[2 * node] + tree[2 * node + 1];
    }

    void range_add(int ql, int qr, long long val, int node, int l, int r) {
        push(node, l, r);
        if (ql > r || qr < l) return;
        if (ql <= l && r <= qr) {
            lazy[node] += val;
            push(node, l, r);
            return;
        }
        int mid = (l + r) / 2;
        range_add(ql, qr, val, 2 * node, l, mid);
        range_add(ql, qr, val, 2 * node + 1, mid + 1, r);
        tree[node] = tree[2 * node] + tree[2 * node + 1];
    }

    long long query(int ql, int qr, int node, int l, int r) {
        push(node, l, r);
        if (ql > r || qr < l) return 0;
        if (ql <= l && r <= qr) return tree[node];
        int mid = (l + r) / 2;
        return query(ql, qr, 2 * node, l, mid) + query(ql, qr, 2 * node + 1, mid + 1, r);
    }
};</long></long></long>