text:
#pragma once
#include <bits/stdc++.h>
using namespace std;

// Computes the maximum weighted matching in a general undirected weighted
//   graph, such that each vertex is incident with at most one edge in the
//   matching, and the selected edges have the maximum sum of weights
// Vertices are 0-indexed
// Template Arguments:
//   T: the type of each edge weight
// Constructor Arguments:
//   matrix: a V x V matrix containing the weight of the edge between two
//     vertices, 0 if an edge doesn't exist, must be symmetric and
//      non negative
//   INF: a value for infinity
// Fields:
//   V: the number of vertices in the graph
//   cardinality: the cardinality of the maximum matching
//   INF: a value for infinity
//   cost: the maximum cost of the weighted matching
//   mate: the other vertex in the matching, or -1 if unmatched
// In practice, has a moderate constant
// Time Complexity:
//   constructor: O(V^3)
// Memory Complexity: O(V^2)
// Tested:
//   https://uoj.ac/problem/81
//   https://judge.yosupo.jp/problem/general_weighted_matching
template <class T> struct GeneralWeightedMaxMatch {
  struct Edge {
    int v, w; T weight; Edge() : v(0), w(0), weight(T()) {}
  };
  int V, VX, cardinality, curStamp; T INF, cost;
  vector<int> mate, match, slack, st, par, S, stamp; vector<T> lab;
  vector<vector<int>> flo, floFrom; vector<vector<Edge>> G; queue<int> q;
  T eDelta(const Edge &e) {
    return lab[e.v] + lab[e.w] - G[e.v][e.w].weight * T(2);
  }
  void updateSlack(int v, int x) {
    if (!slack[x] || eDelta(G[v][x]) < eDelta(G[slack[x]][x])) slack[x] = v;
  }
  void setSlack(int x) {
    slack[x] = 0; for (int v = 1; v <= V; v++)
      if (G[v][x].weight > T() && st[v] != x && S[st[v]] == 0)
        updateSlack(v, x);
  }
  void qPush(int x) {
    if (x <= V) q.push(x);
    else for (int t : flo[x]) qPush(t);
  }
  void setSt(int x, int b) {
    st[x] = b; if (x > V) for (int t : flo[x]) setSt(t, b);
  }
  int getPr(int b, int xr) {
    int pr = find(flo[b].begin(), flo[b].end(), xr) - flo[b].begin();
    if (pr % 2) {
      reverse(flo[b].begin() + 1, flo[b].end());
      return int(flo[b].size()) - pr;
    }
    return pr;
  }
  void setMatch(int v, int w) {
    Edge &e = G[v][w]; match[v] = e.w; if (v <= V) return;
    int xr = floFrom[v][e.v], pr = getPr(v, xr);
    for (int i = 0; i < pr; i++) setMatch(flo[v][i], flo[v][i ^ 1]);
    setMatch(xr, w); rotate(flo[v].begin(), flo[v].begin() + pr, flo[v].end());
  }
  void augment(int v, int w) {
    while (true) {
      int xnw = st[match[v]]; setMatch(v, w); if (!xnw) return;
      setMatch(xnw, v = st[par[w = xnw]]);
    }
  }
  int lca(int v, int w) {
    for (curStamp++; v || w; swap(v, w)) {
      if (!v) continue;
      if (stamp[v] == curStamp) return v;
      stamp[v] = curStamp; v = st[match[v]]; if (v) v = st[par[v]];
    }
    return 0;
  }
  void addBlossom(int v, int anc, int w) {
    int b = V + 1; while (b <= VX && st[b]) b++;
    if (b > VX) VX++;
    lab[b] = T(); S[b] = 0; match[b] = match[anc]; flo[b] = vector<int>{anc};
    auto blossom = [&] (int x) {
      for (int y; x != anc; x = st[par[y]]) {
        flo[b].push_back(x); flo[b].push_back(y = st[match[x]]); qPush(y);
      }
    };
    blossom(v); reverse(flo[b].begin() + 1, flo[b].end()); blossom(w);
    setSt(b, b);
    for (int x = 1; x <= VX; x++) G[b][x].weight = G[x][b].weight = T();
    for (int x = 1; x <= V; x++) floFrom[b][x] = 0;
    for (int xs : flo[b]) {
      for (int x = 1; x <= VX; x++)
        if (G[b][x].weight == T(0) || eDelta(G[xs][x]) < eDelta(G[b][x])) {
          G[b][x] = G[xs][x]; G[x][b] = G[x][xs];
        }
      for (int x = 1; x <= V; x++) if (floFrom[xs][x]) floFrom[b][x] = xs;
    }
    setSlack(b);
  }
  void expandBlossom(int b) {
    for (int t : flo[b]) setSt(t, t);
    int xr = floFrom[b][G[b][par[b]].v], pr = getPr(b, xr);
    for (int i = 0; i < pr; i += 2) {
      int xs = flo[b][i], xns = flo[b][i + 1]; par[xs] = G[xns][xs].v;
      S[xs] = 1; S[xns] = slack[xs] = slack[xns] = 0; qPush(xns);
    }
    S[xr] = 1; par[xr] = par[b];
    for (int i = pr + 1; i < int(flo[b].size()); i++) {
      int xs = flo[b][i]; S[xs] = -1; setSlack(xs);
    }
    st[b] = 0;
  }
  bool onFoundEdge(const Edge &e) {
    int v = st[e.v], w = st[e.w]; if (S[w] == -1) {
      par[w] = e.v; S[w] = 1; slack[w] = 0;
      int nv = st[match[w]]; S[nv] = slack[nv] = 0; qPush(nv);
    } else if (S[w] == 0) {
      int anc = lca(v, w);
      if (!anc) { augment(v, w); augment(w, v); return true; }
      addBlossom(v, anc, w);
    }
    return false;
  }
  bool matching() {
    q = queue<int>(); for (int x = 1; x <= VX; x++) {
      S[x] = -1; slack[x] = 0;
      if (st[x] == x && !match[x]) { par[x] = S[x] = 0; qPush(x); }
    }
    if (q.empty()) return false;
    while (true) {
      while (!q.empty()) {
        int v = q.front(); q.pop(); if (S[st[v]] == 1) continue;
        for (int w = 1; w <= V; w++)
          if (G[v][w].weight > T() && st[v] != st[w]) {
            if (eDelta(G[v][w]) == T()) {
              if (onFoundEdge(G[v][w])) return true;
            } else updateSlack(v, st[w]);
          }
      }
      T d = INF; for (int b = V + 1; b <= VX; b++)
        if (st[b] == b && S[b] == 1) d = min(d, lab[b] / T(2));
      for (int x = 1; x <= VX; x++) if (st[x] == x && slack[x]) {
        if (S[x] == -1) d = min(d, eDelta(G[slack[x]][x]));
        else if (S[x] == 0) d = min(d, eDelta(G[slack[x]][x]) / T(2));
      }
      for (int v = 1; v <= V; v++) {
        if (S[st[v]] == 0) {
          if (lab[v] <= d) return false;
          lab[v] -= d;
        } else if (S[st[v]] == 1) lab[v] += d;
      }
      for (int b = V + 1; b <= VX; b++) if (st[b] == b && S[b] != -1) {
        if (S[b] == 0) lab[b] += d * T(2);
        else lab[b] -= d * T(2);
      }
      q = queue<int>(); for (int x = 1; x <= VX; x++)
        if (st[x] == x && slack[x] && st[slack[x]] != x
            && eDelta(G[slack[x]][x]) == T() && onFoundEdge(G[slack[x]][x]))
          return true;
      for (int b = V + 1; b <= VX; b++)
        if (st[b] == b && S[b] == 1 && lab[b] == T()) expandBlossom(b);
    }
    return false;
  }
  GeneralWeightedMaxMatch(const vector<vector<T>> &matrix,
                          T INF = numeric_limits<T>::max())
      : V(matrix.size()), VX(V), cardinality(0), curStamp(0), INF(INF),
        cost(T()), mate(V, -1), match(V * 2 + 1, 0), slack(V * 2 + 1, 0),
        st(V * 2 + 1, 0), par(V * 2 + 1, 0), S(V * 2 + 1, 0),
        stamp(V * 2 + 1, 0), lab(V * 2 + 1, T()), flo(V * 2 + 1),
        floFrom(V * 2 + 1, vector<int>(V + 1, 0)),
        G(V * 2 + 1, vector<Edge>(V * 2 + 1)) {
    iota(st.begin(), st.begin() + V + 1, 0); T mx = T();
    for (int v = 1; v <= V; v++) for (int w = 1; w <= V; w++) {
      G[v][w].v = v; G[v][w].w = w; floFrom[v][w] = (v == w ?  v : 0);
      assert(G[v][w].weight >= T());
      mx = max(mx, G[v][w].weight = matrix[v - 1][w - 1]);
    }
    fill(lab.begin() + 1, lab.begin() + V + 1, mx);
    while (matching()) cardinality++;
    for (int v = 1; v <= V; v++) if (match[v] && match[v] < v) {
      mate[mate[v - 1] = G[v][match[v]].w - 1] = v - 1;
      cost += G[v][match[v]].weight;
    }
  }
};


text:
#pragma once
#include <bits/stdc++.h>
using namespace std;

// Given N people of type A and M people of type B, and a list of their ranked
//   preferences for partners, the goal is to arrange min(N, M) pairs such that
//   if a person x of type A prefers a person y of type B more than their
//   current partner, then person y prefers their current partner more than x
// It is guaranteed that a solution always exists
// Constructor Arguments:
//   aPrefs: a matrix of size N x M with aPrefs[i][j] representing the jth
//     preferred choice for the partner of type B for person i of type A;
//     aPrefs[i] must be a permutation from 0 to M - 1
//   bPrefs: a matrix of size M x N with bPrefs[j][i] representing the ith
//     preferred choice for the partner of type A for person j of type B;
//     bPrefs[j] must be a permutation from 0 to N - 1
// Fields:
//   N: the number of people of type A
//   M: the number of people of type B
//   bForA: a vector representing the partner of type B for each person of
//     type A, or -1 if that person of type A is unmatched
//   aForB: a vector representing the partner of type A for each person of
//     type B, or -1 if that person of type B is unmatched
// In practice, has a moderate constant
// Time Complexity:
//   constructor O(NM)
// Memory Complexity: O(NM)
// Tested:
//   https://www.spoj.com/problems/STABLEMP/
//   https://open.kattis.com/problems/jealousyoungsters
struct StableMarriage {
  int N, M; vector<int> bForA, aForB;
  StableMarriage(const vector<vector<int>> &aPrefs,
                 const vector<vector<int>> &bPrefs)
      : N(aPrefs.size()), M(bPrefs.size()), bForA(N, -1), aForB(M, -1) {
    bool rev = N > M; if (rev) { swap(N, M); bForA.swap(aForB); }
    auto &A = rev ? bPrefs : aPrefs, &B = rev ? aPrefs : bPrefs;
    vector<vector<int>> bRankOfA(M, vector<int>(N));
    for (int b = 0; b < M; b++) for (int a = 0; a < N; a++)
      bRankOfA[b][B[b][a]] = a;
    queue<int> q; for (int a = 0; a < N; a++) q.push(a);
    vector<int> cur(N, 0); while (!q.empty()) {
      int a = q.front(); q.pop(); while (true) {
        int b = A[a][cur[a]++]; if (aForB[b] == -1) { aForB[b] = a; break; }
        else if (bRankOfA[b][a] < bRankOfA[b][aForB[b]]) {
          q.push(aForB[b]); aForB[b] = a; break;
        }
      }
    }
    for (int b = 0; b < M; b++) if (aForB[b] != -1) bForA[aForB[b]] = b;
    if (rev) { swap(N, M); bForA.swap(aForB); }
  }
};

repo_name: wesley-a-leung/Resources

mime_type: text/x-c++


text:
#pragma once
#include <bits/stdc++.h>
using namespace std;

// Solves the assignment problem of matching N workers to M jobs with the
//   minimum cost where each worker can only be assigned to at most 1 job and
//   each job can only be assigned to at most 1 worker and there are exactly
//   min(N, M) assignments
// Maximum cost assignment can be found by negating all non infinity values in
//   the matrix and taking the negative of the minimum cost
// Maximum cost matching can be found by negating all non infinity values in
//   the matrix, setting all infinity values to 0 and taking the negative of
//   the minimum cost
// Template Arguments:
//   T: the type of the cost
// Constructor Arguments:
//   A: a matrix of size N x M with A[i][j] representing the cost of assigning
//     worker i to job j
//   INF: a value for infinity
// Fields:
//   N: the number of workers
//   M: the number of jobs
//   cost: the minimum cost for a valid assignment
//   jobForWorker: a vector representing the job assigned to each worker, or -1
//     if that worker is not assigned any job
//   workerForJob: a vector representing the worker assigned to each job, or -1
//     if that job is not assigned any worker
// In practice, has a very small constant
// Time Complexity:
//   constructor: O(N^2 M)
// Memory Complexity: O(NM)
// Tested:
//   https://judge.yosupo.jp/problem/assignment
//   https://open.kattis.com/problems/cordonbleu
//   https://www.spoj.com/problems/BABY/
//   https://dmoj.ca/problem/tle17c7p5
//   https://open.kattis.com/problems/workers
template <class T> struct HungarianAlgorithm {
  int N, M; T cost; vector<int> jobForWorker, workerForJob;
  HungarianAlgorithm(vector<vector<T>> A, T INF = numeric_limits<T>::max())
      : N(A.size()), M(A.empty() ? 0 : A[0].size()), cost(T()) {
    bool rev = N > M; if (rev) {
      swap(N, M); vector<vector<T>> B(N, vector<T>(M));
      for (int i = 0; i < N; i++) for (int j = 0; j < M; j++)
        B[i][j] = A[j][i];
      A.swap(B);
    }
    jobForWorker.assign(N, -1); workerForJob.assign(M + 1, N);
    auto add = [&] (pair<int, T> &a, const pair<int, T> &b) {
      a.first += b.first; a.second += b.second;
    };
    auto sub = [&] (pair<int, T> &a, const pair<int, T> &b) {
      a.first -= b.first; a.second -= b.second;
    };
    vector<pair<int, T>> d1(N + 1, make_pair(0, T()));
    vector<pair<int, T>> d2(M + 1, make_pair(0, T()));
    for (int i = 0; i < N; i++) {
      int j0 = M; workerForJob[j0] = i;
      vector<pair<int, T>> dist(M + 1, make_pair(1, T()));
      vector<int> par(M + 1, -1); vector<bool> done(M + 1, false); do {
        done[j0] = true; int i0 = workerForJob[j0], j1 = M;
        pair<int, T> delta = make_pair(1, T());
        for (int j = 0; j < M; j++) if (!done[j]) {
          pair<int, T> d = A[i0][j] >= INF ? make_pair(1, T())
                                           : make_pair(0, A[i0][j]);
          sub(d, d1[i0]); sub(d, d2[j]);
          if (dist[j].first > 0 || d < dist[j]) { dist[j] = d; par[j] = j0; }
          if (dist[j] < delta) delta = dist[j1 = j];
        }
        j0 = j1; for (int j = 0; j <= M; j++) {
          if (done[j]) { add(d1[workerForJob[j]], delta); sub(d2[j], delta); }
          else sub(dist[j], delta);
        }
      } while (workerForJob[j0] != N);
      for (; j0 != M; j0 = par[j0]) workerForJob[j0] = workerForJob[par[j0]];
    }
    workerForJob.pop_back(); for (int j = 0; j < M; j++) {
      if (workerForJob[j] == N) workerForJob[j] = -1;
      else jobForWorker[workerForJob[j]] = j;
    }
    cost = d2[M].first < 0 ? INF : -d2[M].second;
    if (rev) { swap(N, M); jobForWorker.swap(workerForJob); }
  }
};


