Download Cheatsheet for Competitive Programming and more Cheat Sheet Web Programming and Technologies in PDF only on Docsity!
Cheatsheet for competitive programming
bubblesolve (Elias Riedel Gårding)
March 21, 2020
Contents
1 Setup 1 1.1 Caps Lock as Escape............................... 1 1.2 .vimrc..................................... 1 1.3 Compilation.................................. 1
2 Graph algorithms 1 2.1 DFS...................................... 1 2.2 Dijkstra.................................... 3 2.3 Bellman–Ford.................................. 3 2.4 Floyd–Warshall (all-pairs shortest path).................... 4 2.5 Maximum flow (Ford–Fulkerson).......................... 4 2.6 Minimum spanning tree.............................. 6
3 Data structures 6 3.1 Union-Find................................... 6 3.2 Fenwick Tree.................................. 7
4 String algorithms 7 4.1 KMP...................................... 7
5 Number theory 8 5.1 Combinatorics.................................. 8 5.2 Extended Euclidean algorithm.......................... 8 5.3 Chinese remainder theorem............................ 9
6 Standard problems 9 6.1 Knapsack.................................... 9 6.2 Longest increasing subsequence......................... 10 6.3 Interval cover................................. 10
7 Miscellaneous algorithms 11 7.1 Binary search.................................. 11 7.2 Ternary search (unimodal optimization)..................... 12 7.3 Gaussian elimination.............................. 12 7.4 Simplex algorithm (linear programming)..................... 14 7.5 Basis conversion................................ 16
8 Mathematical objects 16 8.1 Fraction.................................... 16 8.2 Bignum..................................... 18 8.3 Matrix..................................... 21 8.4 Polynomial................................... 22
1 Setup
1.1 Caps Lock as Escape
xmodmap -e "clear Lock" # If that doesn't work, xmodmap -e "remove Lock = Caps_Lock" xmodmap -e "keycode 66 = Escape NoSymbol Escape"
1.2 .vimrc
set nocompatible filetype plugin indent on set autoindent set tabstop=4 shiftwidth=4 expandtab set foldmethod=marker set number relativenumber syntax on
1.3 Compilation
compile() { g++ -O2 -Wall -Wextra -Wfloat-equal -Wunreachable-code -Wfatal-errors
-Wformat=2 -std=gnu++17 -D_GLIBCXX_DEBUG "$1" -o "$(basename "$1" .cpp)" }
run() { compile "$1" "$(realpath "$(basename "$1" .cpp)")" }
2 Graph algorithms
2.1 DFS
#pragma once
#include #include #include #include
using NodeFunc = std::function<void (int)>; using EdgeFunc = std::function<void (int, int)>;
void node_nop(int) {} void edge_nop(int, int) {}
_/* Depth-first search in a directed graph.
- Conceptually, a node has one of three colours:
- · White: untouched
- · Gray: in the queue
- · Black: has been popped from the queue
- This class offers three callbacks:
- · explore(from, to): called when exploring the edge
from -> to - · discover(node): called when
node becomes gray - · finish(node): called when
node becomes black */_ class DFS { public: int const N; std::vector<std::vector> const &adj; std::vector visited; EdgeFunc explore; NodeFunc discover, finish;
DFS(std::vector<std::vector> const &adj) : N(adj.size()), adj(adj), visited(N, false), explore(edge_nop), discover(node_nop), finish(node_nop) {}
void dfs_from(int node) { if (visited[node]) return; visited[node] = true;
discover(node);
for (int child : adj[node]) { explore(node, child); dfs_from(child); }
finish(node); }
void run() { for (int node = 0; node < N; ++node) dfs_from(node); }
DFS &on_explore(EdgeFunc f) { explore = f; return *this; }
DFS &on_discover(NodeFunc f) { discover = f; return *this; }
DFS &on_finish(NodeFunc f) { finish = f; return *this; } };
std::vector topological_sort(std::vector<std::vector> const &adj) _/* Given a directed graph in the form of an adjacency list, return a topological
- ordering of the graph's DFS forest: an ordering of the nodes such that a
- parent always appears before all of its descendants. If the graph is acyclic,
- this is a topological ordering of the graph itself. */_ { std::vector result; result.reserve(adj.size());
DFS(adj).on_finish([&] (int node) { result.push_back(node); }).run();
std::reverse(result.begin(), result.end()); return result; }
std::vector<std::vector> scc_kosaraju( std::vector<std::vector> const &adj) _/* Given a directed graph in the form of an adjacency list, return a vector of
- its strongly connected components (where each component is a list of nodes). */_ { int const N = adj.size();
std::vector order = topological_sort(adj);
std::vector<std::vector> adj_T(N); for (int i = 0; i < N; ++i) for (int j : adj[i]) adj_T[j].push_back(i);
std::vector<std::vector> comps(1); DFS dfs(adj_T); dfs.on_finish([&] (int node) { (comps.end() - 1)->push_back(node); });
for (int node : order) { if (!(comps.end() - 1)->empty()) comps.emplace_back();
dfs.dfs_from(node); }
if ((comps.end() - 1)->empty()) comps.erase(comps.end() - 1);
return comps; }
_* negative_cycle is set if there are negative cycles in the graph.
- Nodes numbered from 0 to N - 1 where N = edges.size().
- Distances are INF if there is no path, -INF if arbitrarily short paths exist.
- If inf_ptr is not null, puts INF into *inf_ptr. */_ { constexpr D INF = std::numeric_limits::max() / 2; if (inf_ptr != nullptr) *inf_ptr = INF; int const N = edges.size();
prev.assign(N, -1); std::vector dist(N, INF); dist[source] = 0;
// N - 1 phases for finding shortest paths, // N phases for finding negative cycles. negative_cycle = false; for (int ph = 0; ph < 2 * N - 1; ++ph) // Iterate over all edges for (int u = 0; u < N; ++u) if (dist[u] < INF) // prevent catching negative INF -> INF edges for (std::pair<int, D> const &edge : edges[u]) { int v; D w; std::tie(v, w) = edge; if (dist[v] > dist[u] + w) { if (ph < N - 1) { dist[v] = dist[u] + w; prev[v] = u; } else { negative_cycle = true; dist[v] = -INF; } } }
return dist; }
2.4 Floyd–Warshall (all-pairs shortest path)
#pragma once
#include #include #include
template std::vector<std::vector> floyd_warshall( std::vector<std::vector<std::pair<int, D>>> const &edges, bool &negative_cycle, D inf_ptr = nullptr) _/ Returns a matrix of the shortest distances between all nodes.
- negative_cycle is set if there are negative cycles in the graph._
_* Nodes numbered from 0 to N - 1 where N = edges.size().
- Distances are INF if there is no path, -INF if arbitrarily short paths exist.
- If inf_ptr is not null, puts INF into *inf_ptr. */_ { constexpr D INF = std::numeric_limits::max() / 2; if (inf_ptr != nullptr) *inf_ptr = INF; int const N = edges.size();
// Initialize distance matrix std::vector<std::vector> dist(N, std::vector(N, INF)); for (int u = 0; u < N; ++u) { dist[u][u] = 0; for (std::pair<int, D> const &edge : edges[u]) { int v; D w; std::tie(v, w) = edge; dist[u][v] = std::min(dist[u][v], w); } }
// Main loop for (int k = 0; k < N; ++k) for (int i = 0; i < N; ++i) for (int j = 0; j < N; ++j) if (dist[i][k] < INF && dist[k][j] < INF) dist[i][j] = std::min(dist[i][j], dist[i][k] + dist[k][j]);
// Propagate negative cycles negative_cycle = false; for (int i = 0; i < N; ++i) for (int j = 0; j < N; ++j) for (int k = 0; k < N; ++k) if (dist[i][k] < INF && dist[k][j] < INF && dist[k][k] < 0) { negative_cycle = true; dist[i][j] = -INF; }
return dist; }
2.5 Maximum flow (Ford–Fulkerson)
#pragma once
#include #include #include
template struct Edge { int from, to; D cap; D flow;
Edge(int from, int to, D cap) : from(from), to(to), cap(cap), flow(0) {} Edge() = default;
int other(int origin) const { return origin == from? to : from; }
D max_increase(int origin) const { return origin == from? cap - flow : flow; }
void increase(int origin, D inc) { flow += (origin == from? inc : -inc); }
bool saturated(int origin) const { return max_increase(origin) == 0; } };
template std::vector<std::vector> adj_list(int N, std::vector<Edge> const &edges) { std::vector<std::vector> adj(N); for (int i = 0; i < (int) edges.size(); ++i) { Edge const &e = edges[i]; adj[e.from].push_back(i); adj[e.to].push_back(i); }
return adj; }
template D max_flow_body(int N, std::vector<Edge> &edges, std::vector<std::vector> const &adj, int source, int sink) { for (Edge &e : edges) e.flow = 0;
auto augment = [&] () -> D { std::vector pred(N, -1);
std::queue Q; Q.push(source);
bool found = false; while (!found && !Q.empty()) { int node = Q.front(); Q.pop();
for (int i : adj[node]) { Edge const &e = edges[i]; int other = e.other(node);
if (pred[other] == -1 && !e.saturated(node)) { Q.push(other); pred[other] = i; if (other == sink) { found = true; break; } } } }
if (found) { D max_inc = std::numeric_limits::max(); int node = sink; while (node != source) { Edge const &e = edges[pred[node]]; max_inc = std::min(max_inc, e.max_increase(e.other(node))); node = e.other(node); }
node = sink; while (node != source) { Edge &e = edges[pred[node]]; e.increase(e.other(node), max_inc); node = e.other(node); }
return max_inc; }
return 0; };
D max_flow = 0; D inc; do { inc = augment(); max_flow += inc; } while (inc > 0);
return max_flow; }
template D max_flow(int N, std::vector<Edge> &edges, int source, int sink) _/* Given a directed graph with edge capacities, computes the maximum flow
- from source to sink.
- Returns the maximum flow. Mutates input: assigns flows to the edges. */_ { std::vector<std::vector> adj = adj_list(N, edges); return max_flow_body(N, edges, adj, source, sink); }
void set_parent(int child_root, int parent_root) { parent[child_root] = parent_root; setsize[parent_root] += setsize[child_root]; } public: DisjointSets(int size) : size(size), setsize(size, 1), n_disjoint(size), rank(size, 0) { parent.assign(size, 0); std::iota(parent.begin(), parent.end(), 0); }
bool same_set(int i, int j) { return find(i) == find(j); }
int find(int i) { if (parent[i] != i) parent[i] = find(parent[i]); return parent[i]; }
void unite(int i, int j) { i = find(i); j = find(j); if (i != j) { --n_disjoint; if (rank[i] > rank[j]) set_parent(j, i); else if (rank[j] > rank[i]) set_parent(i, j); else { set_parent(j, i); ++rank[i]; } } } };
3.2 Fenwick Tree
#pragma once
#include
template<typename T, typename F> struct FenwickTree _/* Given an associative binary operation op (e.g. +) with identity element id
- (e.g. 0), the Fenwick tree represents an array x[0], ..., x[N - 1] and
- allows "adding" a value to any x[i], as well as computing the prefix "sum"
- x[0] + ... + x[i - 1], both in time O(n log n).
- While the implementation uses 1-indexing for bit-twiddling purposes,
- the API is 0-indexed._
std::vector data; F op; // F is e.g. T (*)(T, T) or std::function<T (T, T)> T id;
FenwickTree(int N, F op, T id) : data(N, 0), op(op), id(id) {}
// Internal one-indexing: // Parent of i: i - LSB(i) // Successor of parent of i: i + LSB(i) // where LSB(i) = i & -i T &a(int i) { return data[i - 1]; }
// Sum of the first i elements, from 0 to i - 1 T sum(int i) { // --i; ++i; (1-indexing cancels) T acc = id; while (i) { acc = op(acc, a(i)); i -= i & -i; }
return acc; }
void add(int i, T val) { ++i; // 1-indexing while (i <= (int) data.size()) { a(i) = op(a(i), val); i += i & -i; } } };
// Specialization for prefix sums: op = +, id = 0 template struct StandardFenwickTree : FenwickTree<T, T ()(T, T)> { StandardFenwickTree(int N) : FenwickTree<T, T ()(T, T)>( N, [] (T a, T b) { return a + b; }, T(0)) {} };
4 String algorithms
4.1 KMP
#include #include
class KMP { public: int const m; // Pattern length
std::string const P; // Pattern
// Prefix function // pf[q] = max {k | k < q, P[0..=k] suffix of P[0..=q]} // pf[q] {-1, 0, 1, 2, ...} std::vector pf;
KMP(std::string const &pattern) : m(pattern.length()), P(pattern), pf(m) { // Compute prefix function pf[0] = -1; for (int q = 0; q < m - 1; ++q) { int k = pf[q];
while (k != -1 && P[k + 1] != P[q + 1]) k = pf[k];
pf[q + 1] = P[k + 1] == P[q + 1]? k + 1 : -1; } }
std::vector match(std::string const &T) // Returns the index in T of all matches of P. { int const n = T.length(); std::vector matches;
int q = -1; // index in P of last matched letter for (int i = 0; i < n; ++i) { while (q != -1 && P[q + 1] != T[i]) q = pf[q];
if (P[q + 1] == T[i]) ++q;
if (q == m - 1) { matches.push_back(i - m + 1); q = pf[q]; } }
return matches; } };
5 Number theory
5.1 Combinatorics
#pragma once
template T int_pow(T x, int n) { T ans = 1;
while (n > 0) if (n % 2 == 0) x *= x, n /= 2; else ans *= x, --n;
return ans; }
template T fact(int n) { T ans = 1; for (int k = 1; k <= n; ++k) ans *= k;
return ans; }
template T nPr(int n, int k) { T ans = 1; for (int i = 0; i < k; ++i) ans *= n - i;
return ans; }
template T nCr(int n, int k) { if (k > n / 2) k = n - k;
T ans = 1; for (int i = 0; i < k; ++i) ans = (ans * (n - i)) / (i + 1);
return ans; }
5.2 Extended Euclidean algorithm
#pragma once
#include #include
template T sign(T x) { return (x > T(0)) - (x < T(0)); }
template T mod(T x, T m) { return (x % m + m) % m; }
template std::pair<T, T> extended_euclidean(T a, T b, T gcd = nullptr) _/ Solve a x + b y = gcd(a, b)._
A[k + 1] = A[k]; for (int c = w[k]; c <= capacity; ++c) A[k + 1][c] = std::max(A[k][c], A[k][c - w[k]] + v[k]); }
// Find best capacity int best_c = 0; for (int c = 0; c <= capacity; ++c) if (A[N][c] > A[N][best_c]) best_c = c;
// Backtrack int c = best_c; for (int k = N - 1; k >= 0; --k) // Either keep or remove element k if (w[k] <= c && A[k + 1][c] == A[k][c - w[k]] + v[k]) { ans.push_back(k); c -= w[k]; }
return A[N][best_c]; }
6.2 Longest increasing subsequence
#pragma once
#include "binary_search.hpp"
#include #include
template int longest_increasing_subsequence(std::vector const &x, std::vector &ans) _/* Given a sequence x, constructs the longest strictly increasing subsequence;
- i.e. an index sequence ans such that
- x[ans[0]] < x[ans[1]] < ...
- Returns the length of this sequence. */_ { int const N = x.size(); std::vector pred(N, -1);
// At each time i: A[j] = index k (0 <= k < i) of smallest endpoint x[k] of // an increasing subsequence of length j + 1. std::vector A;
for (int i = 0; i < N; ++i) { int const J = A.size(); // Length of longest subsequence so far
// Binary search for the location j to insert x: // x[A[j]] < x[i] <= x[A[j + 1]] auto p = [&] (int j) { return x[A[j]] < x[i]; };
int j = int_binsearch_last(p, 0, J);
if (j == J - 1) // x[i] is the endpoint of a new sequence of length J + 1 A.push_back(i); else // x[i] < x[A[j + 1]] // x[i] is a smaller endpoint of a sequence of length j A[j + 1] = i;
if (j != -1) pred[i] = A[j]; }
// Backtrack ans.assign(A.size(), -1); int i = A.empty()? -1 : A[A.size() - 1]; for (auto it = ans.rbegin(); it != ans.rend(); ++it) { *it = i; i = pred[i]; }
return A.size(); }
6.3 Interval cover
#pragma once
#include "binary_search.hpp"
#include #include #include #include
template int interval_cover(std::vector<std::pair<T, T>> const &intervals, T lo, T hi, std::vector &ans) _/*
- Given a vector of closed intervals [a_i, b_i], select a smallest subset
- that covers the target interval [lo, hi].
- Returns the minimal number of such intervals, and populates ans with their
- indices in intervals. */_ { T const NEG_INF = std::numeric_limits::lowest(); int const N = intervals.size();
// Sort intervals by starting point // retaining information about original order std::vector idx(N); std::iota(idx.begin(), idx.end(), 0); std::sort(idx.begin(), idx.end(), [&] (int i, int j) {
return intervals[i].first < intervals[j].first; });
// To save typing: I(i) = intervals[idx[i]] auto I = [&] (int i) -> std::pair<T, T> const & { return intervals[idx[i]]; };
// Invariant: (lo, hi] is the interval that is not yet covered, // except for before the first interval is added - then [lo, hi] is not yet // covered. T last_lo = NEG_INF; do { // Find the intervals with starting points in [last_lo, lo] with a // binary search, then find the one with the largest endpoint // by linear search T best_endpoint = lo; int best_i = -1;
int i = int_binsearch_first( [&] (int i) { return I(i).first >= last_lo; }, 0, N); for (; i < N && I(i).first <= lo; ++i) if ((ans.empty() && I(i).second >= best_endpoint) || I(i).second > best_endpoint) { best_endpoint = I(i).second; best_i = i; }
if (best_i == -1) { ans.clear(); return -1; // Impossible to cover }
ans.push_back(idx[best_i]); last_lo = lo; lo = best_endpoint; } while (lo < hi);
return ans.size(); }
7 Miscellaneous algorithms
7.1 Binary search
#pragma once
#include
template int int_binsearch_last(P p, int lo, int hi) _/* Takes a predicate p: int -> bool.
- Assumes that there is some I (lo <= I < hi) such that p(i) = (i <= I),_
_* i.e. p starts out true and then switches to false after I.
- Finds and returns I (the last number i such that p(i) is true),
- or lo - 1 if no such I exists (including if lo = hi). */_ { while (hi - lo > 1) { int mid = lo + (hi - lo) / 2; if (p(mid)) lo = mid; // mid <= I else hi = mid; // I < mid }
return (lo != hi && p(lo))? lo : lo - 1; }
template int int_binsearch_first(P p, int lo, int hi) _/* Takes a predicate p: int -> bool.
- Assumes that there is some I (lo <= I < hi) such that p(i) = (i >= I).
- i.e. p starts out false and then switches to true at I.
- Finds and returns I (the first number i such that p(i) is true),
- or hi if no such I exists (including if lo = hi). */_ { while (hi - lo > 1) { int mid = lo + (hi - lo - 1) / 2; if (p(mid)) hi = mid + 1; // mid >= I, search [lo, mid + 1) else lo = mid + 1; // I > mid, search [mid + 1, hi) }
return (lo != hi && p(lo))? lo : hi; }
template double double_binsearch_last(P p, double lo, double hi, int n_it) _/* Takes a predicate P: double -> bool.
- Assumes that there is some X (lo <= X <= hi) such that p(x) = (x <= X),
- i.e. p starts out true and then switches to false after X.
- Finds and returns X (the last number x such that p(x) is true),
- or -infinity if no such X exists.
- This version runs a fixed number of iterations; to get results with a given
- precision, use
while (hi - lo > eps) instead (but this is probably a bad - idea because of rounding errors). */_ { while (n_it--) { double mid = (lo + hi) / 2; if (p(mid)) lo = mid; // mid <= X else hi = mid; // X < mid }
return p(lo)? lo : -std::numeric_limits::infinity(); }
template double double_binsearch_first(P p, double lo, double hi, int n_it)
// vector x; // GaussianSolver solver(A, b, x); // solver.solve() // Solution to Ax = b is now in x. // A and b are modified into the reduced row echelon form of the system. // Other variables: solver.consistent, solver.rank, solver.determinant
template class GaussianSolver { private: // Add row i to row j void add(int i, int j, T factor) { b[j] += factor * b[i]; for (int k = 0; k < M; ++k) A[j][k] += factor * A[i][k]; }
void swap_rows(int i, int j) { std::swap(b[i], b[j]); for (int k = 0; k < M; ++k) std::swap(A[i][k], A[j][k]);
determinant = -determinant; }
void scale(int i, T factor) { b[i] *= factor; for (int k = 0; k < M; ++k) A[i][k] *= factor;
determinant /= factor; }
public: std::vector<std::vector> &A; std::vector &b; std::vector &x; int const N, M; // sqrt(std::numeric_limits::epsilon()) by default, epsilon is too small! T eps; int rank; T determinant; bool consistent; std::vector uniquely_determined;
GaussianSolver(std::vector<std::vector> &A, std::vector &b, std::vector &x, T eps = sqrt(std::numeric_limits::epsilon())) : A(A), b(b), x(x), N(A.size()), M(A[0].size()), eps(eps), rank(0), determinant(1), consistent(true) { uniquely_determined.assign(M, false); }
void solve() { // pr, pc: pivot row and column for (int pr = 0, pc = 0; pr < N && pc < M; ++pr, ++pc) { // Find pivot with largest absolute value int best_r = -1; T best = 0; for (int r = pr; r < N; ++r) if (std::abs(A[r][pc]) > best) { best_r = r; best = std::abs(A[r][pc]); }
if (std::abs(best) <= eps) { // No pivot found --pr; // only increase pc in the next iteration continue; }
// Rank = number of pivots ++rank;
// Move pivot to top and scale to 1 swap_rows(pr, best_r); scale(pr, (T) 1 / A[pr][pc]); A[pr][pc] = 1; // for numerical stability
// Eliminate entries below pivot for (int r = pr + 1; r < N; ++r) { add(pr, r, -A[r][pc]); A[r][pc] = 0; // for numerical stability } }
// Eliminate entries above pivots for (int pr = rank - 1; pr >= 0; --pr) { // Find pivot int pc = 0; while (std::abs(A[pr][pc]) <= eps) ++pc;
for (int r = pr - 1; r >= 0; --r) { add(pr, r, -A[r][pc]); A[r][pc] = 0; // for numerical stability } }
// Check for inconsistency: an equation of the form 0 = 1 for (int r = N - 1; r >= rank; --r) if (std::abs(b[r]) > eps) { consistent = false; return; }
// Calculate a solution for x // One solution is setting all non-pivot variables to 0 x.assign(M, 0); for (int pr = 0; pr < rank; ++pr)
for (int pc = 0; pc < M; ++pc) if (std::abs(A[pr][pc]) > eps) { // Pivot; A[pr][pc] == 1 x[pc] = b[pr]; break; }
// Mark variables as uniquely determined or not for (int pr = 0; pr < rank; ++pr) { int nonzero_count = 0; int pc = -1; for (int c = 0; c < M; ++c) if (std::abs(A[pr][c]) > eps) { if (nonzero_count == 0) pc = c; // Pivot ++nonzero_count; } if (nonzero_count == 1) uniquely_determined[pc] = true; } } };
7.4 Simplex algorithm (linear programming)
// Simplex algorithm for linear programming. // Written using the theory from // Cormen, Leiserson, Rivest, Stein: Introduction to Algorithms
#pragma once
#include #include #include #include #include #include
// For debug() #include #include
template class SimplexSolver { public: // Standard form: Maximize c x subject to A x <= b; x >= 0. // Slack form: Basic variables xb and nonbasic variables xn. x = [xb, xn] // Maximize z = nu + c xn subject to xb = b - A xn; x >= 0.
T const eps;
int n; // number of nonbasic variables int m; // number of constraints std::vector<std::vector> A; std::vector b; std::vector c;
T nu;
// Holds INDICES of all variables (ranging from 0 to n + m - 1) std::vector nonbasic_vars, basic_vars;
bool feasible;
SimplexSolver(std::vector<std::vector> A, std::vector b, std::vector c, T eps = sqrt(std::numeric_limits::epsilon())) : eps(eps), n(c.size()), m(b.size()), A(A), b(b), c(c), nu(0), nonbasic_vars(n), basic_vars(m), feasible(true) { assert((int) A.size() == m && (int) A[0].size() == n);
// Transform from standard form to slack form // Initially: nonbasic_vars: 0 to n - 1, basic_vars: n to n + m - 1 std::iota(nonbasic_vars.begin(), nonbasic_vars.end(), 0); std::iota(basic_vars.begin(), basic_vars.end(), n); }
// xn_e: "entering" variable: nonbasic -> basic // xb_l: "leaving" variable: basic -> nonbasic void pivot(int e, int l) { std::swap(nonbasic_vars[e], basic_vars[l]); int const e_new = l, l_new = e; // Just to avoid confusion
std::vector<std::vector> A_new(A); std::vector b_new(b); std::vector c_new(c); T nu_new;
// New constraint for xn_e: replace // xb_l = b_l - A_lj xn_j // with // (*) xn_e = b_l / A_le - xb_l / A_le - A_lj xn_j / A_le (j ≠ e) b_new[e_new] = b[l] / A[l][e]; A_new[e_new][l_new] = (T) 1 / A[l][e]; for (int j = 0; j < n; ++j) if (j != e) A_new[e_new][j] = A[l][j] / A[l][e];
// Substitute (*) in the other constraint equations: // In each xb_i = b_i - A_ij xn_j (i ≠ l), replace A_ie xn_e // with A_ie (b_l / A_le - xb_l / A_le - A_lj xn_j / A_le (j ≠ e)) for (int i = 0; i < m; ++i) if (i != l) { b_new[i] -= A[i][e] / A[l][e] * b[l]; A_new[i][l_new] = -A[i][e] / A[l][e]; for (int j = 0; j < n; ++j) if (j != e) A_new[i][j] -= A[i][e] * A[l][j] / A[l][e]; }
// Find the index of x0, now a non-basic variable int j = std::find(nonbasic_vars.begin(), nonbasic_vars.end(), n + m - 1) - nonbasic_vars.begin();
// Erase x0 from the constraints and the list of variables for (std::vector &row : A) row.erase(row.begin() + j); nonbasic_vars.erase(nonbasic_vars.begin() + j); --n;
// Restore original objective function, substituting basic variables // with their RHS nu = original_nu; c.assign(n, 0); // Loop over all originally non-basic variables for (int var = 0; var < n; ++var) { int j = std::find(nonbasic_vars.begin(), nonbasic_vars.end(), var) - nonbasic_vars.begin(); if (j != n) c[j] += original_c[var]; else { int i = std::find(basic_vars.begin(), basic_vars.end(), var) - basic_vars.begin(); // Substitute xb_i = b_i - A_ij xn_j nu += original_c[var] * b[i]; for (int j = 0; j < n; ++j) c[j] -= original_c[var] * A[i][j]; } } } else { --n; feasible = false; } }
void debug() const { printf("Nonbasic vars: "); for (int j : nonbasic_vars) printf("x[%d] ", j); std::cout << std::endl; printf("Basic vars: "); for (int i : basic_vars) printf("x[%d] ", i); std::cout << std::endl;
std::cout << "Optimize " << nu; for (int j = 0; j < n; ++j) { std::cout << " + (" << c[j] << ") * "; printf("x[%d]", nonbasic_vars[j]); } puts(""); for (int i = 0; i < m; ++i) { printf("x[%d] = ", basic_vars[i]); std::cout << "(" << b[i] << ")"; for (int j = 0; j < n; ++j) { std::cout << " - (" << A[i][j] << ") * "; printf("x[%d]", nonbasic_vars[j]);
puts(""); }
puts(""); } };
7.5 Basis conversion
#pragma once
#include #include #include
std::string const DIGITS = "0123456789ABCDEF";
unsigned long long int basis_string_to_number(std::string &s, int b) { unsigned long long int result = 0ULL; for (char d : s) { result = b * result
- (std::find(DIGITS.begin(), DIGITS.end(), d) - DIGITS.begin()); }
return result; }
std::string number_to_basis_string(unsigned long long int n, int b) { std::vector ds; do { ds.push_back(DIGITS[n % b]); n = n / b; } while (n != 0);
return std::string(ds.rbegin(), ds.rend()); }
8 Mathematical objects
8.1 Fraction
#pragma once
#include #include #include #include #include
template struct Fraction { T n, d;
Fraction(T n, T d) : n(n), d(d) { reduce(); };
Fraction() : Fraction(0) {} Fraction(T n) : Fraction(n, 1) {} Fraction(Fraction const &) = default; Fraction &operator=(Fraction const &) = default;
void reduce() { T gcd = std::__gcd(n, d); n /= gcd; d /= gcd;
if (d < 0) { n = -n; d = -d; } }
bool operator==(Fraction const &other) const { return n * other.d == other.n * d; }
bool operator<(Fraction const &other) const { assert(d > 0 && other.d > 0); return n * other.d < other.n * d; }
bool operator>(Fraction const &other) const { assert(d > 0 && other.d > 0); return n * other.d > other.n * d; }
bool operator<=(Fraction const &other) const { return !(*this > other); }
bool operator>=(Fraction const &other) const { return !(*this < other); }
Fraction operator+(Fraction const &other) const { return Fraction(n * other.d + other.n * d, d * other.d); }
Fraction operator-(Fraction const &other) const { return Fraction(n * other.d - other.n * d, d * other.d); }
Fraction operator*(Fraction const &other) const { return Fraction(n * other.n, d * other.d); }
Fraction operator/(Fraction const &other) const {
return Fraction(n * other.d, d * other.n); }
Fraction &operator+=(Fraction const &other) { return *this = *this + other; }
Fraction &operator-=(Fraction const &other) { return *this = *this - other; }
Fraction &operator*=(Fraction const &other) { return *this = *this * other; }
Fraction &operator/=(Fraction const &other) { return *this = *this / other; }
Fraction operator+() const { return *this; }
Fraction operator-() const { return Fraction(-n, d); }
void print(FILE *f) const { fprintf(f, "%lld / %lld", (long long int) n, (long long int) d); }
friend std::ostream &operator<<(std::ostream &s, Fraction const &frac) { return s << frac.n << " / " << frac.d; } };
namespace std { template class numeric_limits<Fraction> { public: static constexpr Fraction max() { return Fraction(numeric_limits::max()); } static constexpr Fraction lowest() { return Fraction(numeric_limits::lowest()); } static constexpr Fraction epsilon() { return Fraction(0); } };
template Fraction abs(Fraction const &f) { return Fraction(abs(f.n), abs(f.d)); }
bool operator>=(Bignum const &other) const { return !(*this < other); }
Bignum operator-() const { Bignum ans = *this; ans.sign = -ans.sign;
return ans.trim(); }
Bignum abs() const { Bignum ans = *this; ans.sign = 1; return ans; }
Bignum &operator+=(Bignum const &other) { check_basis(other);
if (sign != other.sign) return *this -= -other;
digits.resize(1 + std::max(n_digits(), other.n_digits())); int carry = 0; for (int i = 0; i < n_digits(); ++i) { int next_digit = get_digit(i) + other.get_digit(i) + carry;
carry = next_digit / basis; next_digit %= basis;
digits[i] = next_digit; } assert(carry == 0);
return trim(); }
Bignum &operator-=(Bignum const &other) { check_basis(other);
if (sign != other.sign) return *this += -other;
if (abs() < other.abs()) return *this = -(other - *this);
digits.resize(1 + std::max(n_digits(), other.n_digits())); bool borrow = false; for (int i = 0; i < n_digits() - 1; ++i) { int next_digit = get_digit(i) - other.get_digit(i);
borrow = next_digit < 0; if (borrow) { --digits[i + 1]; next_digit += basis; }
digits[i] = next_digit; }
return trim(); }
Bignum operator+(Bignum const &other) const { check_basis(other);
Bignum ans = *this; return ans += other; }
Bignum operator-(Bignum const &other) const { check_basis(other);
Bignum ans = *this; return ans -= other; }
Bignum &operator+=(int n) { return *this += Bignum(n, basis); }
Bignum &operator-=(int n) { return *this -= Bignum(n, basis); }
Bignum &operator++() { return *this += Bignum(1, basis); }
Bignum &operator--() { return *this -= Bignum(1, basis); } Bignum operator+(int n) { return *this + Bignum(n, basis); }
Bignum operator-(int n) { return *this - Bignum(n, basis); }
friend Bignum operator+(int n, Bignum const &b) { return b + n; }
friend Bignum operator-(int n, Bignum const &b) { return Bignum(n, b.basis) - b; }
Bignum &operator*=(int n) { if (n < 0) { sign = -sign;
return *this *= -n; } if (n == 0) return *this = Bignum(0, basis); if (n == 1) return *this;
if (n % 2 == 0) { *this += *this; return *this *= n / 2; } else { Bignum tmp = this; return (this *= n - 1) += tmp; } }
Bignum operator*(int n) const { Bignum ans = *this; return ans *= n; }
friend Bignum operator*(int n, Bignum const &b) { return b * n; }
Bignum operator*(Bignum const &other) const { return (std::max(n_digits(), other.n_digits()) < KARATSUBA_CUTOFF) ? naive_mult(other) : karatsuba_mult(other); }
Bignum naive_mult(Bignum const &other) const { check_basis(other);
Bignum ans(0, basis); for (int i = 0; i < n_digits(); ++i) ans += (other << i) * digits[i];
return ans; }
Bignum karatsuba_mult(Bignum const &other) const { check_basis(other);
if (n_digits() == 1 || other.n_digits() == 1) return naive_mult(other);
if (sign == -1) return -abs().karatsuba_mult(other); if (other.sign == -1) return -karatsuba_mult(other.abs());
int k = std::max(n_digits(), other.n_digits()) / 2; Bignum a, b, c, d; std::tie(a, b) = split(k); std::tie(c, d) = other.split(k);
// *this = a << k + b
// other = c << k + d // *this * other = ac << 2k + (ad + bc) << k + bd // ad + bc = (a + b)(c + d) - ac - bd Bignum ac = a * c; Bignum bd = b * d; Bignum ad_plus_bc = (a + b) * (c + d) - ac - bd;
return (ac << (2 * k)) + (ad_plus_bc << k) + bd; }
Bignum &operator*=(Bignum const &other) { check_basis(other);
return *this = *this * other; }
Bignum &operator/=(int n) { assert(std::abs(n) < basis);
if (n < 0) { sign = -sign; return *this /= -n; } if (n == 0) throw std::runtime_error("Division by zero"); if (n == 1) return *this;
int carry = 0; for (auto it = digits.rbegin(); it != digits.rend(); ++it) { long long numerator = *it + ((long long) carry) * basis; *it = numerator / n; carry = numerator % n; }
return trim(); }
Bignum operator/(int n) const { Bignum ans = *this; return ans /= n; }
Bignum operator/(Bignum const &n) const { if (n < 0) return -(*this / -n); if (sign < 0) return -(abs() / n); if (n == 0) throw std::runtime_error("Division by zero"); if (n == 1) return *this;
// Binary search for last number x such that n * x <= *this // Total time: O(log(n) d log(d)) = O(d log²(d)) // where d is the number of digits in n Bignum lo(0, basis), hi(*this); while (hi - lo > 1) { Bignum mid = lo + (hi - lo) / 2; if (n * mid <= *this) lo = mid; else hi = mid;
