CSP-Sの点数を上げるためのテンプレート集

数学

高速累乗

int pow_mod(int base, int exp) {
    int result = 1;
    while (exp > 0) {
        if (exp & 1) result = (result * base) % MOD;
        base = (base * base) % MOD;
        exp >>= 1;
    }
    return result;
}

ルーカスの定理（逆元の線形計算付き）

int inv[N], fact[N], fact_inv[N];
void prepare() {
    inv[1] = 1;
    for (int i = 2; i < N; i++) inv[i] = (MOD - MOD / i) * inv[MOD % i] % MOD;
    fact_inv[0] = fact_inv[1] = 1;
    for (int i = 2; i < N; i++) fact_inv[i] = fact_inv[i - 1] * inv[i] % MOD;
    fact[0] = fact[1] = 1;
    for (int i = 2; i < N; i++) fact[i] = fact[i - 1] * i % MOD;
}
int comb(int n, int r) {
    if (r > n) return 0;
    return fact[n] * fact_inv[n - r] % MOD * fact_inv[r] % MOD;
}
int lucas(int n, int r) {
    if (r == 0) return 1;
    return comb(n % MOD, r % MOD) * lucas(n / MOD, r / MOD) % MOD;
}

拡張ユークリッドアルゴリズム

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

拡張ユークリッドによる逆元計算

long long mod_inverse(long long a, long long mod) {
    long long x, y;
    long long gcd = extended_gcd(a, mod, x, y);
    return (gcd == 1) ? (x % mod + mod) % mod : -1;
}

フェルマの小定理による逆元

long long mod_inverse(long long a, long long mod) {
    return power(a, mod - 2, mod);
}

線形逆元の計算

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

CRT（中国剰余定理）

long long chinese_remainder_theorem(int k, long long *a, long long *r) {
    long long product = 1, result = 0;
    for (int i = 1; i <= k; i++) product *= r[i];
    for (int i = 1; i <= k; i++) {
        long long m = product / r[i], b, y;
        extended_gcd(m, r[i], b, y);
        result = (result + a[i] * m * b % product) % product;
    }
    return (result % product + product) % product;
}

拡張CRT

long long extended_chinese_remainder() {
    long long mod = b[1];
    long long res = a[1];
    for (long long i = 2; i <= n; ++i) {
        long long diff = (a[i] - res % b[i] + b[i]) % b[i];
        long long g = extended_gcd(mod, b[i], x, y);
        long long p = b[i] / g;
        x = (x * diff / g % p + p) % p;
        res += x * mod;
        mod *= p;
        res = (res % mod + mod) % mod;
    }
    return res;
}

オイラーの篩

bool is_prime[N];
int primes[N], phi[N], count = 1;
void sieve() {
    phi[1] = 1;
    for (int i = 2; i < N; i++) {
        if (!is_prime[i]) {
            primes[count++] = i;
            phi[i] = i - 1;
        }
        for (int j = 1; j < count && i * primes[j] < N; j++) {
            is_prime[i * primes[j]] = true;
            if (i % primes[j] == 0) {
                phi[i * primes[j]] = phi[i] * primes[j];
                break;
            } else {
                phi[i * primes[j]] = phi[i] * phi[primes[j]];
            }
        }
    }
}

ミラーラビン素数判定

bool miller_rabin(long long n, int rounds) {
    if (n == 2) return true;
    if (n < 2 || (n & 1) == 0) return false;
    long long u = n - 1, x;
    int t = 0;
    while ((u & 1) == 0) {
        t++;
        u >>= 1;
    }
    for (int i = 0; i < rounds; i++) {
        long long a = rand() % (n - 1) + 1;
        x = power(a, u, n);
        for (int j = 0; j < t; j++) {
            long long y = (x * x) % n;
            if (y == 1 && x != 1 && x != n - 1) return false;
            x = y;
        }
        if (x != 1) return false;
    }
    return true;
}

オイラー関数（単体）

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

除算ブロック