text:
#pragma once
#include <bits/stdc++.h>
using namespace std;

// Computes the maximum matching in a general undirected graph, such that each
//   vertex is incident with at most one edge in the matching
// Vertices are 0-indexed
// Constructor Arguments:
//   G: a generic undirected graph structure
//     Required Functions:
//       operator [v] const: iterates over the adjacency list of vertex v
//         (which is a list of ints)
//       size() const: returns the number of vertices in the graph
// Fields:
//   V: the number of vertices in the graph
//   cardinality: the cardinality of the maximum matching
//   mate: the other vertex in the matching, or -1 if unmatched
// In practice, has a very small constant
// Time Complexity:
//   constructor: O(VE log V)
// Memory Complexity: O(V)
// Tested:
//   https://judge.yosupo.jp/problem/general_matching
//   https://codeforces.com/contest/1089/problem/B
//   https://uoj.ac/problem/79
struct GeneralUnweightedMaxMatching {
  int V, cardinality, front, back; vector<int> mate, first, q;
  vector<pair<int, int>> label;
  void rematch(int v, int w) {
    int t = mate[v]; mate[v] = w; if (mate[t] != v) return;
    if (label[v].second == -1) rematch(mate[t] = label[v].first, t);
    else { int x, y; tie(x, y) = label[v]; rematch(x, y); rematch(y, x); }
  }
  int find(int v) {
    return label[first[v]].first < 0 ? first[v] : first[v] = find(first[v]);
  }
  void relabel(int x, int y) {
    int r = find(x), s = find(y), join = 0; if (r == s) return;
    pair<int, int> h = label[r] = label[s] = make_pair(~x, y); while (true) {
      if (s != V) swap(r, s);
      r = find(label[mate[r]].first); if (label[r] == h) { join = r; break; }
      else label[r] = h;
    }
    for (int v : {first[x], first[y]})
      for (; v != join; v = first[label[mate[v]].first]) {
        label[v] = make_pair(x, y); first[q[back++] = v] = join;
      }
  }
  template <class Graph> bool augment(const Graph &G, int s) {
    front = back = 0; label[q[back++] = s] = make_pair(first[s] = V, -1);
    while (front < back) {
      int v = q[front++]; for (int w : G[v]) {
        if (mate[w] == V && w != s) { rematch(mate[w] = v, w); return true; }
        else if (label[w].first >= 0) relabel(v, w);
        else if (label[mate[w]].first == -1)
          label[mate[first[q[back++] = mate[w]] = w]].first = v;
      }
    }
    return false;
  }
  template <class Graph> GeneralUnweightedMaxMatching(const Graph &G)
      : V(G.size()), cardinality(0), front(0), back(0), mate(V + 1, V),
        first(V + 1, V), q(V), label(V + 1, make_pair(-1, -1)) {
    for (int v = 0; v < V; v++) if (mate[v] == V && augment(G, v)) {
      cardinality++; for (int i = 0; i < back; i++)
        label[q[i]] = label[mate[q[i]]] = make_pair(-1, -1);
      label[V] = make_pair(-1, -1);
    }
    mate.pop_back(); for (auto &&m : mate) if (m == V) m = -1;
  }
};


text:
#pragma once
#include <bits/stdc++.h>
using namespace std;

// Given N people, and a list of their ranked preferences for roommates,
//   the goal is to arrange N / 2 pairs such that if a person x prefers
//   a person y more than their current roommate, then person y prefers
//   their current roommate more than x
// Constructor Arguments:
//   prefs: a matrix of size N x (N - 1) with prefs[i][j] representing the jth
//     preferred choice for the person i
//     prefs[i] must be a permutation of 0 to N - 1, excluding i
// Fields:
//   N: the number of people
//   mate: a vector representing the roommate of each person; all -1 if there
//     is no stable matching
// In practice, has a moderate constant
// Time Complexity:
//   constructor O(N^2)
// Memory Complexity: O(N^2)
// Tested:
//   https://codeforces.com/contest/1423/problem/A
struct StableRoommates {
  struct NoMatch {}; int N; vector<int> mate;
  StableRoommates(vector<vector<int>> prefs) : N(prefs.size()), mate(N, -1) {
    if (N % 2 == 1 || N <= 0) return;
    vector<vector<int>> rnk(N, vector<int>(N, 0));
    vector<int> fr(N, 0), bk(N, N - 1), proposed(N, -1);
    vector<vector<bool>> active(N, vector<bool>(N, true)); queue<int> q;
    auto rem = [&] (int i, int j) { active[i][j] = active[j][i] = false; };
    auto clip = [&] (int i) {
      while (fr[i] < bk[i] && !active[i][prefs[i][fr[i]]]) fr[i]++;
      while (fr[i] < bk[i] && !active[i][prefs[i][bk[i] - 1]]) bk[i]--;
      if (fr[i] >= bk[i]) throw NoMatch();
    };
    auto add = [&] (int i, int j) {
      proposed[mate[i] = j] = i; while (true) {
        clip(j); if (prefs[j][bk[j] - 1] != i) rem(j, prefs[j][bk[j] - 1]);
        else break;
      }
    };
    auto nxt = [&] (int i) {
      clip(i); int j = prefs[i][fr[i]++];
      clip(i); prefs[i][--fr[i]] = j; return proposed[prefs[i][fr[i] + 1]];
    };
    for (int i = 0; i < N; i++) {
      q.push(i); for (int j = 0; j < N - 1; j++) rnk[i][prefs[i][j]] = j;
    }
    try {
      while (!q.empty()) {
        int i = q.front(); q.pop(); while (true) {
          clip(i); int j = prefs[i][fr[i]], i2 = proposed[j];
          if (i2 != -1 && rnk[j][i2] < rnk[j][i]) { rem(i, j); continue; }
          if (i2 != -1) { mate[i2] = proposed[j] = -1; q.push(i2); }
          add(i, j); break;
        }
      }
      int cur = 0; while (true) {
        for (; cur < N; cur++) {
          clip(cur); if (bk[cur] - fr[cur] > 1) break;
        }
        if (cur == N) break;
        vector<int> cyc1, cyc2; int i = cur, j = i;
        do { i = nxt(i); j = nxt(nxt(j)); } while (i != j);
        do { cyc1.push_back(j); j = nxt(j); } while (i != j);
        for (int k : cyc1) {
          j = mate[k]; cyc2.push_back(j); mate[k] = proposed[j] = -1;
          rem(k, j);
        }
        for (int k = 0; k < int(cyc1.size()); k++)
          add(cyc1[k], cyc2[(k + 1) % cyc2.size()]);
      }
    } catch (NoMatch &) { fill(mate.begin(), mate.end(), -1); }
  }
};


text:
#pragma once
#include <bits/stdc++.h>
using namespace std;

// Computes the maximum matching (and minimum vertex cover) on
//   an unweighted bipartite graph
// The maximum independent set is any vertex not in the minimum vertex cover
// Vertices are 0-indexed
// Constructor Arguments:
//   G: a generic undirected bipartite graph structure
//     Required Functions:
//       operator [v] const: iterates over the adjacency list of vertex v
//         (which is a list of ints)
//       size() const: returns the number of vertices in the graph
//   color: a vector of booleans of size V indicating the color of each vertex
// Fields:
//   V: the number of vertices in the graph
//   cardinality: the cardinality of the maximum matching
//   mate: the other vertex in the matching, or -1 if unmatched
//   inCover: a vector of booleans indicating whether a vertex is in the
//    minimum vertex cover or not
// In practice, has a very small constant
// Time Complexity:
//   constructor: O((V + E) sqrt V)
// Memory Complexity: O(V)
// Tested:
//   https://www.spoj.com/problems/MATCHING/
//   https://judge.yosupo.jp/problem/bipartitematching
//   https://dmoj.ca/problem/coci19c6p4
struct HopcroftKarpMaxMatch {
  int V, cardinality;
  vector<int> mate, lvl, q, type0; vector<bool> color, inCover, vis;
  template <class BipartiteGraph> bool bfs(const BipartiteGraph &G) {
    int front = 0, back = 0; for (int v : type0) {
      if (mate[v] == -1) lvl[q[back++] = v] = 0;
      else lvl[v] = -1;
    }
    while (front < back) {
      int v = q[front++]; for (int w : G[v]) {
        if (mate[w] == -1) return true;
        else if (lvl[mate[w]] == -1) lvl[q[back++] = mate[w]] = lvl[v] + 1;
      }
    }
    return false;
  }
  template <class BipartiteGraph> bool dfs(const BipartiteGraph &G, int v) {
    for (int w : G[v])
      if (mate[w] == -1 || (lvl[mate[w]] == lvl[v] + 1 && dfs(G, mate[w]))) {
        mate[mate[v] = w] = v; return true;
      }
    lvl[v] = -1; return false;
  }
  template <class BipartiteGraph>
  void dfsVertexCover(const BipartiteGraph &G, int v) {
    if (vis[v]) return;
    vis[v] = true;
    for (int w : G[v]) if ((mate[v] == w) == color[v]) dfsVertexCover(G, w);
  }
  template <class BipartiteGraph>
  HopcroftKarpMaxMatch(const BipartiteGraph &G, const vector<bool> &color)
      : V(G.size()), cardinality(0), mate(V, -1), lvl(V), q(V),
        color(color), inCover(V, false), vis(V, false) {
    for (int v = 0; v < V; v++) if (!color[v]) type0.push_back(v);
    while (bfs(G)) for (int v : type0) if (mate[v] == -1 && dfs(G, v))
      cardinality++;
    for (int v : type0) if (mate[v] == -1) dfsVertexCover(G, v);
    for (int v = 0; v < V; v++) inCover[v] = vis[v] == color[v];
  }
};


text:
#pragma once
#include <bits/stdc++.h>
using namespace std;

// Floyd Warshall's all pairs shortest path algorithm for weighted graphs
//   with negative weight
// Able to detect negative cycles
// Vertices are 0-indexed
// Template Arguments:
//   T: the type of the weight of the edges in the graph
// Constructor Arguments:
//   matrix: a V x V matrix containing the minimum weight of a directed edge
//     between two vertices, INF if an edge doesn't exist
//   INF: a value for infinity, must be negatable
// Fields:
//   dist: a vector of vectors of the shortest distance between each pair
//     of vertices, INF if there is no path, -INF if the shortest path
//     has no lower bound
//   par: a vector of vectors of the parent vertex for each vertex in the
//     shortest path tree for each source vertex (par[v][w] is the parent
//     of vertex w in the shortest path tree from vertex v), or -1 if there is
//     no parent
//   hasNegativeCycle: a boolean that is true if there is a negative cycle
//     in the graph and false otherwise
// Functions:
//   getPath(v, w): returns the list of directed edges on the path from
//     vertex v to vertex w
// In practice, has a very small constant
// Time Complexity:
//   constructor: O(V^3)
//   getPath: O(V)
// Memory Complexity: O(V^2)
// Tested:
//   Fuzz and Stress Tested
//   https://open.kattis.com/problems/allpairspath
template <class T> struct FloydWarshallAPSP {
  using Edge = tuple<int, int, T>; vector<vector<T>> dist;
  vector<vector<int>> par; T INF; bool hasNegativeCycle;
  FloydWarshallAPSP(const vector<vector<T>> &matrix,
                    T INF = numeric_limits<T>::max())
      : dist(matrix), par(dist.size(), vector<int>(dist.size(), -1)),
        INF(INF), hasNegativeCycle(false) {
    int V = dist.size(); for (int v = 0; v < V; v++) {
      for (int w = 0; w < V; w++) if (dist[v][w] < INF) par[v][w] = v;
      dist[v][v] = min(T(), dist[v][v]); par[v][v] = -1;
    }
    for (int u = 0; u < V; u++) for (int v = 0; v < V; v++)
      if (dist[v][u] < INF) for (int w = 0; w < V; w++)
        if (dist[u][w] < INF && dist[v][w] > dist[v][u] + dist[u][w]) {
          dist[v][w] = dist[v][u] + dist[u][w]; par[v][w] = par[u][w];
        }
    for (int u = 0; u < V; u++) for (int v = 0; v < V; v++)
      for (int w = 0; w < V; w++)
        if (dist[w][w] < T() && dist[u][w] < INF && dist[w][v] < INF) {
          dist[u][v] = -INF; hasNegativeCycle = true; break;
        }
  }
  vector<Edge> getPath(int v, int w) {
    vector<Edge> path; for (; par[v][w] != -1; w = par[v][w])
      path.emplace_back(par[v][w], w, dist[v][w] - dist[v][par[v][w]]);
    reverse(path.begin(), path.end()); return path;
  }
};