for (int left = 1, right; left <= n; left = right + 1) {
    right = n / (n / left);
    cout << (n / left) << ' ' << (right - left + 1) << endl;
}

スキャンライン

#include <bits/stdc++.h>
using namespace std;
const int MAXN = 1e6 + 10;
typedef long long ll;

int n, cnt = 0;
ll x1, y1, x2, y2, X[MAXN << 1];

struct ScanLine {
    ll l, r, h;
    int mark;
    bool operator < (const ScanLine &rhs) const {
        return h < rhs.h;
    }
} lines[MAXN << 1];

struct SegmentTree {
    int l, r, sum;
    ll len;
} tree[MAXN << 2];

void build_tree(int x, int l, int r) {
    tree[x].l = l;
    tree[x].r = r;
    tree[x].len = 0;
    tree[x].sum = 0;
    if (l == r) return;
    int mid = (l + r) >> 1;
    build_tree(x << 1, l, mid);
    build_tree(x << 1 | 1, mid + 1, r);
}

void pushup(int x) {
    if (tree[x].sum) {
        tree[x].len = X[tree[x].r + 1] - X[tree[x].l];
    } else {
        tree[x].len = tree[x << 1].len + tree[x << 1 | 1].len;
    }
}

void update_tree(int x, ll L, ll R, int c) {
    int l = tree[x].l, r = tree[x].r;
    if (X[r + 1] <= L || R <= X[l]) return;
    if (L <= X[l] && X[r + 1] <= R) {
        tree[x].sum += c;
        pushup(x);
        return;
    }
    update_tree(x << 1, L, R, c);
    update_tree(x << 1 | 1, L, R, c);
    pushup(x);
}

int main() {
    scanf("%d", &n);
    for (int i = 1; i <= n; i++) {
        scanf("%lli %lli %lli %lli", &x1, &y1, &x2, &y2);
        X[2 * i - 1] = x1;
        X[2 * i] = x2;
        lines[2 * i - 1] = (ScanLine) {x1, x2, y1, 1};
        lines[2 * i] = (ScanLine) {x1, x2, y2, -1};
    }
    n <<= 1;
    sort(lines + 1, lines + n + 1);
    sort(X + 1, X + n + 1);
    int total = unique(X + 1, X + n + 1) - X - 1;
    build_tree(1, 1, total - 1);
    ll area = 0;
    for (int i = 1; i < n; i++) {
        update_tree(1, lines[i].l, lines[i].r, lines[i].mark);
        area += tree[1].len * (lines[i + 1].h - lines[i].h);
    }
    printf("%lli", area);
    return 0;
}

ガウス消去法

#include <bits/stdc++.h>
using namespace std;
const int N = 105;
int n;
double matrix[N][N], epsilon = 1e-6;

int main() {
    cin >> n;
    for (int i = 1; i <= n; i++) {
        for (int j = 1; j <= n + 1; j++) {
            cin >> matrix[i][j];
        }
    }
    int current_row = 1;
    for (int col = 1; col <= n; col++) {
        int pivot_row = current_row;
        for (int row = current_row; row <= n; row++) {
            if (abs(matrix[row][col]) > abs(matrix[pivot_row][col])) {
                pivot_row = row;
            }
        }
        if (abs(matrix[pivot_row][col]) < epsilon) continue;
        swap(matrix[current_row], matrix[pivot_row]);
        for (int j = n + 1; j >= col; j--) {
            matrix[current_row][j] /= matrix[current_row][col];
        }
        for (int row = current_row + 1; row <= n; row++) {
            for (int j = n + 1; j >= col; j--) {
                matrix[row][j] -= matrix[current_row][j] * matrix[row][col];
            }
        }
        current_row++;
    }
    if (current_row <= n) {
        for (int i = current_row; i <= n; i++) {
            if (abs(matrix[i][n + 1]) > epsilon) {
                cout << -1;
                return 0;
            }
        }
        cout << 0;
        return 0;
    }
    for (int i = n; i >= 1; i--) {
        for (int j = i + 1; j <= n; j++) {
            matrix[i][n + 1] -= matrix[i][j] * matrix[j][n + 1];
        }
    }
    for (int i = 1; i <= n; i++) {
        cout << "x" << i << "=" << matrix[i][n + 1] << '\n';
    }
}