text:
#pragma once
#include <bits/stdc++.h>
using namespace std;

// Bellman Ford's single source shortest path algorithm for weighted graphs
//   with negative weights
// Able to detect negative cycles
// Vertices are 0-indexed
// Template Arguments:
//   T: the type of the weight of the edges in the graph
// Constructor Arguments:
//   V: the number of vertices in the graph
//   edges: a vector of tuples in the form (v, w, weight) representing
//     a directed edge in the graph from v to w with weight of weight
//   s: a single source vertex
//   src: a vector of source vertices
//   INF: a value for infinity, must be negatable
// Fields:
//   dist: vector of shortest distance from the closest source vertex to each
//     vertex, or INF if there is no path, -INF if the shortest path
//     has no lower bound
//   par: the parent vertex for each vertex in the shortest path tree, or
//     -1 if there is no parent
//   hasNegativeCycle: a boolean that is true if there is a negative cycle
//     in the graph and false otherwise
// Functions:
//   getPath(v): returns the list of directed edges on the path from the
//     closest source vertex to vertex v
// In practice, has a very small constant
// Time Complexity:
//   constructor: O(E (V + log E))
//   getPath: O(V)
// Memory Complexity: O(V + E)
// Tested:
//   Stress Tested
//   https://open.kattis.com/problems/shortestpath3
//   https://dmoj.ca/problem/sssp
template <class T> struct BellmanFordSSSP {
  using Edge = tuple<int, int, T>; vector<T> dist; vector<int> par; T INF;
  bool hasNegativeCycle;
  BellmanFordSSSP(int V, vector<Edge> edges, const vector<int> &srcs,
                  T INF = numeric_limits<T>::max())
      : dist(V, INF), par(V, -1), INF(INF), hasNegativeCycle(false) {
    sort(edges.begin(), edges.end()); for (int s : srcs) dist[s] = T();
    for (int i = 0; i < V - 1; i++) for (auto &&e : edges) {
      int v, w; T weight; tie(v, w, weight) = e;
      if (dist[v] < INF && dist[w] > dist[v] + weight)
        dist[w] = dist[par[w] = v] + weight;
    }
    for (bool inCycle = true; inCycle;) {
      inCycle = false; for (auto &&e : edges) {
        int v, w; T weight; tie(v, w, weight) = e;
        if (dist[v] < INF && dist[w] > -INF
            && (dist[v] <= -INF || dist[w] > dist[v] + weight)) {
          dist[w] = -INF; inCycle = hasNegativeCycle = true;
        }
      }
    }
  }
  BellmanFordSSSP(int V, const vector<Edge> &edges, int s,
                  T INF = numeric_limits<T>::max())
      : BellmanFordSSSP(V, edges, vector<int>{s}, INF) {}
  vector<Edge> getPath(int v) {
    vector<Edge> path; for (; par[v] != -1; v = par[v])
      path.emplace_back(par[v], v, dist[v] - dist[par[v]]);
    reverse(path.begin(), path.end()); return path;
  }
};


text:
#pragma once
#include <bits/stdc++.h>
#include "../../datastructures/trees/heaps/PersistentRandomizedHeap.h"
using namespace std;

// Computes the K shortest walks from s to t in a directed graph
// Template Arguments:
//   T: the type of the weight of the edges
// Function Arguments:
//   V: number of vertices in the directed graph
//   edges: a vector of tuples in the form (v, w, weight) representing
//     a directed edge in the graph from vertex v to w with
//     weight of weight
//   K: the number of shortest walks to compute
//   s: the source vertex
//   t: the destination vertex
//   INF: a value for infinity
// Return Value: a vector of type T of size K with the K shortest walks
// In practice, has a moderate constant
// Time Complexity: O((V + E + K) log V)
// Memory Complexity: O((V + E) log V + K)
// Tested:
//   https://judge.yosupo.jp/problem/k_shortest_walk
template <class T>
vector<T> kShortestWalks(int V, const vector<tuple<int, int, T>> &edges, int K,
                         int s, int t, T INF = numeric_limits<T>::max()) {
  if (K == 0) return vector<T>();
  using Heap = PersistentRandomizedHeap<pair<T, int>, greater<pair<T, int>>>;
  using ptr = typename Heap::ptr; vector<Heap> cand(V);
  struct Node1 {
    T d; int v; Node1(T d, int v) : d(d), v(v) {}
    bool operator < (const Node1 &o) const { return d > o.d; }
  };
  struct Node2 {
    T d; int v; ptr p; Node2(T d, int v, ptr p) : d(d), v(v), p(p) {}
    bool operator < (const Node2 &o) const { return d > o.d; }
  };
  vector<T> dist(V, INF), ret; ret.reserve(K); vector<int> ord; ord.reserve(V);
  vector<pair<int, int>> par(V, make_pair(-1, -1));
  vector<vector<tuple<int, T, int>>> G(V), H(V);
  std::priority_queue<Node1> PQ;
  std::priority_queue<Node2> ans;
  for (int i = 0; i < int(edges.size()); i++) {
    int v, w; T weight; tie(v, w, weight) = edges[i];
    G[v].emplace_back(w, weight, i); H[w].emplace_back(v, weight, i);
  }
  PQ.emplace(dist[t] = T(), t); while (!PQ.empty()) {
    T d = PQ.top().d; int v = PQ.top().v; PQ.pop(); if (d > dist[v]) continue;
    ord.push_back(v);
    for (auto &&e : H[v]) if (dist[get<0>(e)] > dist[v] + get<1>(e)) {
      PQ.emplace(dist[get<0>(e)] = dist[v] + get<1>(e), get<0>(e));
      par[get<0>(e)] = make_pair(get<2>(e), v);
    }
  }
  for (int v : ord) {
    for (auto &&e : G[v])
      if (get<2>(e) != par[v].first && dist[get<0>(e)] < INF) {
        T d = dist[get<0>(e)] - dist[v] + get<1>(e);
        cand[v].push(make_pair(d, get<0>(e)));
      }
    if (par[v].first != -1) cand[v].merge(cand[par[v].second]);
    if (v == s) {
      ret.push_back(dist[v]); if (cand[v].root)
        ans.emplace(dist[v] + cand[v].top().first, v, cand[v].root);
    }
  }
  while (int(ret.size()) < K) {
    if (ans.empty()) { ret.push_back(INF); continue; }
    T d = ans.top().d; int v = ans.top().v; ptr p = ans.top().p; ans.pop();
    ret.push_back(d);
    if (p->l) ans.emplace(d + p->l->val.first - p->val.first, v, p->l);
    if (p->r) ans.emplace(d + p->r->val.first - p->val.first, v, p->r);
    v = p->val.second;
    if (cand[v].root) ans.emplace(d + cand[v].top().first, v, cand[v].root);
  }
  return ret;
}


text:
#pragma once
#include <bits/stdc++.h>
using namespace std;

// Computes the shortest Hamiltonian cycle
// All vertices are visited exactly once before returning to the start vertex
// Vertices are 0-indexed
// Template Arguments:
//   T: the type of the weight of the edges in the graph
// Constructor Arguments:
//   matrix: a V x V matrix containing the minimum weight of a directed edge
//     between two vertices, INF if an edge doesn't exist
// Fields:
//   shortestCycleDist: the distance of the shortest Hamiltonian cycle, INF if
//     no cycle exists (V == 0 implies shortestCycleDist == 0,
//     V == 1 implies shortestCycleDist == INF)
//   ord: the order the vertices are visited in the shortest Hamiltonian cycle,
//     the start vertex is not repeated
// In practice, has a small constant
// Time Complexity:
//   constructor: O(V^2 2^V)
// Memory Complexity: O(V 2^V)
// Tested:
//   Fuzz Tested
template <class T> struct ShortestHamiltonianCycle {
  vector<vector<T>> dp; T INF, shortestCycleDist; vector<int> ord;
  ShortestHamiltonianCycle(const vector<vector<T>> &matrix,
                           T INF = numeric_limits<T>::max())
      : dp(1 << matrix.size(), vector<T>(matrix.size(), INF)), INF(INF),
        shortestCycleDist(matrix.empty() ? T() : INF), ord(matrix.size(), 0) {
    int V = matrix.size(); if (V > 0) dp[1][0] = T();
    for (int mask = 1; mask < (1 << V); mask += 2) for (int i = 1; i < V; i++)
      if ((mask >> i) & 1) {
        for (int o = mask ^ (1 << i), j = 0; j < V; j++) if ((mask >> j) & 1) {
          if (matrix[j][i] < INF && dp[o][j] < INF
              && dp[mask][i] > dp[o][j] + matrix[j][i])
            dp[mask][i] = dp[o][j] + matrix[j][i];
        }
      }
    int cur = (1 << V) - 1;
    for (int i = 1; i < V; i++) if (dp.back()[i] < INF && matrix[i][0] < INF)
      shortestCycleDist = min(shortestCycleDist, dp.back()[i] + matrix[i][0]);
    if (shortestCycleDist >= INF) return;
    auto get = [&] (int mask, int j, int last) {
      if (dp[mask][j] >= INF || matrix[j][last] >= INF) return INF;
      return dp[mask][j] + matrix[j][last];
    };
    for (int last = 0, i = V - 2; i >= 0; i--) {
      int bj = -1; for (int j = 1; j < V; j++) if ((cur >> j) & 1)
        if (bj == -1 || get(cur, bj, last) > get(cur, j, last)) bj = j;
      cur ^= 1 << (last = ord[i] = bj);
    }
  }
};


text:
#pragma once
#include <bits/stdc++.h>
#include "../search/TransitiveClosureSCC.h"
using namespace std;

// Johnson's all pairs shortest path algorithm for weighted graphs
//   with negative weights
// Able to detect negative cycles
// Vertices are 0-indexed
// Template Arguments:
//   T: the type of the weight of the edges in the graph
// Constructor Arguments:
//   G: a generic weighted graph structure
//     Required Functions:
//       operator [v] const: iterates over the adjacency list of vertex v
//         (which is a list of pair<int, T> with weights of type T)
//       size() const: returns the number of vertices in the graph
//   MAXV: the maximum number of vertices in the graph
//   INF: a value for infinity, must be negatable
// Fields:
//   dist: a vector of vectors of the shortest distance between each pair
//     of vertices, INF if there is no path, -INF if the shortest path
//     has no lower bound
//   par: a vector of vectors of the parent vertex for each vertex in the
//     shortest path tree for each source vertex (par[v][w] is the parent
//     of vertex w in the shortest path tree from vertex v), or -1 if there is
//     no parent
//   hasNegativeCycle: a boolean that is true if there is a negative cycle
//     in the graph and false otherwise
// Functions:
//   getPath(v, w): returns the list of directed edges on the path from
//     vertex v to vertex w
// In practice, has a small constant
// Time Complexity:
//   constructor: O(VE log E)
//   getPath: O(V)
// Memory Complexity: O(V + E + MAXV V / 64)
// Tested:
//   https://dmoj.ca/problem/apsp
template <class T, const int MAXV> struct JohnsonAPSP {
  using Edge = tuple<int, int, T>; int V; vector<vector<T>> dist;
  vector<vector<int>> par; T INF; bool hasNegativeCycle;
  struct Node {
    T d; int v; Node(T d, int v) : d(d), v(v) {}
    bool operator < (const Node &o) const { return d > o.d; }
  };
  template <class WeightedGraph>
  JohnsonAPSP(const WeightedGraph &G, T INF = numeric_limits<T>::max())
      : V(G.size()), dist(V, vector<T>(V, INF)), par(V, vector<int>(V, -1)),
        INF(INF), hasNegativeCycle(false) {
    bool hasNegativeWeight = false; for (int v = 0; v < V; v++)
      for (auto &&e : G[v]) hasNegativeWeight |= e.second < T();
    vector<bool> inNegCyc(V, false); vector<T> h(V, T());
    vector<int> id; vector<bitset<MAXV>> neg; if (hasNegativeWeight) {
      TransitiveClosureSCC<MAXV> tc(G); id = tc.scc.id;
      vector<bool> compInNegCyc(tc.dp.size(), false);
      for (int i = 0; i < V - 1; i++) for (int v = 0; v < V; v++)
        for (auto &&e : G[v]) if (id[v] == id[e.first])
          h[e.first] = min(h[e.first], h[v] + e.second);
      for (int v = 0; v < V; v++) for (auto &&e : G[v])
        if (id[v] == id[e.first] && h[e.first] > h[v] + e.second)
          compInNegCyc[id[v]] = true;
      for (int v = 0; v < V; v++)
        hasNegativeCycle |= (inNegCyc[v] = compInNegCyc[id[v]]);
      fill(h.begin(), h.end(), T());
      for (int i = 0; i < V - 1; i++) for (int v = 0; v < V; v++)
        for (auto &&e : G[v]) if (!inNegCyc[v] && !inNegCyc[e.first])
          h[e.first] = min(h[e.first], h[v] + e.second);
      neg.resize(tc.dp.size()); for (auto &&e : tc.DAG) {
        if (compInNegCyc[e.second]) neg[e.first] |= tc.dp[e.second];
        else neg[e.first] |= neg[e.second];
      }
      for (int i = 0; i < int(neg.size()); i++) if (compInNegCyc[i])
        neg[i] = tc.dp[i];
    }
    for (int s = 0; s < V; s++) {
      if (!inNegCyc[s]) {
        std::priority_queue<Node> PQ;
        PQ.emplace(dist[s][s] = T(), s); while (!PQ.empty()) {
          T d = PQ.top().d; int v = PQ.top().v; PQ.pop();
          if (d > dist[s][v]) continue;
          for (auto &&e : G[v]) {
            int w = e.first; T weight = e.second + h[v] - h[w];
            if (!inNegCyc[w] && dist[s][w] > dist[s][v] + weight)
              PQ.emplace(dist[s][w] = dist[s][par[s][w] = v] + weight, w);
          }
        }
      }
      if (hasNegativeWeight) for (int v = 0; v < V; v++)
        if (neg[id[s]][id[v]]) dist[s][v] = -INF;
    }
    if (hasNegativeWeight) for (int v = 0; v < V; v++)
      for (int w = 0; w < V; w++) if (dist[v][w] != INF && dist[v][w] != -INF)
        dist[v][w] = dist[v][w] - h[v] + h[w];
  }
  vector<Edge> getPath(int v, int w) {
    vector<Edge> path; for (; par[v][w] != -1; w = par[v][w])
      path.emplace_back(par[v][w], w, dist[v][w] - dist[v][par[v][w]]);
    reverse(path.begin(), path.end()); return path;
  }
};