データ構造

単調スタック

void push(int value) {
    while (!stack.empty() && stack.top() < value) stack.pop();
    stack.push(value);
}

単調キュー

deque<int> dq;
for (int i = 1; i <= n; i++) {
    int input;
    cin >> input;
    while (!dq.empty() && dq.back() > input) dq.pop_back();
    dq.push_back(input);
}

01トライ

int trie[N * 31][2], node_count = 1;
void insert(int pos, int depth, int value) {
    bool bit = (value & (1 << depth));
    if (!trie[pos][bit]) trie[pos][bit] = ++node_count;
    if (depth) insert(trie[pos][bit], depth - 1, value);
}
int query(int pos, int depth, int value) {
    if (depth < 0) return 0;
    bool bit = (value & (1 << depth));
    if (trie[pos][!bit]) return (1 << depth) + query(trie[pos][!bit], depth - 1, value);
    if (trie[pos][bit]) return query(trie[pos][bit], depth - 1, value);
    return 0;
}

ST表

#include <bits/stdc++.h>
using namespace std;
const int LOGN = 21;
const int MAXN = 2000001;
int st_table[MAXN][LOGN + 1], log_values[MAXN + 1];

int read() {
    char c = getchar();
    int x = 0, f = 1;
    while (c < '0' || c > '9') {
        if (c == '-') f = -1;
        c = getchar();
    }
    while (c >= '0' && c <= '9') {
        x = x * 10 + c - '0';
        c = getchar();
    }
    return x * f;
}

void preprocess() {
    log_values[1] = 0;
    log_values[2] = 1;
    for (int i = 3; i < MAXN; i++) {
        log_values[i] = log_values[i / 2] + 1;
    }
}

int main() {
    int n = read(), m = read();
    for (int i = 1; i <= n; i++) st_table[i][0] = read();
    preprocess();
    for (int j = 1; j <= LOGN; j++) {
        for (int i = 1; i + (1 << j) - 1 <= n; i++) {
            st_table[i][j] = max(st_table[i][j - 1], st_table[i + (1 << (j - 1))][j - 1]);
        }
    }
    for (int i = 1; i <= m; i++) {
        int x = read(), y = read();
        int s = log_values[y - x + 1];
        printf("%d\n", max(st_table[x][s], st_table[y - (1 << s) + 1][s]));
    }
    return 0;
}

フィックステータス木

#define low_bit(x) ((x) & -(x))
using namespace std;
int n, m, tree[200000 + 10];
void update(int index, int delta) {
    for (int i = index; i <= n; i += low_bit(i)) tree[i] += delta;
}
int query(int range) {
    int result = 0;
    for (int i = range; i >= 1; i -= low_bit(i)) result += tree[i];
    return result;
}

線分木

void add_tag(int node, int value) {
    seg[node].tag += value;
    seg[node].sum += (seg[node].r - seg[node].l + 1) * value;
}
void push_up(int node) {
    seg[node].sum = seg[left_child(node)].sum + seg[right_child(node)].sum;
}
void push_down(int node) {
    add_tag(left_child(node), seg[node].tag);
    add_tag(right_child(node), seg[node].tag);
    seg[node].tag = 0;
}
void build(int node, int left, int right) {
    seg[node].l = left;
    seg[node].r = right;
    if (left == right) {
        seg[node].sum = arr[left];
        seg[node].tag = 0;
        return;
    }
    int mid = left + (right - left) / 2;
    build(left_child(node), left, mid);
    build(right_child(node), mid + 1, right);
    push_up(node);
}
void modify(int node, int left, int right, int value) {
    if (left <= seg[node].l && seg[node].r <= right) {
        add_tag(node, value);
        return;
    }
    push_down(node);
    if (left <= seg[left_child(node)].r) modify(left_child(node), left, right, value);
    if (seg[right_child(node)].l <= right) modify(right_child(node), left, right, value);
    push_up(node);
}
int query(int node, int left, int right) {
    if (left <= seg[node].l && seg[node].r <= right) return seg[node].sum;
    push_down(node);
    int sum = 0;
    if (left <= seg[left_child(node)].r) sum += query(left_child(node), left, right);
    if (seg[right_child(node)].l <= right) sum += query(right_child(node), left, right);
    return sum;
}