text:
#pragma once
#include <bits/stdc++.h>
#include "../search/TopologicalOrder.h"
using namespace std;

// Computes the shortest path (based on the comparator) on a directed
//   acyclic graph
// Vertices are 0-indexed
// Template Arguments:
//   T: the type of the weight of the edges in the graph
//   Cmp: the comparator to compare two distances; less<T> will compute
//       the shortest path while greater<T> will compute the longest path
//     Required Functions:
//       operator (a, b): returns true if and only if a compares less than b
// Constructor Arguments:
//   G: a generic directed acyclic graph structure (weighted or unweighted)
//     Required Functions:
//       operator [v] const: iterates over the adjacency list of vertex v
//         (which is a list of ints for an unweighted graph, or a list of
//         pair<int, T> for a weighted graph with weights of type T)
//       size() const: returns the number of vertices in the graph
//   s: a single source vertex
//   src: a vector of source vertices
//   INF: a value for infinity for shortest path, or negative infinity
//     for longest path; (Cmp()(INF, x)) must return false for all
//     values of x
//   cmp: an instance of the Cmp struct
// Fields:
//   dist: vector of distance from the closest source vertex to each vertex,
//     or INF if unreachable, and is also the shortest distance for
//     an unweighted graph
//   par: the parent vertex for each vertex in the shortest path tree,
//     or -1 if there is no parent
// Functions:
//   getPath(v): returns the list of directed edges on the path from the
//     closest source vertex to vertex v
// In practice, has a moderate constant
// Time Complexity:
//   constructor: O(V + E)
//   getPath: O(V)
// Memory Complexity: O(V)
// Tested:
//   https://atcoder.jp/contests/dp/tasks/dp_g
template <class T, class Cmp = less<T>> struct DAGSP {
  using Edge = tuple<int, int, T>; vector<T> dist; vector<int> par; T INF;
  TopologicalOrder ord;
  int getTo(int e) { return e; }
  T getWeight(int) { return 1; }
  int getTo(const pair<int, T> &e) { return e.first; }
  T getWeight(const pair<int, T> &e) { return e.second; }
  template <class DAG>
  DAGSP(const DAG &G, const vector<int> &srcs, T INF, Cmp cmp = Cmp())
      : dist(G.size(), INF), par(G.size(), -1), INF(INF), ord(G) {
    for (int s : srcs) dist[s] = T();
    for (int v : ord.ord) for (auto &&e : G[v]) {
      int w = getTo(e); T weight = getWeight(e);
      if (cmp(dist[v], INF) && cmp(dist[v] + weight, dist[w]))
        dist[w] = dist[par[w] = v] + weight;
    }
  }
  template <class DAG> DAGSP(const DAG &G, int s, T INF, Cmp cmp = Cmp())
      : DAGSP(G, vector<int>{s}, INF, cmp) {}
  vector<Edge> getPath(int v) {
    vector<Edge> path; for (; par[v] != -1; v = par[v])
      path.emplace_back(par[v], v, dist[v] - dist[par[v]]);
    reverse(path.begin(), path.end()); return path;
  }
};


text:
#pragma once
#include <bits/stdc++.h>
using namespace std;

// Classical Dijkstra's single source shortest path algorithm for
//   weighted graphs without negative weights
// Vertices are 0-indexed
// Template Arguments:
//   T: the type of the weight of the edges in the graph
// Constructor Arguments:
//   G: a generic weighted graph structure
//     Required Functions:
//       operator [v] const: iterates over the adjacency list of vertex v
//         (which is a list of pair<int, T> with weights of type T)
//       size() const: returns the number of vertices in the graph
//   s: a single source vertex
//   src: a vector of source vertices
//   INF: a value for infinity
// Fields:
//   dist: vector of shortest distance from the closest source vertex to each
//     vertex, or INF if unreachable
//   par: the parent vertex for each vertex in the shortest path tree, or
//     -1 if there is no parent
// Functions:
//   getPath(v): returns the list of directed edges on the path from the
//     closest source vertex to vertex v
// In practice, has a small constant
// Time Complexity:
//   constructor: O(V^2 + E)
//   getPath: O(V)
// Memory Complexity: O(V)
// Tested:
//   Stress Tested
//   https://dmoj.ca/problem/sssp
template <class T> struct ClassicalDijkstraSSSP {
  using Edge = tuple<int, int, T>; vector<T> dist; vector<int> par; T INF;
  template <class WeightedGraph>
  ClassicalDijkstraSSSP(const WeightedGraph &G, const vector<int> &srcs,
                        T INF = numeric_limits<T>::max())
      : dist(G.size(), INF), par(G.size(), -1), INF(INF) {
    vector<bool> done(G.size(), false); for (int s : srcs) dist[s] = T();
    for (int i = 0; i < int(G.size()) - 1; i++) {
      int v = -1; for (int w = 0; w < int(G.size()); w++)
        if (!done[w] && (v == -1 || dist[v] > dist[w])) v = w;
      if (dist[v] >= INF) break;
      done[v] = true;
      for (auto &&e : G[v]) if (dist[e.first] > dist[v] + e.second)
        dist[e.first] = dist[par[e.first] = v] + e.second;
    }
  }
  template <class WeightedGraph>
  ClassicalDijkstraSSSP(const WeightedGraph &G, int s,
                        T INF = numeric_limits<T>::max())
      : ClassicalDijkstraSSSP(G, vector<int>{s}, INF) {}
  vector<Edge> getPath(int v) {
    vector<Edge> path; for (; par[v] != -1; v = par[v])
      path.emplace_back(par[v], v, dist[v] - dist[par[v]]);
    reverse(path.begin(), path.end()); return path;
  }
};


text:
#pragma once
#include <bits/stdc++.h>
using namespace std;

// Computes the shortest Hamiltonian path
// All vertices are visited exactly once, does not return to the start
// Vertices are 0-indexed
// Template Arguments:
//   T: the type of the weight of the edges in the graph
// Constructor Arguments:
//   matrix: a V x V matrix containing the minimum weight of a directed edge
//     between two vertices, INF if an edge doesn't exist
// Fields:
//   shortestPathDist: the distance of the shortest Hamiltonian path, INF if
//     no path exists (V <= 1 implies shortestPathDist == 0)
//   ord: the order the vertices are visited in the shortest Hamiltonian path
// In practice, has a small constant
// Time Complexity:
//   constructor: O(V^2 2^V)
// Memory Complexity: O(V 2^V)
// Tested:
//   Fuzz Tested
template <class T> struct ShortestHamiltonianPath {
  vector<vector<T>> dp; T INF, shortestPathDist; vector<int> ord;
  ShortestHamiltonianPath(const vector<vector<T>> &matrix,
                          T INF = numeric_limits<T>::max())
      : dp(1 << matrix.size(), vector<T>(matrix.size(), INF)), INF(INF),
        shortestPathDist(matrix.empty() ? T() : INF), ord(matrix.size(), 0) {
    int V = matrix.size(); for (int i = 0; i < V; i++) dp[1 << i][i] = T();
    for (int mask = 0; mask < (1 << V); mask++) for (int i = 0; i < V; i++)
      if ((mask >> i) & 1) {
        for (int o = mask ^ (1 << i), j = 0; j < V; j++) if ((mask >> j) & 1) {
          if (matrix[j][i] < INF && dp[o][j] < INF
              && dp[mask][i] > dp[o][j] + matrix[j][i])
            dp[mask][i] = dp[o][j] + matrix[j][i];
        }
      }
    int cur = (1 << V) - 1; for (int i = 0; i < V; i++)
      shortestPathDist = min(shortestPathDist, dp.back()[i]);
    if (shortestPathDist >= INF) return;
    auto get = [&] (int mask, int j, int last) {
      T weight = (last == -1 ? T() : matrix[j][last]);
      if (dp[mask][j] < INF && weight < INF) return dp[mask][j] + weight;
      else return INF;
    };
    for (int last = -1, i = V - 1; i >= 0; i--) {
      int bj = -1; for (int j = 0; j < V; j++) if ((cur >> j) & 1)
        if (bj == -1 || get(cur, bj, last) > get(cur, j, last)) bj = j;
      cur ^= 1 << (last = ord[i] = bj);
    }
  }
};


text:
#pragma once
#include <bits/stdc++.h>
using namespace std;

// Dijkstra's single source shortest path algorithm for weighted graphs
//   without negative weights
// Vertices are 0-indexed
// Template Arguments:
//   T: the type of the weight of the edges in the graph
// Constructor Arguments:
//   G: a generic weighted graph structure
//     Required Functions:
//       operator [v] const: iterates over the adjacency list of vertex v
//         (which is a list of pair<int, T> with weights of type T)
//       size() const: returns the number of vertices in the graph
//   s: a single source vertex
//   src: a vector of source vertices
//   INF: a value for infinity
// Fields:
//   dist: vector of shortest distance from the closest source vertex to each
//     vertex, or INF if unreachable
//   par: the parent vertex for each vertex in the shortest path tree, or
//     -1 if there is no parent
// Functions:
//   getPath(v): returns the list of directed edges on the path from the
//     closest source vertex to vertex v
// In practice, has a small constant
// Time Complexity:
//   constructor: O((V + E) log E)
//   getPath: O(V)
// Memory Complexity: O(V + E)
// Tested:
//   Stress Tested
//   https://judge.yosupo.jp/problem/shortest_path
//   https://open.kattis.com/problems/shortestpath1
//   https://dmoj.ca/problem/sssp
template <class T> struct DijkstraSSSP {
  using Edge = tuple<int, int, T>; vector<T> dist; vector<int> par; T INF;
  struct Node {
    T d; int v; Node(T d, int v) : d(d), v(v) {}
    bool operator < (const Node &o) const { return d > o.d; }
  };
  template <class WeightedGraph>
  DijkstraSSSP(const WeightedGraph &G, const vector<int> &srcs,
               T INF = numeric_limits<T>::max())
      : dist(G.size(), INF), par(G.size(), -1), INF(INF) {
    std::priority_queue<Node> PQ;
    for (int s : srcs) PQ.emplace(dist[s] = T(), s);
    while (!PQ.empty()) {
      T d = PQ.top().d; int v = PQ.top().v; PQ.pop();
      if (d > dist[v]) continue;
      for (auto &&e : G[v]) if (dist[e.first] > dist[v] + e.second)
        PQ.emplace(dist[e.first] = dist[par[e.first] = v] + e.second, e.first);
    }
  }
  template <class WeightedGraph> DijkstraSSSP(const WeightedGraph &G, int s,
                                              T INF = numeric_limits<T>::max())
      : DijkstraSSSP(G, vector<int>{s}, INF) {}
  vector<Edge> getPath(int v) {
    vector<Edge> path; for (; par[v] != -1; v = par[v])
      path.emplace_back(par[v], v, dist[v] - dist[par[v]]);
    reverse(path.begin(), path.end()); return path;
  }
};


text:
#pragma once
#include <bits/stdc++.h>
#include "Point3D.h"
#include "Sphere3D.h"
using namespace std;

// Functions for a 3D polyhedron

// Determines twice the vector area of a face of coplanar points
// Function Arguments:
//   face: the coplanar points of the face
// Return Value: twice the vector area of the face, points are in ccw order
//   when looking from above (where up is the direction of the return vector)
// Time Complexity: O(N)
// Memory Complexity: O(1)
// Tested:
//   https://www.spoj.com/problems/CH3D/
//   https://open.kattis.com/problems/threedprinter
pt3 vectorArea2(const vector<pt3> &face) {
  pt3 ret(0, 0, 0); for (int i = 0, n = face.size(); i < n; i++)
    ret += face[i] * face[(i + 1) % n];
  return ret;
}

// Returns the absolute area of a face of coplanar points
// Function Arguments:
//   face: the coplanar points of the face
// Return Value: the area of the face
// Time Complexity: O(N)
// Memory Complexity: O(1)
// Tested:
//   https://www.spoj.com/problems/CH3D/
T faceArea(const vector<pt3> &face) { return abs(vectorArea2(face)) / T(2); }

// Flips some of the faces of a polyhedron such that they are oriented
//   consistently (they will either all be pointed outwards or inwards)
// Function Arguments:
//   faces: a reference to a vector of vectors of points representing
//     the faces of the polyhedron to be reoriented
// Time Complexity: O(N) for N total points
// Memory Complexity: O(1)
// Tested:
//   https://open.kattis.com/problems/threedprinter
void reorient(vector<vector<pt3>> &faces) {
  int n = faces.size(); vector<vector<pair<int, bool>>> G(n);
  map<pair<pt3, pt3>, int> seen; for (int v = 0; v < n; v++)
    for (int i = 0, m = faces[v].size(); i < m; i++) {
      pt3 a = faces[v][i], b = faces[v][(i + 1) % m];
      auto it = seen.find(make_pair(a, b)); bool f = true;
      if (it == seen.end()) { it = seen.find(make_pair(b, a)); f = false; }
      if (it == seen.end()) seen[make_pair(a, b)] = v;
      else {
        int w = it->second; G[v].emplace_back(w, f); G[w].emplace_back(v, f);
      }
    }
  vector<char> flip(n, -1); vector<int> q(n); int front = 0, back = 0;
  flip[q[back++] = 0] = 0; while (front < back) {
    int v = q[front++]; for (auto &&e : G[v]) if (flip[e.first] == -1)
      flip[q[back++] = e.first] = flip[v] ^ e.second;
  }
  for (int v = 0; v < n; v++) {
    assert(flip[v] != -1);
    if (flip[v]) reverse(faces[v].begin(), faces[v].end());
  }
}

// Returns the surface area of a polyhedron
// Function Arguments:
//   faces: a vector of vectors of points representing
//     the faces of the polyhedron with consistent orientation
// Return Value: the surface area of the polyhedron
// Time Complexity: O(N) for N total points
// Memory Complexity: O(1)
// Tested:
//   https://www.spoj.com/problems/CH3D/
T getSurfaceArea(const vector<vector<pt3>> &faces) {
  T sa = 0; for (auto &&face : faces) sa += faceArea(face);
  return sa;
}

// Returns 6 times the signed volume of a polyhedron
// Function Arguments:
//   faces: a vector of vectors of points representing
//     the faces of the polyhedron with consistent orientation
// Return Value: 6 times the signed volume of the polyhedron, positive if
//   all vector areas point outwards, negative if inwards
// Time Complexity: O(N) for N total points
// Memory Complexity: O(1)
// Tested:
//   https://www.spoj.com/problems/CH3D/
//   https://open.kattis.com/problems/threedprinter
T getVolume6(const vector<vector<pt3>> &faces) {
  T vol6 = 0; for (auto &&face : faces) vol6 += (vectorArea2(face) | face[0]);
  return vol6;
}

// Determines if a point is inside a polyhedron or not
// Function Arguments:
//   faces: a vector of vectors of points representing
//     the faces of the polyhedron with consistent orientation
//   p: the point to check
// Return Value: 1 if inside the polyhedron, 0 if on the face, -1 if outside
// Time Complexity: O(N) for N total points
// Memory Complexity: O(1)
int isInPolyhedron(const vector<vector<pt3>> &faces, pt3 p) {
  T sum = 0, PI = acos(T(-1)); Sphere3D s(p, 1); for (auto &&face : faces) {
    pt3 a = face[0], b = face[1], c = face[2], n = (b - a) * (c - a);
    if (eq((n | p) - (n | a), 0)) return 0;
    sum += remainder(s.surfaceAreaOnSph(face), 4 * PI);
  }
  return eq(sum, 0) ? -1 : 1;
}


text:
#pragma once
#include <bits/stdc++.h>
#include "../../utils/EpsCmp.h"
using namespace std;

// Functions for a 3D point
// * operator between 2 points is cross product
// | operator between 2 points is dot product
#define OP(op, U, a, x, y, z) \
  pt3 operator op (U a) const { return pt3(x, y, z); } \
  pt3 &operator op##= (U a) { return *this = *this op a; }
#define CMP(op, body) bool operator op (pt3 p) const { return body; }
struct pt3 {
  T x, y, z; constexpr pt3(T x = 0, T y = 0, T z = 0) : x(x), y(y), z(z) {}
  pt3 operator + () const { return *this; }
  pt3 operator - () const { return pt3(-x, -y, -z); }
  OP(+, pt3, p, x + p.x, y + p.y, z + p.z)
  OP(-, pt3, p, x - p.x, y - p.y, z - p.z)
  OP(*, T, a, x * a, y * a, z * a) OP(/, T, a, x / a, y / a, z / a)
  friend pt3 operator * (T a, pt3 p) { return pt3(a * p.x, a * p.y, a * p.z); }
  bool operator < (pt3 p) const {
    return eq(x, p.x) ? (eq(y, p.y) ? lt(z, p.z) : lt(y, p.y)) : lt(x, p.x);
  }
  CMP(<=, !(p < *this)) CMP(>, p < *this) CMP(>=, !(*this < p))
  CMP(==, !(*this < p) && !(p < *this)) CMP(!=, *this < p || p < *this)
  T operator | (pt3 p) const { return x * p.x + y * p.y + z * p.z; }
  OP(*, pt3, p, y * p.z - z * p.y, z * p.x - x * p.z, x * p.y - y * p.x)
};
#undef OP
#undef CMP
istream &operator >> (istream &stream, pt3 &p) {
  return stream >> p.x >> p.y >> p.z;
}
ostream &operator << (ostream &stream, pt3 p) {
  return stream << p.x << ' ' << p.y << ' ' << p.z;
}
T norm(pt3 p) { return p | p; }
T abs(pt3 p) { return sqrt(norm(p)); }
pt3 unit(pt3 p) { return p / abs(p); }
T distSq(pt3 a, pt3 b) { return norm(b - a); }
T dist(pt3 a, pt3 b) { return abs(b - a); }
// returns an angle in the range [0, PI]
T ang(pt3 a, pt3 b, pt3 c) {
  a = unit(a - b); c = unit(c - b); return 2 * atan2(abs(a - c), abs(a + c));
}
pt3 rot(pt3 a, pt3 axis, T theta) {
  return a * cos(theta) + (unit(axis) * a * sin(theta))
      + (unit(axis) * (unit(axis) | a) * (1 - cos(theta)));
}
// sign of volume6 and above: 1 if d is above the plane abc with
//   normal ab x ac, 0 if on the plane, -1 if below the plane
// 6 times the signed area of the tetrahedron abcd
T volume6(pt3 a, pt3 b, pt3 c, pt3 d) {
  return (b - a) * (c - a) | (d - a);
}
int above(pt3 a, pt3 b, pt3 c, pt3 d) { return sgn(volume6(a, b, c, d)); }
// Converts a position based on radius (r >= 0), inclination/latitude
//   (-PI / 2 <= theta <= PI / 2), and azimuth/longitude (-PI < phi <= PI)
// Convention is that the x axis passes through the meridian (phi = 0), and
//   the z axis passes through the North Pole (theta = Pi / 2)
pt3 sph(T r, T theta, T phi) {
  return pt3(r * cos(theta) * cos(phi), r * cos(theta) * sin(phi),
             r * sin(theta));
}
T inc(pt3 p) { return atan2(p.z, T(sqrt(p.x * p.x + p.y * p.y))); }
T az(pt3 p) { return atan2(p.y, p.x); }


text:
#pragma once
#include <bits/stdc++.h>
#include "../../utils/Random.h"
#include "Point3D.h"
using namespace std;

// Computes the convex hull of a set of N points 3D points (convex set of
//   minimum points with the minimum volume)
// Function Arguments:
//   P: the vector of points
// Return Value: a vector of faces in the convex hull, with each face being a
//   vector of exactly 3 points, all facing outwards
// In practice, has a large constant
// Time Complexity: O(N log N) expected
// Memory Complexity: O(N log N)
// Tested:
//   https://www.spoj.com/problems/CH3D/
//   https://open.kattis.com/problems/worminapple
//   https://open.kattis.com/problems/starsinacan
vector<vector<pt3>> convexHull3D(vector<pt3> P) {
  vector<array<int, 3>> hullInd; shuffle(P.begin(), P.end(), rng);
  int n = P.size(); for (int i = 1, num = 1; i < n; i++) {
    if (num == 1) {
      if (P[0] != P[i]) { swap(P[1], P[i]); num++; }
    } else if (num == 2) {
      if (!eq(norm((P[1] - P[0]) * (P[i] - P[0])), 0)) {
        swap(P[2], P[i]); num++;
      }
    } else if (num == 3) {
      if (above(P[0], P[1], P[2], P[i]) != 0) { swap(P[3], P[i]); num++; }
    }
  }
  vector<bool> active; vector<vector<int>> vis(n), rvis;
  vector<array<pair<int, int>, 3>> other; vector<int> label(n, -1);
  auto addFace = [&] (int a, int b, int c) {
    hullInd.push_back({a, b, c}); active.push_back(true);
    rvis.emplace_back(); other.emplace_back();
  };
  auto addEdge = [&] (int a, int b) {
    vis[b].push_back(a); rvis[a].push_back(b);
  };
  auto abv = [&] (int a, int b) {
    array<int, 3> f = hullInd[a];
    return above(P[f[0]], P[f[1]], P[f[2]], P[b]) > 0;
  };
  auto edge = [&] (int f, int s) {
    return make_pair(hullInd[f][s], hullInd[f][(s + 1) % 3]);
  };
  auto glue = [&] (int af, int as, int bf, int bs) {
    other[af][as] = make_pair(bf, bs); other[bf][bs] = make_pair(af, as);
  };
  addFace(0, 1, 2); addFace(0, 2, 1); if (abv(1, 3)) swap(P[1], P[2]);
  for (int i = 0; i < 3; i++) glue(0, i, 1, 2 - i);
  for (int i = 3; i < n; i++) addEdge(abv(1, i), i);
  for (int i = 3; i < n; i++) {
    vector<int> rem; for (auto &&t : vis[i])
      if (active[t]) { active[t] = false; rem.push_back(t); }
    if (rem.empty()) continue;
    int st = -1; for (auto &&r : rem) for (int j = 0; j < 3; j++) {
      int o = other[r][j].first; if (active[o]) {
        int a, b; tie(a, b) = edge(r, j); addFace(a, b, i); st = a;
        int cur = int(rvis.size()) - 1; label[a] = cur; vector<int> tmp;
        set_union(rvis[r].begin(), rvis[r].end(),
                  rvis[o].begin(), rvis[o].end(), back_inserter(tmp));
        for (auto &&x : tmp) if (abv(cur, x)) addEdge(cur, x);
        glue(cur, 0, other[r][j].first, other[r][j].second);
      }
    }
    for (int x = st, y; ; x = y) {
      int lx = label[x]; glue(lx, 1, label[y = hullInd[lx][1]], 2);
      if (y == st) break;
    }
  }
  vector<vector<pt3>> hull;
  for (int i = 0; i < int(hullInd.size()); i++) if (active[i])
    hull.push_back(vector<pt3>{P[hullInd[i][0]], P[hullInd[i][1]],
                               P[hullInd[i][2]]});
  return hull;
}


text:
#pragma once
#include <bits/stdc++.h>
#include "Point3D.h"
using namespace std;

// Affine Transformations in 3D
// Functions:
//   prependMatrix(m2, b2): sets m = m2 * m and b = m2 * b + b2
//   transform(t): applies the AffineTransformation t to this
//   scale(p): scales the x, y, z coordinates by p.x, p.y, p.z respectively
//   translate(p): translates the point by p
//   rotate(axis, theta): rotates the point theta radians around the line from
//     the origin in the direction of axis
//   reflect(norm): reflects the point across the plane passing through the
//     origin with the normal vector norm
//   project(norm): projects the point across the plane passing through the
//     origin with the normal vector norm
//   applyTransform(p): applies the transformation to the point p
//   inverse(): returns the inverse of this transformation,
//     determinant of m must nonzero
// Time Complexity:
//   constructor, prependMatrix, transform, scale, translate, rotate, reflect,
//     project, applyTransform, inverse: O(1)
// Memory Complexity: O(1)
// Tested:
//   https://codeforces.com/problemsets/acmsguru/problem/99999/265
struct AffineTransform3D {
  static array<array<T, 3>, 3> inverse(const array<array<T, 3>, 3> &A) {
    T a = A[0][0], b = A[0][1], c = A[0][2];
    T d = A[1][0], e = A[1][1], f = A[1][2];
    T g = A[2][0], h = A[2][1], k = A[2][2];
    T det = a * (e * k - f * h) - b * (d * k - f * g) + c * (d * h - e * g);
    assert(!eq(det, 0)); array<array<T, 3>, 3> ret;
    ret[0][0] = (e * k - f * h) / det;
    ret[0][1] = -(b * k - c * h) / det;
    ret[0][2] = (b * f - c * e) / det;
    ret[1][0] = -(d * k - f * g) / det;
    ret[1][1] = (a * k - c * g) / det;
    ret[1][2] = -(a * f - c * d) / det;
    ret[2][0] = (d * h - e * g) / det;
    ret[2][1] = -(a * h - b * g) / det;
    ret[2][2] = (a * e - b * d) / det;
    return ret;
  }
  array<array<T, 3>, 3> m; array<T, 3> b;
  AffineTransform3D() {
    for (int i = 0; i < 3; i++) {
      b[i] = T(0); for (int j = 0; j < 3; j++) m[i][j] = i == j ? T(1) : T(0);
    }
  }
  void prependMatrix(const array<array<T, 3>, 3> &m2, const array<T, 3> &b2) {
    array<array<T, 3>, 3> resm; array<T, 3> resb;
    for (int i = 0; i < 3; i++) for (int j = 0; j < 3; j++) {
      resm[i][j] = T(0);
      for (int k = 0; k < 3; k++) resm[i][j] += m2[i][k] * m[k][j];
    }
    for (int i = 0; i < 3; i++) {
      resb[i] = b2[i]; for (int j = 0; j < 3; j++) resb[i] += m2[i][j] * b[j];
    }
    m = resm; b = resb;
  }
  void transform(const AffineTransform3D &t) { prependMatrix(t.m, t.b); }
  void scale(pt3 p) {
    prependMatrix({array<T, 3>{p.x, T(0), T(0)},
                              {T(0), p.y, T(0)},
                              {T(0), T(0), p.z}},
                  {T(0), T(0), T(0)});
  }
  void translate(pt3 p) {
    prependMatrix({array<T, 3>{T(1), T(0), T(0)},
                              {T(0), T(1), T(0)},
                              {T(0), T(0), T(1)}},
                  {p.x, p.y, p.z});
  }
  void rotate(pt3 axis, T theta) {
    axis = unit(axis);
    T x = axis.x, y = axis.y, z = axis.z, c = cos(theta), s = sin(theta);
    prependMatrix({array<T, 3>{x * x * (1 - c) + c,
                               x * y * (1 - c) - z * s,
                               x * z * (1 - c) + y * s},
                              {y * x * (1 - c) + z * s,
                               y * y * (1 - c) + c,
                               y * z * (1 - c) - x * s},
                              {z * x * (1 - c) - y * s,
                               z * y * (1 - c) + x * s,
                               z * z * (1 - c) + c}},
                  {T(0), T(0), T(0)});
  }
  void reflect(pt3 norm) {
    norm = unit(norm); T a = norm.x, b = norm.y, c = norm.z;
    prependMatrix({array<T, 3>{1 - 2 * a * a, -2 * a * b, -2 * a * c},
                              {-2 * b * a, 1 - 2 * b * b, -2 * b * c},
                              {-2 * c * a, -2 * c * b, 1 - 2 * c * c}},
                  {T(0), T(0), T(0)});
  }
  void project(pt3 norm) {
    norm = unit(norm); T a = norm.x, b = norm.y, c = norm.z;
    prependMatrix({array<T, 3>{b * b + c * c, -a * b, -a * c},
                              {-b * a, a * a + c * c, -b * c},
                              {-c * a, -c * b, a * a + b * b}},
                  {T(0), T(0), T(0)});
  }
  pt3 applyTransform(pt3 p) {
    return pt3(m[0][0] * p.x + m[0][1] * p.y + m[0][2] * p.z + b[0],
               m[1][0] * p.x + m[1][1] * p.y + m[1][2] * p.z + b[1],
               m[2][0] * p.x + m[2][1] * p.y + m[2][2] * p.z + b[2]);
  }
  AffineTransform3D inverse() const {
    AffineTransform3D ret;
    ret.prependMatrix({array<T, 3>{T(1), T(0), T(0)},
                                  {T(0), T(1), T(0)},
                                  {T(0), T(0), T(1)}},
                      {-b[0], -b[1], -b[2]});
    ret.prependMatrix(inverse(m), {T(0), T(0), T(0)}); return ret;
  }
};