FHQ Treap

mt19937 rng(time(0));
const int MAX_N = 1e6 + 10;
class Node {
public:
    int lson = 0, rson = 0;
    int value, priority, size = 0;
    Node() {
        lson = rson = size = 0;
    }
};
class Treap {
public:
    int size, root;
    Node nodes[MAX_N];
    int new_node(int val) {
        nodes[++size].value = val;
        nodes[size].priority = rng();
        nodes[size].lson = nodes[size].rson = 0;
        nodes[size].size = 1;
        return size;
    }
    Treap() {
        size = 0;
        root = 0;
    }
    void push_up(int node) {
        nodes[node].size = nodes[nodes[node].lson].size + nodes[nodes[node].rson].size + 1;
    }
    void split(int node, int key, int &left, int &right) {
        if (!node) {
            left = right = 0;
            return;
        } else if (nodes[node].value <= key) {
            left = node;
            split(nodes[node].rson, key, nodes[node].rson, right);
        } else {
            right = node;
            split(nodes[node].lson, key, left, nodes[node].lson);
        }
        push_up(node);
    }
    int merge(int left, int right) {
        if (!left || !right) return left + right;
        if (nodes[left].priority > nodes[right].priority) {
            nodes[left].rson = merge(nodes[left].rson, right);
            push_up(left);
            return left;
        } else {
            nodes[right].lson = merge(left, nodes[right].lson);
            push_up(right);
            return right;
        }
    }
    void insert(int val) {
        int left, right;
        split(root, val, left, right);
        root = merge(merge(left, new_node(val)), right);
    }
    void remove(int val) {
        int left, middle, right;
        split(root, val, left, right);
        split(left, val - 1, left, middle);
        middle = merge(nodes[middle].lson, nodes[middle].rson);
        root = merge(merge(left, middle), right);
    }
    int rank(int val) {
        int left, right;
        split(root, val - 1, left, right);
        int result = nodes[left].size;
        root = merge(left, right);
        return result + 1;
    }
    int get_value_by_rank(int rnk) {
        int current = root;
        while (current) {
            if (nodes[nodes[current].lson].size + 1 == rnk) break;
            else if (nodes[nodes[current].lson].size + 1 > rnk)
                current = nodes[current].lson;
            else {
                rnk -= nodes[nodes[current].lson].size + 1;
                current = nodes[current].rson;
            }
        }
        return nodes[current].value;
    }
    int predecessor(int val) {
        int left, right;
        split(root, val - 1, left, right);
        int current = left;
        while (nodes[current].rson) current = nodes[current].rson;
        int result = nodes[current].value;
        root = merge(left, right);
        return result;
    }
    int successor(int val) {
        int left, right;
        split(root, val, left, right);
        int current = right;
        while (nodes[current].lson) current = nodes[current].lson;
        int result = nodes[current].value;
        root = merge(left, right);
        return result;
    }
};

除算ブロック

for (int left = 1, right; left <= n; left = right + 1) {
    right = n / (n / left);
    cout << (n / left) << ' ' << (right - left + 1) << endl;
}

文字列

AC自動機（/KMP）

struct TrieNode {
    int children[26], fail, flag, count;
    void clear() {
        memset(children, 0, sizeof(children));
        fail = flag = count = 0;
    }
} trie[MAXN];
queue<int> q;