text:
#pragma once
#include <bits/stdc++.h>
#include "Point3D.h"
#include "Line3D.h"
using namespace std;

// Functions for a 3D plane
struct Plane3D {
  pt3 n; T d;
  // ax + by + cz = d, above is ax + by + cz > d
  Plane3D(T a = 0, T b = 0, T c = 0, T d = 0) : n(a, b, c), d(d) {}
  // normal n, offset d
  Plane3D(pt3 n, T d) : n(n), d(d) {}
  // normal n, point p
  Plane3D(pt3 n, pt3 p) : n(n), d(n | p) {}
  // 3 non-collinear points p, q, r
  Plane3D(pt3 p, pt3 q, pt3 r) : Plane3D((q - p) * (r - p), p) {}
  T eval(pt3 p) const { return (n | p) - d; }
  // sign of isAbove, dist: 1 if above plane, 0 if on plane,
  //   -1 if below plane
  int isAbove(pt3 p) const { return sgn(eval(p)); }
  T dist(pt3 p) const { return eval(p) / abs(n); }
  T distSq(pt3 p) const { T e = eval(p); return e * e / norm(n); }
  Plane3D translate(pt3 p) const { return Plane3D(n, d + (n | p)); }
  Plane3D shiftUp(T e) const { return Plane3D(n, d + e * abs(n)); }
  pt3 proj(pt3 p) const { return p - n * eval(p) / norm(n); }
  pt3 refl(pt3 p) const { return p - n * T(2) * eval(p) / norm(n); }
  // returns 1 if the projections of a, b, c on this plane form a ccw turn
  //   looking from above, 0 if collinear, -1 if cw
  int ccwProj(pt3 a, pt3 b, pt3 c) const { return sgn((b - a) * (c - a) | n); }
  // returns a tuple (a, b, c) of 3 non-collinear points on the plane,
  //   guaranteed that ccwProj(a, b, c) == 1
  tuple<pt3, pt3, pt3> getPts() const {
    pt3 v = pt3(1, 0, 0); if (eq(abs(n * v), 0)) v = pt3(0, 1, 0);
    pt3 v1 = n * v, v2 = n * v1; pt3 a = proj(pt3(0, 0, 0));
    return make_tuple(a, a + v1, a + v2);
  }
};
Line3D perpThrough(Plane3D pi, pt3 o) {
  Line3D ret; ret.o = o; ret.d = pi.n; return ret;
}
Plane3D perpThrough(Line3D l, pt3 o) { return Plane3D(l.d, o); }

// Transforms points to a new coordinate system where the x and y axes are
//   on the plane, with the z axis being the normal vector (positive z is in
//   the direction of the normal vector)
// Z coordinate is guaranteed to be the distance to the plane (positive if
//   above plane, negative if below, 0 if on)
// Constructor Arguments:
//   pi: a plane
//   p, q, r: 3 non-collinear points
// Functions:
//   transform(p): transforms p into the new coordinate system
// Time Complexity:
//   constructor, transform: O(1)
// Memory Complexity: O(1)
// Tested:
//   https://dmoj.ca/problem/utso15p6
//   https://open.kattis.com/problems/starsinacan
struct CoordinateTransformation {
  pt3 o, dx, dy, dz;
  CoordinateTransformation(Plane3D pi) {
    pt3 p, q, r; tie(p, q, r) = pi.getPts(); o = p;
    dx = unit(q - p); dz = unit(dx * (r - p)); dy = dz * dx;
  }
  CoordinateTransformation(pt3 p, pt3 q, pt3 r) : o(p) {
    dx = unit(q - p); dz = unit(dx * (r - p)); dy = dz * dx;
  }
  pt3 transform(pt3 p) const {
    return pt3((p - o) | dx, (p - o) | dy, (p - o) | dz);
  }
};

// Intersection of plane and line
// Function Arguments:
//   pi: the plane
//   l: the line
//   res: a reference to the point to store the intersection if it exists
// Return Value: 0 if no intersection, 1 if point of intersection, 2 otherwise
// Time Complexity: O(1)
// Memory Complexity: O(1)
int planeLineIntersection(Plane3D pi, Line3D l, pt3 &res) {
  T a = pi.n | l.d; if (eq(norm(a), 0)) return pi.isAbove(l.o) == 0 ? 2 : 0;
  res = l.o - l.d * pi.eval(l.o) / a; return 1;
}

// Intersection of 2 planes
// Function Arguments:
//   pi1: the first plane
//   pi2: the second plane
//   res: a reference to the line to store the intersection if it exists
// Return Value: 0 if no intersection, 1 if line of intersection, 2 otherwise
// Time Complexity: O(1)
// Memory Complexity: O(1)
// Tested:
//   https://dmoj.ca/problem/utso15p6
int planePlaneIntersection(Plane3D pi1, Plane3D pi2, Line3D &res) {
  pt3 d = pi1.n * pi2.n; if (eq(norm(d), 0))
    return eq(abs(pi1.d / abs(pi1.n)), abs(pi2.d / abs(pi2.n))) ? 2 : 0;
  res.o = (pi2.n * pi1.d - pi1.n * pi2.d) * d / norm(d); res.d = d; return 1;
}


text:
#pragma once
#include <bits/stdc++.h>
#include "Point3D.h"
#include "Line3D.h"
#include "Plane3D.h"
using namespace std;

// Functions for a 3D sphere
struct Sphere3D {
  pt3 o; T r; Sphere3D(T r = 0) : o(0, 0), r(r) {}
  Sphere3D(pt3 o, T r) : o(o), r(r) {}
  // 1 if p is inside this sphere, 0 if on this sphere,
  //   -1 if outside this sphere
  int contains(pt3 p) const { return sgn(r * r - distSq(o, p)); }
  // 1 if s is strictly inside this sphere, 0 if inside and touching this
  //   sphere, -1 otherwise
  int contains(Sphere3D s) const {
    T dr = r - s.r; return lt(dr, 0) ? -1 : sgn(dr * dr - distSq(o, s.o));
  }
  // 1 if s is strictly outside this sphere, 0 if outside and touching this
  //   sphere, -1 otherwise
  int disjoint(Sphere3D s) const {
    T sr = r + s.r; return sgn(distSq(o, s.o) - sr * sr);
  }
  pt3 proj(pt3 p) const { return o + (p - o) * r / dist(o, p); }
  pt3 inv(pt3 p) const { return o + (p - o) * r * r / distSq(o, p); }
  // Shortest distance on the sphere between the projections of a and b onto
  //   this sphere
  T greatCircleDist(pt3 a, pt3 b) const { return r * ang(a, o, b); }
  // Is a-b a valid great circle segment (not on opposite sides)
  bool isGreatCircleSeg(pt3 a, pt3 b) const {
    assert(contains(a) == 0 && contains(b) == 0);
    a -= o; b -= o; return !eq(norm(a * b), 0) || lt(0, (a | b));
  }
  // Returns whether p is on the great circle segment a-b
  bool onSphSeg(pt3 p, pt3 a, pt3 b) const {
    assert(isGreatCircleSeg(a, b)); p -= o; a -= o; b -= o; pt3 n = a * b;
    if (eq(norm(n), 0)) return eq(norm(a * p), 0) && lt(0, (a | p));
    return eq((n | p), 0) && !lt((n | a * p), 0) && !lt(0, (n | b * p));
  }
  // Returns the points of intersection (or segment of intersection) between
  //   the great circle segments a-b and p-q
  vector<pt3> greatCircleSegIntersection(pt3 a, pt3 b, pt3 p, pt3 q) const {
    assert(isGreatCircleSeg(a, b) && isGreatCircleSeg(p, q));
    pt3 ab = (a - o) * (b - o), pq = (p - o) * (q - o);
    int oa = sgn(pq | (a - o)), ob = sgn(pq | (b - o));
    int op = sgn(ab | (p - o)), oq = sgn(ab | (q - o));
    if (oa != ob && op != oq && oa != op)
      return vector<pt3>{proj(o + ab * pq * op)};
    vector<pt3> ret; if (onSphSeg(p, a, b)) ret.push_back(p);
    if (onSphSeg(q, a, b)) ret.push_back(q);
    if (onSphSeg(a, p, q)) ret.push_back(a);
    if (onSphSeg(b, p, q)) ret.push_back(b);
    sort(ret.begin(), ret.end());
    ret.erase(unique(ret.begin(), ret.end()), ret.end()); return ret;
  }
  // Returns the angle between 3 points projected onto the sphere,
  //   positive angle if a-b-c forms a ccw turn, negative if cw, 0 if collinear
  T angSph(pt3 a, pt3 b, pt3 c) const {
    a -= o; b -= o; c -= o; T theta = ang(b * a, pt3(0, 0, 0), b * c);
    return (a * b | c) < 0 ? -theta : theta;
  }
  // Returns the surface area of a polygon projected onto the sphere,
  //   inside area if points are in ccw order, outside area if points are
  //   in cw order
  T surfaceAreaOnSph(const vector<pt3> &poly) const {
    int n = poly.size(); T PI = acos(T(-1)), a = -(n - 2) * PI;
    for (int i = 0; i < n; i++) {
      T ang = angSph(poly[i], poly[(i + 1) % n], poly[(i + 2) % n]);
      if (ang < 0) ang += 2 * PI;
      a += ang;
    }
    return r * r * a;
  }
};

// Determine the intersection of a sphere and a line
// Function Arguments:
//   s: the sphere
//   l: the line
// Return Value: the points of intersection (if any) of the sphere and
//   the line, guaranteed to be sorted based on projection on the line
// Time Complexity: O(1)
// Memory Complexity: O(1)
vector<pt3> sphereLineIntersection(Sphere3D s, Line3D l) {
  vector<pt3> ret; T h2 = s.r * s.r - l.distSq(s.o); pt3 p = l.proj(s.o);
  if (eq(h2, 0)) ret.push_back(p);
  else if (lt(0, h2)) {
    pt3 h = l.d * sqrt(h2) / abs(l.d);
    ret.push_back(p - h); ret.push_back(p + h);
  }
  return ret;
}

// Determine the intersection of a sphere and a plane
// Function Arguments:
//   s: the sphere
//   pi: the plane
//   res: a pair of pt3 and T, representing the center of the circle, and the
//     radius of the circle of intersection if it exists, guaranteed to be on
//     the plane pi
// Return Value: a boolean indicating whether an intersection exists or not
// Time Complexity: O(1)
// Memory Complexity: O(1)
bool spherePlaneIntersection(Sphere3D s, Plane3D pi, pair<pt3, T> &res) {
  T d2 = s.r * s.r - pi.distSq(s.o); if (lt(d2, 0)) return false;
  res.first = pi.proj(s.o); res.second = sqrt(max(d2, T(0))); return true;
}

// Determine the surface area and volume of the sphere cap above the
//   intersection of a sphere and a half-space defined by the space above
//   a plane (surface area does not include the base of the cap)
// Function Arguments:
//   s: the sphere
//   pi: the plane with the half-space defined as the space above the plane
// Return Value: a pair containing the surface area and the volume of the
//   sphere above the intersection of the sphere and the half-space
// Time Complexity: O(1)
// Memory Complexity: O(1)
pair<T, T> sphereHalfSpaceIntersectionSAV(Sphere3D s, Plane3D pi) {
  T d2 = s.r * s.r - pi.distSq(s.o);
  T h = lt(d2, 0) ? T(0) : s.r - abs(pi.dist(s.o));
  if (pi.isAbove(s.o) > 0) h = s.r * 2 - h;
  T PI = acos(T(-1));
  return make_pair(PI * 2 * s.r * h, PI * h * h / 3 * (3 * s.r - h));
}

// Determine the intersection of two spheres
// Function Arguments:
//   c1: the first circle
//   c2: the second circle
//   c: a tuple containing the plane the circle lies on (pointing away
//     from s1), the center of the circle, and the radius
// Return Value: 0 if no intersection, 2 if identical spheres, 1 otherwise
// Time Complexity: O(1)
// Memory Complexity: O(1)
int sphereSphereIntersection(Sphere3D s1, Sphere3D s2,
                             tuple<Plane3D, pt3, T> &c) {
  pt3 d = s2.o - s1.o; T d2 = norm(d);
  if (eq(d2, 0)) return eq(s1.r, s2.r) ? 2 : 0;
  T pd = (d2 + s1.r * s1.r - s2.r * s2.r) / 2, h2 = s1.r * s1.r - pd * pd / d2;
  if (lt(h2, 0)) return 0;
  pt3 o = s1.o + d * pd / d2;
  c = make_tuple(Plane3D(d, o), o, sqrt(max(h2, T(0)))); return 1;
}

// Determine the surface area and volume of the intersection of two spheres
// Function Arguments:
//   s1: the first sphere
//   s2: the second sphere
// Return Value: a pair containing the surface area and volume of the
//   intersection of the two spheres
// Time Complexity: O(1)
// Memory Complexity: O(1)
pair<T, T> sphereSphereIntersectionSAV(Sphere3D s1, Sphere3D s2) {
  pt3 d = s2.o - s1.o; T d2 = norm(d), dr = abs(s1.r - s2.r), PI = acos(T(-1));
  if (!lt(dr * dr, d2)) {
    T r = min(s1.r, s2.r);
    return make_pair(PI * 4 * r * r, PI * 4 * r * r * r / 3);
  }
  T sr = s1.r + s2.r; if (lt(sr * sr, d2)) return make_pair(T(0), T(0));
  T pd = (d2 + s1.r * s1.r - s2.r * s2.r) / 2;
  Plane3D pi = Plane3D(d, s1.o + d * pd / d2);
  pair<T, T> a = sphereHalfSpaceIntersectionSAV(s1, pi);
  pair<T, T> b = sphereHalfSpaceIntersectionSAV(s2, Plane3D(-pi.n, -pi.d));
  a.first += b.first; a.second += b.second; return a;
}


text:
#pragma once
#include <bits/stdc++.h>
#include "Point3D.h"
using namespace std;

// Functions for a 3D plane (represented in parametric form o + kd for real k)
struct Line3D {
  pt3 o, d;
  // points p and q
  Line3D(pt3 p = pt3(0, 0, 0), pt3 q = pt3(0, 0, 0)) : o(p), d(q - p) {}
  bool onLine(pt3 p) const { return eq(norm(d * (p - o)), 0); }
  T distSq(pt3 p) const { return norm(d * (p - o)) / norm(d); }
  T dist(pt3 p) const { return sqrt(distSq(p)); }
  Line3D translate(pt3 p) const { return Line3D(o + p, d); }
  pt3 proj(pt3 p) const { return o + d * (d | (p - o)) / norm(d); }
  pt3 refl(pt3 p) const { return proj(p) * T(2) - p; }
  // compares points by projection (3 way comparison)
  int cmpProj(pt3 p, pt3 q) const { return sgn((d | p) - (d | q)); }
};

// Closest distance between 2 lines
// Function Arguments:
//   l1: the first line
//   l2: the second line
// Return Value: the closest distance between 2 lines
// Time Complexity: O(1)
// Memory Complexity: O(1)
T lineLineDist(Line3D l1, Line3D l2) {
  pt3 n = l1.d * l2.d; if (eq(norm(n), 0)) return l1.dist(l2.o);
  return abs((l2.o - l1.o) | n) / abs(n);
}

// Closest point on the line l1 to the line l2
// If l1 and l2 are parallel, it returns the projection of (0, 0, 0) on l1
// Function Arguments:
//   l1: the first line
//   l2: the second line
// Return Value: the closest point on the line l1 to the line l2
// Time Complexity: O(1)
// Memory Complexity: O(1)
pt3 closestOnL1ToL2(Line3D l1, Line3D l2) {
  pt3 n = l1.d * l2.d; if (eq(norm(n), 0)) return l1.proj(pt3(0, 0, 0));
  pt3 n2 = l2.d * n; return l1.o + l1.d * ((l2.o - l1.o) | n2) / (l1.d | n2);
}

// Intersection of 2 lines
// Function Arguments:
//   l1: the first line
//   l2: the second line
//   res: a reference to the point to store the intersection if it exists
// Return Value: 0 if no intersection, 1 if point of intersection, 2 otherwise
// Time Complexity: O(1)
// Memory Complexity: O(1)
int lineLineIntersection(Line3D l1, Line3D l2, pt3 &res) {
  pt3 n = l1.d * l2.d; if (eq(norm(n), 0)) return eq(l1.dist(l2.o), 0) ? 2 : 0;
  res = closestOnL1ToL2(l1, l2); return 1;
}

// Determines if a point is on a line segment
// Function Arguments:
//   p: the point to check if on the line segment
//   a: one endpoint of the line segment
//   b: the other endpoint of the line segment
// Return Value: true if p is on the line segment a-b, false otherwise
// Time Complexity: O(1)
// Memory Complexity: O(1)
bool onSeg(pt3 p, pt3 a, pt3 b) {
  if (a == b) return p == a;
  Line3D l(a, b);
  return l.onLine(p) && l.cmpProj(a, p) <= 0 && l.cmpProj(p, b) <= 0;
}

// Determine the intersection of two line segments
// Function Arguments:
//   a: one endpoint of the first line segment
//   b: the other endpoint of the first line segment
//   p: one endpoint of the second line segment
//   q: the other endpoint of the second line segment
// Return Value: if the line segments do not intersect, an empty vector
//   of points; if the line segments intersect at a point, a vector containing
//   the point of intersection; if the line segments have a line segment of
//   intersection, a vector containing the two endpoints of the
//   line segment intersection (it can return more if there is precision error)
vector<pt3> segSegIntersection(pt3 a, pt3 b, pt3 p, pt3 q) {
  vector<pt3> ret; if (a == b) {
    if (onSeg(a, p, q)) ret.push_back(a);
  } else if (p == q) {
    if (onSeg(p, a, b)) ret.push_back(p);
  } else {
    pt3 inter; Line3D l1(a, b), l2(p, q);
    int cnt = lineLineIntersection(l1, l2, inter);
    if (cnt == 1) ret.push_back(inter);
    else if (cnt == 2) {
      if (onSeg(p, a, b)) ret.push_back(p);
      if (onSeg(q, a, b)) ret.push_back(q);
      if (onSeg(a, p, q)) ret.push_back(a);
      if (onSeg(b, p, q)) ret.push_back(b);
    }
  }
  sort(ret.begin(), ret.end());
  ret.erase(unique(ret.begin(), ret.end()), ret.end()); return ret;
}

// Finds the closest point on a line segment to another point
// Function Arguments
//   p: the reference point
//   a: one endpoint of the line segment
//   b: the other endpoint of the line segment
// Return Value: the closest point to p on the line segment a-b
// Time Complexity: O(1)
// Memory Complexity: O(1)
pt3 closestPtOnSeg(pt3 p, pt3 a, pt3 b) {
  if (a != b) {
    Line3D l(a, b);
    if (l.cmpProj(a, p) < 0 && l.cmpProj(p, b) < 0) return l.proj(p);
  }
  return lt(dist(p, a), dist(p, b)) ? a : b;
}

// Finds the distance to the closest point on a line segment to another point
// Function Arguments
//   p: the reference point
//   a: one endpoint of the line segment
//   b: the other endpoint of the line segment
// Return Value: the distance to the closest point to p on the line segment a-b
// Time Complexity: O(1)
// Memory Complexity: O(1)
T ptSegDist(pt3 p, pt3 a, pt3 b) {
  if (a != b) {
    Line3D l(a, b);
    if (l.cmpProj(a, p) < 0 && l.cmpProj(p, b) < 0) return l.dist(p);
  }
  return min(dist(p, a), dist(p, b));
}

// Finds the closest distance between two line segments
// Function Arguments
//   a: one endpoint of the first line segment
//   b: the other endpoint of the first line segment
//   p: one endpoint of the second line segment
//   q: the other endpoint of the second line segment
// Return Value: the closest distance between the two line
// Time Complexity: O(1)
// Memory Complexity: O(1)
T segSegDist(pt3 a, pt3 b, pt3 p, pt3 q) {
  return !segSegIntersection(a, b, p, q).empty()
      ? 0
      : min({ptSegDist(p, a, b), ptSegDist(q, a, b),
             ptSegDist(a, p, q), ptSegDist(b, p, q)});
}


text:
#pragma once
#include <bits/stdc++.h>
#include "../../utils/Random.h"
#include "Point.h"
#include "Circle.h"
using namespace std;

// Determines the minimum enclosing circle of a set of points
// Function Arguments:
//   P: the points
// Return Value: the minimum enclosing circle of P
// In practice, has a large constant
// Time Complexity: O(N) expected
// Memory Complexity: O(N)
// Tested:
//   https://www.spoj.com/problems/QCJ4/
//   https://open.kattis.com/problems/starsinacan
Circle minEnclosingCircle(vector<pt> P) {
  shuffle(P.begin(), P.end(), rng); Circle c(P[0], 0);
  for (int i = 0; i < int(P.size()); i++) if (lt(c.r, dist(P[i], c.o))) {
    c = Circle(P[i], 0);
    for (int j = 0; j < i; j++) if (lt(c.r, dist(P[j], c.o))) {
      pt p = (P[i] + P[j]) / T(2); c = Circle(p, dist(P[i], p));
      for (int k = 0; k < j; k++) if (lt(c.r, dist(P[k], c.o)))
        c = circumcircle(P[i], P[j], P[k]);
    }
  }
  return c;
}


text:
#pragma once
#include <bits/stdc++.h>
#include "Point.h"
#include "Angle.h"
using namespace std;

// Computes the minimum area from 3 distinct points out of a set of N points
// Angle::pivot is modified
// Constructor Arguments:
//   P: the vector of points
// Fields:
//   PA: the first point in the minimum area triangle
//   PB: the second point in the minimum area triangle
//   PC: the third point in the minimum area triangle
//   minArea2: twice the area of the minimum area triangle
//   hull: the points in the convex hull of P
// In practice, has a moderate constant
// Time Complexity:
//   constructor: O(N^2 log N)
// Memory Complexity: O(N)
// Tested:
//   Fuzz Tested
struct MinTriangleArea {
  pt PA, PB, PC; T minArea2;
  MinTriangleArea(const vector<pt> &P) : minArea2(numeric_limits<T>::max()) {
    int n = P.size(); if (n < 3) return;
    vector<Angle> A(n); for (int i = 0; i < n; i++) A[i] = Angle(P[i]);
    for (int i = 0; i < n; i++) {
      Angle::setPivot(P[i]); sort(A.begin(), A.end());
      for (int j = 0; j < n - 1; j++) {
        pt b = A[j].p, c = A[(j + 1) % (n - 1)].p;
        T a2 = area2(P[i], b, c); if (a2 < 0) { swap(b, c); a2 = -a2; }
        if (a2 < minArea2) { PA = P[i]; PB = b; PC = c; minArea2 = a2; }
      }
    }
  }
};


text:
#pragma once
#include <bits/stdc++.h>
#include "Point.h"
using namespace std;

// Computes the Delaunay Triangulation of a set of distinct points
// If all points are collinear, there is no triangulation
// If there are 4 or more points on the same circle, the triangulation is
//   ambiguous, otherwise it is unique
// Each circumcircle does not completely contain any of the input points
// Constructor Arguments:
//   P: the distinct points
// Fields:
//   tri: a vector of arrays of 3 points, representing a triangle in
//     the triangulation, points given in ccw order
// In practice, has a moderate constant
// Time Complexity:
//   construction: O(N log N)
// Memory Complexity: O(N)
// Tested:
//   https://dmoj.ca/problem/cco08p6
//   https://codeforces.com/problemsets/acmsguru/problem/99999/383
struct DelaunayTriangulation {
  static constexpr const pt inf = pt(numeric_limits<T>::max(),
                                     numeric_limits<T>::max());
  struct Quad {
    static vector<Quad> V;
    int rot, o; pt p; Quad(int rot) : rot(rot), o(-1), p(inf) {}
    int &r() { return V[rot].rot; }
    pt &f() { return V[r()].p; }
    int prv() { return V[V[rot].o].rot; }
    int nxt() { return V[r()].prv(); }
  };
  int head; vector<Quad> &V = Quad::V; vector<array<pt, 3>> tri;
  int makeQuad(int rot) { V.emplace_back(rot); return int(V.size()) - 1; }
  bool circ(pt p, pt a, pt b, pt c) {
    T p2 = norm(p), A = norm(a) - p2, B = norm(b) - p2, C = norm(c) - p2;
    return lt(0, area2(p, a, b) * C + area2(p, b, c) * A + area2(p, c, a) * B);
  }
  int makeEdge(pt a, pt b) {
    int r = ~head ? head : makeQuad(makeQuad(makeQuad(makeQuad(-1))));
    head = V[r].o; V[V[r].r()].r() = r; for (int i = 0; i < 4; i++) {
      r = V[r].rot; V[r].o = i & 1 ? r : V[r].r();
    }
    V[r].p = a; V[r].f() = b; return r;
  }
  void splice(int a, int b) {
    swap(V[V[V[a].o].rot].o, V[V[V[b].o].rot].o); swap(V[a].o, V[b].o);
  }
  int connect(int a, int b) {
    int q = makeEdge(V[a].f(), V[b].p);
    splice(q, V[a].nxt()); splice(V[q].r(), b); return q;
  }
  bool valid(int a, int base) {
    return ccw(V[a].f(), V[base].f(), V[base].p) > 0;
  }
  pair<int, int> rec(const vector<pt> &P, int lo, int hi) {
    int k = hi - lo + 1; if (k <= 3) {
      int a = makeEdge(P[lo], P[lo + 1]), b = makeEdge(P[lo + 1], P[hi]);
      if (k == 2) return make_pair(a, V[a].r());
      splice(V[a].r(), b); int side = sgn(ccw(P[lo], P[lo + 1], P[hi]));
      int c = side ? connect(b, a) : -1;
      return make_pair(side < 0 ? V[c].r() : a, side < 0 ? c : V[b].r());
    }
    int a, b, ra, rb, mid = lo + (hi - lo) / 2;
    tie(ra, a) = rec(P, lo, mid); tie(b, rb) = rec(P, mid + 1, hi);
    while ((ccw(V[b].p, V[a].f(), V[a].p) < 0 && (a = V[a].nxt()))
        || (ccw(V[a].p, V[b].f(), V[b].p) > 0 && (b = V[V[b].r()].o)));
    int base = connect(V[b].r(), a); if (V[a].p == V[ra].p) ra = V[base].r();
    if (V[b].p == V[rb].p) rb = base;
    while (true) {
      int l = V[V[base].r()].o; if (valid(l, base))
        while (circ(V[V[l].o].f(), V[base].f(), V[base].p, V[l].f())) {
          int t = V[l].o; splice(l, V[l].prv());
          splice(V[l].r(), V[V[l].r()].prv()); V[l].o = head; head = l; l = t;
        }
      int r = V[base].prv(); if (valid(r, base))
        while (circ(V[V[r].prv()].f(), V[base].f(), V[base].p, V[r].f())) {
          int t = V[r].prv(); splice(r, V[r].prv());
          splice(V[r].r(), V[V[r].r()].prv()); V[r].o = head; head = r; r = t;
        }
      if (!valid(l, base) && !valid(r, base)) break;
      if (!valid(l, base) || (valid(r, base)
          && circ(V[r].f(), V[r].p, V[l].f(), V[l].p)))
        base = connect(r, V[base].r());
      else base = connect(V[base].r(), V[l].r());
    }
    return make_pair(ra, rb);
  }
  DelaunayTriangulation(vector<pt> P) : head(-1) {
    sort(P.begin(), P.end()); assert(unique(P.begin(), P.end()) == P.end());
    V.reserve(P.size() * 16); if (int(P.size()) < 2) return;
    int e = rec(P, 0, int(P.size()) - 1).first, qi = 0, ind = 0;
    vector<bool> vis(V.size(), false); vector<int> q{e}; q.reserve(V.size());
    while (ccw(V[V[e].o].f(), V[e].f(), V[e].p) < 0) e = V[e].o;
    auto add = [&] {
      int c = e; do {
        if (ind % 3 == 0) tri.emplace_back();
        tri.back()[ind++ % 3] = V[c].p;
        vis[c] = true; q.push_back(V[c].r()); c = V[c].nxt();
      } while (c != e);
    };
    add(); tri.clear(); ind = 0; tri.reserve(V.size() / 3 + 1);
    while (qi < int(q.size())) if (!vis[e = q[qi++]]) add();
    assert(ind % 3 == 0); V = vector<Quad>();
  }
};