void insert(char* str, int id) {
    int current = 1, length = strlen(str);
    for (int i = 0; i < length; i++) {
        int ch = str[i] - 'a';
        if (!trie[current].children[ch]) trie[current].children[ch] = ++node_count;
        current = trie[current].children[ch];
    }
    if (!trie[current].flag) trie[current].flag = id;
    mapping[id] = trie[current].flag;
}

void compute_fail() {
    for (int i = 0; i < 26; i++) trie[0].children[i] = 1;
    q.push(1);
    while (!q.empty()) {
        int current = q.front();
        q.pop();
        int failure = trie[current].fail;
        for (int i = 0; i < 26; i++) {
            int child = trie[current].children[i];
            if (!child) {
                trie[current].children[i] = trie[failure].children[i];
                continue;
            }
            trie[child].fail = trie[failure].children[i];
            in[trie[child].fail]++;
            q.push(child);
        }
    }
}

void topological_sort() {
    for (int i = 1; i <= node_count; i++) {
        if (in[i] == 0) q.push(i);
    }
    while (!q.empty()) {
        int current = q.front();
        q.pop();
        visited[trie[current].flag] = trie[current].count;
        int parent = trie[current].fail;
        in[parent]--;
        trie[parent].count += trie[current].count;
        if (in[parent] == 0) q.push(parent);
    }
}

void search(char* text) {
    int current = 1, length = strlen(text);
    for (int i = 0; i < length; i++) {
        current = trie[current].children[text[i] - 'a'];
        trie[current].count++;
    }
}

文字列ハッシュ

typedef unsigned long long ull;
const ull SEED_POOL[] = {146527, 19260817};
const ull MOD_POOL[] = {1000000009, 998244353};

struct Hash {
    ull seed, mod;
    vector<ull> power, prefix_sum;
    
    Hash() {}
    
    Hash(const string& str, const int seed_idx, const int mod_idx) {
        seed = SEED_POOL[seed_idx];
        mod = MOD_POOL[mod_idx];
        int len = str.length();
        power.resize(len + 1);
        prefix_sum.resize(len + 1);
        power[0] = 1;
        for (int i = 1; i <= len; i++) power[i] = power[i - 1] * seed % mod;
        for (int i = 1; i <= len; i++) prefix_sum[i] = (prefix_sum[i - 1] * seed % mod + str[i - 1]) % mod;
    }
    
    ull get(int start, int end) {
        return (prefix_sum[end] - prefix_sum[start] * power[end - start] % mod + mod) % mod;
    }
    
    ull substring(int start, int length) {
        return get(start, start + length);
    }
};

グラフ理論

隣接リスト表現

struct Edge {
    int to, weight, next;
} edges[M * 2];
int head[N], edge_count;

void add_edge(int from, int to, int weight) {
    edges[++edge_count] = (Edge){to, weight, head[from]};
    head[from] = edge_count;
}

Dijkstra

struct Edge {
    int to, weight;
};
vector<Edge> graph[N];
struct Node {
    int vertex, distance;
    bool operator<(const Node& other) const {
        return distance > other.distance;
    }
};
priority_queue<Node> pq;
int distances[N];
bool visited[N];

void dijkstra(int start) {
    memset(distances, 0x3f, sizeof distances);
    pq.push({start, 0});
    distances[start] = 0;
    while (!pq.empty()) {
        int current = pq.top().vertex;
        pq.pop();
        if (visited[current]) continue;
        visited[current] = true;
        for (int i = 0; i < graph[current].size(); i++) {
            int neighbor = graph[current][i].to;
            int cost = distances[current] + graph[current][i].weight;
            if (distances[neighbor] > cost) {
                distances[neighbor] = cost;
                pq.push({neighbor, cost});
            }
        }
    }
}

木の分割（重み付き）

void dfs(int node, int parent) {
    subtree_size[node] = 1;
    for (int &child : children[node]) {
        if (child != parent) {
            parent_of[child] = node;
            depth[child] = depth[node] + 1;
            dfs(child, node);
            subtree_size[node] += subtree_size[child];
            if (subtree_size[child] >= subtree_size[heavy_child[node]])
                heavy_child[node] = child, swap(child, children[node][0]);
        }
    }
}