vector<DelaunayTriangulation::Quad> DelaunayTriangulation::Quad::V
    = vector<DelaunayTriangulation::Quad>();


text:
#pragma once
#include <bits/stdc++.h>
#include "Point.h"
#include "Line.h"
using namespace std;

// Helper struct for bentleyOttmann
struct Seg {
  static T X; mutable pt p, q; mutable int i; bool isQuery;
  pair<T, T> y() const {
    pt v = q - p; return make_pair(cross(v, p) + v.y * X, v.x);
  }
  Seg(pt p, pt q, int i) : p(p), q(q), i(i), isQuery(false) {
    assert(!eq(p.x, q.x)); if (p > q) swap(this->p, this->q);
  }
  Seg(T Y) : p(X, Y), q(p), i(-1), isQuery(true) {}
  bool operator < (const Seg &s) const {
    if (isQuery) {
      pair<T, T> p2 = s.y(); return lt(p.y * p2.second, p2.first);
    } else if (s.isQuery) {
      pair<T, T> p1 = y(); return lt(p1.first, s.p.y * p1.second);
    }
    pair<T, T> p1 = y(), p2 = s.y();
    return lt(p1.first * p2.second, p2.first * p1.second);
  }
};

T Seg::X = 0;

// Helper struct for bentleyOttmann
struct SegEvent {
  pt p; int type, i, j;
  SegEvent(pt p, int type, int i, int j) : p(p), type(type), i(i), j(j) {}
  bool operator > (const SegEvent &e) const {
    return eq(e.p.x, p.x) ? (eq(e.type, type) ? lt(e.p.y, p.y)
                                              : lt(e.type, type))
                          : lt(e.p.x, p.x);
  }
};

// Helper struct for bentleyOttmann
struct VerticalSeg {
  T lo, hi; int i; VerticalSeg(T lo, T hi, int i) : lo(lo), hi(hi), i(i) {}
  bool operator < (const VerticalSeg &vs) const { return lt(hi, vs.hi); }
};

// Runs a callback on all pairs of intersecting line segments (proper or not)
// Template Arguments:
//   F: the type of the function f
// Function Arguments:
//   segs: a vector of pairs of the endpoints of the line segments
//   f(i, j): the function to run a callback on for each pair of intersecting
//     line segments, where i and j are the indices of the intersecting
//     line segments; each unordered pair (i, j) is passed at most once
// In practice, has a moderate constant
// Time Complexity: O((N + K) log N) for N line segments and K intersections
// Memory Complexity: O(N + K)
// Tested:
//   Fuzz Tested
//   https://open.kattis.com/problems/polygon
//   https://open.kattis.com/problems/scholarslawn
template <class F> void bentleyOttmann(const vector<pair<pt, pt>> &segs, F f) {
  std::priority_queue<SegEvent, vector<SegEvent>, greater<SegEvent>> events;
  int n = segs.size(); for (int i = 0; i < n; i++) {
    pt a, b; tie(a, b) = segs[i]; if (a > b) swap(a, b);
    int isEq = eq(a.x, b.x);
    events.emplace(a, isEq, i, i); events.emplace(b, 4 - isEq, i, i);
  }
  multiset<Seg> active; using iter = decltype(active)::iterator;
  vector<iter> ptr(n, active.end()); set<pair<int, int>> seen;
  auto checkSeen = [&] (int i, int j) {
    if (seen.count(make_pair(i, j)) || seen.count(make_pair(j, i)))
      return true;
    seen.emplace(i, j); return false;
  };
  auto checkInter = [&] (iter a, iter b) {
    vector<pt> p = segSegIntersection(a->p, a->q, b->p, b->q);
    if (int(p.size()) == 1 && cross(a->q - a->p, b->q - b->p) < 0) {
      if (!checkSeen(a->i, b->i)) f(a->i, b->i);
      events.emplace(p[0], 2, a->i, b->i);
    }
  };
  auto checkRange = [&] (int i, T lo, T hi) {
    auto it = active.lower_bound(Seg(lo));
    auto check = [&] {
      pair<T, T> y = it->y(); return !lt(hi * y.second, y.first);
    };
    for (; it != active.end() && check(); it++) {
      pt a, b; tie(a, b) = segs[i];
      if (segSegIntersects(a, b, it->p, it->q) && !checkSeen(i, it->i))
        f(i, it->i);
    }
  };
  multiset<VerticalSeg> vertical; while (!events.empty()) {
    SegEvent e = events.top(); events.pop(); Seg::X = e.p.x; if (e.type == 0) {
      checkRange(e.i, e.p.y, e.p.y);
      ptr[e.i] = active.emplace(segs[e.i].first, segs[e.i].second, e.i);
      if (ptr[e.i] != active.begin()) checkInter(prev(ptr[e.i]), ptr[e.i]);
      if (next(ptr[e.i]) != active.end()) checkInter(ptr[e.i], next(ptr[e.i]));
    } else if (e.type == 1) {
      T lo = segs[e.i].first.y, hi = segs[e.i].second.y;
      if (lo > hi) swap(lo, hi);
      checkRange(e.i, lo, hi);
      auto it = vertical.lower_bound(VerticalSeg(lo, lo, e.i));
      for (; it != vertical.end(); it++) if (!checkSeen(e.i, it->i))
        f(e.i, it->i);
      vertical.emplace(lo, hi, e.i);
    } else if (e.type == 2) {
      if (next(ptr[e.i]) != ptr[e.j]) continue;
      if (cross(ptr[e.i]->q - ptr[e.i]->p, ptr[e.j]->q - ptr[e.j]->p) < 0) {
        swap(ptr[e.i], ptr[e.j]); swap(ptr[e.i]->i, ptr[e.j]->i);
        swap(ptr[e.i]->p, ptr[e.j]->p); swap(ptr[e.i]->q, ptr[e.j]->q);
        if (ptr[e.j] != active.begin()) checkInter(prev(ptr[e.j]), ptr[e.j]);
        if (next(ptr[e.i]) != active.end())
          checkInter(ptr[e.i], next(ptr[e.i]));
      }
    } else if (e.type == 4) {
      if (ptr[e.i] != active.begin() && next(ptr[e.i]) != active.end())
        checkInter(prev(ptr[e.i]), next(ptr[e.i]));
      active.erase(ptr[e.i]); ptr[e.i] = active.end();
      checkRange(e.i, e.p.y, e.p.y);
    } else vertical.clear();
  }
}

// Transforms a set of line segments into another set of line segments
//   such that the union of the line segments remains the same, and any
//   intersection between any pair of segments is an endpoint of both segments
// Function Arguments:
//   segs: a vector of pairs of the endpoints of the line segments
// Return Value: the resulting vector of pairs of the endpoints of the
//   line segments after the transformation
// In practice, has a moderate constant
// Time Complexity: O((N + K) log N) for N line segments and K intersections
// Memory Complexity: O(N + K)
// Tested:
//   https://open.kattis.com/problems/scholarslawn
vector<pair<pt, pt>> segArrangement(const vector<pair<pt, pt>> &segs) {
  int n = segs.size(); vector<vector<pt>> ptsOnSegs(n); set<pair<pt, pt>> ret;
  bentleyOttmann(segs, [&] (int i, int j) {
    for (auto &&p : segSegIntersection(segs[i].first, segs[i].second,
                                       segs[j].first, segs[j].second)) {
      ptsOnSegs[i].push_back(p);
      ptsOnSegs[j].push_back(p);
    }
  });
  for (int i = 0; i < n; i++) {
    ptsOnSegs[i].push_back(segs[i].first);
    ptsOnSegs[i].push_back(segs[i].second);
    sort(ptsOnSegs[i].begin(), ptsOnSegs[i].end());
    ptsOnSegs[i].erase(unique(ptsOnSegs[i].begin(), ptsOnSegs[i].end()),
                       ptsOnSegs[i].end());
    for (int j = 0; j < int(ptsOnSegs[i].size()) - 1; j++) {
      pt a = ptsOnSegs[i][j], b = ptsOnSegs[i][j + 1];
      if (!ret.count(make_pair(a, b)) && !ret.count(make_pair(b, a)))
        ret.emplace(a, b);
    }
  }
  return vector<pair<pt, pt>>(ret.begin(), ret.end());
}


text:
#pragma once
#include <bits/stdc++.h>
using namespace std;

// Supports adding lines in the form f(x) = mx + b, or line segments
//   in the form of f(x) = mx + b over l <= x <= r, and finding
//   the maximum value of f(x) where all possible value of x are
//   known beforehand
// Template Arguments:
//   T: the type of the slope (m) and intercept (b) of the line, as well as
//     the type of the function argument (x)
//   Cmp: the comparator to compare two f(x) values,
//       convention is same as std::priority_queue in STL
//     Required Functions:
//       operator (a, b): returns true if and only if a compares less than b
// Constructor Arguments:
//   xs: a vector of type T of the possible x values for each query
//   cmp: an instance of the Cmp struct
//   inf: a value for positive infinity, must be negatable
// Functions:
//   addLine(m, b): adds a line in the form f(x) = mx + b to the set of lines
//   addLine(m, b, l, r): adds a line segment in the form f(x) = mx + b
//     where l <= x <= r, to the set of lines
//   getMax(x): finds the maximum value of f(x) (based on the comparator)
//     for all inserted lines, x must be in the vector xs
//   clear(): removes all lines from the seg
// In practice, has a moderate constant, performance compared to
//   IncrementalConvexHullTrick (which uses multiset) and
//   IncrementalConvexHullTrickSqrtBuffer can vary, faster than
//   SparseLiChaoTree
// Time Complexity:
//   constructor: O(1)
//   addLine, getMax: O(log N) for the range [0, N)
//   addLineSegment: O((log N) ^ 2) for the range [0, N)
//   clear: O(N) for the range [0, N)
// Memory Complexity: O(N) for the range [0, N)
// Tested:
//   https://judge.yosupo.jp/problem/line_add_get_min
//   https://judge.yosupo.jp/problem/segment_add_get_min
//   https://open.kattis.com/problems/longestlife
//   https://www.spoj.com/problems/CHTPRAC/
template <class T, class Cmp = less<T>> struct LiChaoTree {
  using Line = pair<T, T>; int N; Cmp cmp; T INF; vector<Line> TR; vector<T> X;
  T eval(Line l, int i) const { return l.first * X[i] + l.second; }
  bool majorize(Line a, Line b, int l, int r) {
    return !cmp(eval(a, l), eval(b, l)) && !cmp(eval(a, r), eval(b, r));
  }
  int cInd(T x) const {
    return lower_bound(X.begin(), X.end(), x) - X.begin();
  }
  int fInd(T x) const {
    return upper_bound(X.begin(), X.end(), x) - X.begin() - 1;
  }
  LiChaoTree(const vector<T> &xs, Cmp cmp = Cmp(),
             T inf = numeric_limits<T>::max())
      : cmp(cmp), INF(min(inf, -inf, cmp)), X(xs) {
    sort(X.begin(), X.end()); X.erase(unique(X.begin(), X.end()), X.end());
    N = X.size(); TR.assign(N == 0 ? 0 : 1 << __lg(N * 4 - 1), Line(T(), INF));
  }
  void addLine(int k, int tl, int tr, Line line) {
    if (majorize(line, TR[k], tl, tr)) swap(line, TR[k]);
    if (majorize(TR[k], line, tl, tr)) return;
    if (cmp(eval(TR[k], tl), eval(line, tl))) swap(line, TR[k]);
    int m = tl + (tr - tl) / 2; if (!cmp(eval(line, m), eval(TR[k], m))) {
      swap(line, TR[k]); addLine(k * 2, tl, m, line);
    } else addLine(k * 2 + 1, m + 1, tr, line);
  }
  void addLineSegment(int k, int tl, int tr, int l, int r, Line line) {
    if (r < tl || tr < l) return;
    if (l <= tl && tr <= r) { addLine(k, tl, tr, line); return; }
    int m = tl + (tr - tl) / 2; addLineSegment(k * 2, tl, m, l, r, line);
    addLineSegment(k * 2 + 1, m + 1, tr, l, r, line);
  }
  T getMax(int k, int tl, int tr, int i) const {
    T ret = eval(TR[k], i); if (tl == tr) return ret;
    int m = tl + (tr - tl) / 2;
    if (i <= m) return max(ret, getMax(k * 2, tl, m, i), cmp);
    else return max(ret, getMax(k * 2 + 1, m + 1, tr, i), cmp);
  }
  void addLine(T m, T b) { addLine(1, 0, N - 1, Line(m, b)); }
  void addLineSegment(T m, T b, T l, T r) {
    int li = cInd(l), ri = fInd(r);
    if (li <= ri) addLineSegment(1, 0, N - 1, li, ri, Line(m, b));
  }
  T getMax(T x) const { return getMax(1, 0, N - 1, cInd(x)); }
  void clear() { fill(TR.begin(), TR.end(), Line(T(), INF)); }
};