void dfs2(int node, int parent, int chain_head) {
    dfn[node] = ++counter;
    position[dfn[node]] = node;
    chain_head[node] = chain_head;
    for (int &child : children[node]) {
        if (child != parent) {
            if (child == heavy_child[node])
                dfs2(child, node, chain_head);
            else
                dfs2(child, node, child);
        }
    }
}

int lca(int a, int b) {
    while (chain_head[a] != chain_head[b]) {
        if (depth[chain_head[b]] > depth[chain_head[a]])
            swap(a, b);
        a = parent_of[chain_head[a]];
    }
    return depth[a] > depth[b] ? b : a;
}

倍増LCA

void dfs(int root, int parent) {
    ancestor[root][0] = parent;
    depth[root] = depth[parent] + 1;
    for (int i = 1; i < 31; ++i) {
        ancestor[root][i] = ancestor[ancestor[root][i - 1]][i - 1];
        cost[root][i] = cost[ancestor[root][i - 1]][i - 1] + cost[root][i - 1];
    }
    for (auto child : children[root]) {
        if (child == parent) continue;
        cost[child][0] = weight[root][child];
        dfs(child, root);
    }
}

int find_lca(int a, int b) {
    if (depth[a] > depth[b]) swap(a, b);
    int diff = depth[b] - depth[a], total_cost = 0;
    for (int j = 0; diff; ++j, diff >>= 1)
        if (diff & 1) total_cost += cost[b][j], b = ancestor[b][j];
    if (b == a) return total_cost;
    for (int j = 30; j >= 0 && b != a; --j) {
        if (ancestor[a][j] != ancestor[b][j]) {
            total_cost += cost[a][j] + cost[b][j];
            a = ancestor[a][j];
            b = ancestor[b][j];
        }
    }
    total_cost += cost[a][0] + cost[b][0];
    return total_cost;
}

強連結成分分解（Tarjan）

void tarjan(int node) {
    dfn[node] = low[node] = ++counter;
    visited[node] = true;
    stack.push(node);
    for (auto neighbor : graph[node]) {
        if (!dfn[neighbor]) {
            tarjan(neighbor);
            low[node] = min(low[node], low[neighbor]);
        } else if (visited[neighbor]) {
            low[node] = min(low[node], dfn[neighbor]);
        }
    }
    if (low[node] == dfn[node]) {
        int current_node;
        do {
            current_node = stack.top();
            visited[current_node] = false;
            component[current_node] = component_count;
            components[component_count].push_back(current_node);
            stack.pop();
        } while (current_node != node);
        component_count++;
    }
}

トポロジカルソート

bool topological_sort() {
    vector<int> result;
    queue<int> zero_in_degree;
    for (int i = 1; i <= n; i++) {
        if (in_degree[i] == 0) zero_in_degree.push(i);
    }
    while (!zero_in_degree.empty()) {
        int current = zero_in_degree.front();
        zero_in_degree.pop();
        result.push_back(current);
        for (auto neighbor : graph[current]) {
            if (--in_degree[neighbor] == 0) {
                zero_in_degree.push(neighbor);
            }
        }
    }
    if (result.size() == n) {
        for (auto i : result) cout << i << ' ';
        return true;
    } else {
        return false;
    }
}

最小生成木（Kruskal）

for (int i = 1; i <= n; i++) parent[i] = i;
while (!edges.empty()) {
    auto [weight, u, v] = edges.top();
    edges.pop();
    if (find_parent(u) == find_parent(v)) continue;
    if (!(n--) || !(k--)) break;
    union_sets(u, v);
    max_weight = max(max_weight, weight);
}

Dinic最大流

int n, m, source, sink, level[N], head[N], current[N];
struct Edge {
    int to, capacity, next;
} edges[M * 2];
int edge_count = 1;

void add_edge(int from, int to, int capacity) {
    edges[++edge_count] = (Edge){to, capacity, head[from]};
    head[from] = edge_count;
}

void link_edge(int from, int to, int capacity) {
    add_edge(from, to, capacity);
    add_edge(to, from, 0);
}

bool bfs() {
    queue<int> q;
    for (int i = source; i <= sink; i++) level[i] = -1;
    level[source] = 0;
    q.push(source);
    while (!q.empty()) {
        int current = q.front();
        current_flow[current] = head[current];
        q.pop();
        for (int i = head[current]; i; i = edges[i].next) {
            if (level[edges[i].to] == -1 && edges[i].capacity) {
                level[edges[i].to] = level[current] + 1;
                q.push(edges[i].to);
            }
        }
    }
    return level[sink] != -1;
}

int dfs(int node, int flow) {
    if (node == sink) return flow;
    int total_flow = 0;
    for (int i = current[node]; i && flow; i = edges[i].next) {
        current[node] = i;
        if (level[edges[i].to] == level[node] + 1 && edges[i].capacity) {
            int result = dfs(edges[i].to, min(flow, edges[i].capacity));
            if (result > 0) {
                edges[i].capacity -= result;
                edges[i ^ 1].capacity += result;
                total_flow += result;
                flow -= result;
            }
        }
    }
    return total_flow;
}

最小費用流（MCMF）

int n, m, source, sink, dist[N], head[N], current[N], total_cost, total_flow;
bool visited[N];
struct Edge {
    int to, capacity, cost, next;
} edges[M * 2];
int edge_count = 1;

void add_edge(int from, int to, int capacity, int cost) {
    edges[++edge_count] = (Edge){to, capacity, cost, head[from]};
    head[from] = edge_count;
}

void link_edge(int from, int to, int capacity, int cost) {
    add_edge(from, to, capacity, cost);
    add_edge(to, from, 0, -cost);
}

bool spfa() {
    for (int i = source; i <= sink; i++) dist[i] = 1e9;
    for (int i = source; i <= sink; i++) visited[i] = false;
    queue<int> q;
    q.push(source);
    visited[source] = true;
    dist[source] = 0;
    while (!q.empty()) {
        int current = q.front();
        current_flow[current] = head[current];
        q.pop();
        visited[current] = false;
        for (int i = head[current]; i; i = edges[i].next) {
            if (dist[edges[i].to] > dist[current] + edges[i].cost && edges[i].capacity) {
                dist[edges[i].to] = dist[current] + edges[i].cost;
                if (!visited[edges[i].to]) {
                    q.push(edges[i].to);
                    visited[edges[i].to] = true;
                }
            }
        }
    }
    return dist[sink] != 1e9;
}

int dfs(int node, int flow) {
    if (node == sink) return flow;
    visited[node] = true;
    int total_flow = 0;
    for (int i = current[node]; i && flow; i = edges[i].next) {
        current[node] = i;
        if (!visited[edges[i].to] && dist[edges[i].to] == dist[node] + edges[i].cost && edges[i].capacity) {
            int result = dfs(edges[i].to, min(flow, edges[i].capacity));
            edges[i].capacity -= result;
            edges[i ^ 1].capacity += result;
            total_flow += result;
            flow -= result;
            total_cost += result * edges[i].cost;
        }
    }
    visited[node] = false;
    return total_flow;
}

その他の手法

行列の高速累乗

typedef vector<long long> vec;
typedef vector<vec> mat;

mat multiply(mat& A, mat& B) {
    mat C(A.size(), vec(B[0].size()));
    for (int i = 0; i < A.size(); i++)
        for (int k = 0; k < B.size(); k++)
            if (A[i][k])
                for (int j = 0; j < B[0].size(); j++)
                    C[i][j] = (C[i][j] + A[i][k] * B[k][j]) % MOD;
    return C;
}

mat power(mat A, long long n) {
    mat result(A.size(), vec(A.size()));
    for (int i = 0; i < A.size(); i++) result[i][i] = 1;
    for (; n; n >>= 1, A = multiply(A, A))
        if (n & 1) result = multiply(result, A);
    return result;
}