text:
#include <bits/stdc++.h>
#include "../../../../../Content/C++/datastructures/trees/fenwicktrees/FenwickTree.h"
#include "../../../../../Content/C++/datastructures/trees/fenwicktrees/FenwickTree1D.h"
#include "../../../../../Content/C++/search/BinarySearch.h"
using namespace std;

void test1() {
  const auto start_time = chrono::system_clock::now();
  mt19937_64 rng(0);
  int N = 1e7;
  vector<long long> A(N);
  for (auto &&ai : A) ai = rng() % int(1e6) + 1;
  FenwickTree1D<long long> FT(A);
  int Q = 1e7;
  vector<long long> ans;
  for (int i = 0; i < Q; i++) {
    int t = rng() % 2;
    if (t == 0) {
      int i = rng() % N;
      long long v = rng() % int(1e6) + 1;
      FT.update(i, v);
    } else {
      int l = rng() % N, r = rng() % N;
      if (l > r) swap(l, r);
      ans.push_back(FT.query(l, r));
    }
  }
  const auto end_time = chrono::system_clock::now();
  double sec = ((end_time - start_time).count() / double(chrono::system_clock::period::den));
  cout << "Subtest 1 (1D vector Constructor) Passed" << endl;
  cout << "  N: " << N << endl;
  cout << "  Q: " << Q << endl;
  cout << "  Time: " << fixed << setprecision(3) << sec << "s" << endl;
  long long checkSum = 0;
  for (auto &&a : ans) checkSum = 31 * checkSum + a;
  cout << "  Checksum: " << checkSum << endl;
}

void test2() {
  const auto start_time = chrono::system_clock::now();
  mt19937_64 rng(0);
  int N = 1e7;
  vector<long long> A(N);
  for (auto &&ai : A) ai = rng() % int(1e6) + 1;
  FenwickTree1D<long long> FT(N);
  for (int i = 0; i < N; i++) FT.update(i, A[i]);
  int Q = 1e7;
  vector<long long> ans;
  for (int i = 0; i < Q; i++) {
    int t = rng() % 2;
    if (t == 0) {
      int i = rng() % N;
      long long v = rng() % int(1e6) + 1;
      FT.update(i, v);
    } else {
      int l = rng() % N, r = rng() % N;
      if (l > r) swap(l, r);
      ans.push_back(FT.query(l, r));
    }
  }
  const auto end_time = chrono::system_clock::now();
  double sec = ((end_time - start_time).count() / double(chrono::system_clock::period::den));
  cout << "Subtest 2 (1D default value constructor) Passed" << endl;
  cout << "  N: " << N << endl;
  cout << "  Q: " << Q << endl;
  cout << "  Time: " << fixed << setprecision(3) << sec << "s" << endl;
  long long checkSum = 0;
  for (auto &&a : ans) checkSum = 31 * checkSum + a;
  cout << "  Checksum: " << checkSum << endl;
}

void test3() {
  const auto start_time = chrono::system_clock::now();
  mt19937_64 rng(0);
  int N = 1e7;
  vector<long long> A(N);
  for (auto &&ai : A) ai = rng() % int(1e6) + 1;
  FenwickTree<1, long long> FT(N);
  for (int i = 0; i < N; i++) FT.update(i, A[i]);
  int Q = 1e7;
  vector<long long> ans;
  for (int i = 0; i < Q; i++) {
    int t = rng() % 2;
    if (t == 0) {
      int i = rng() % N;
      long long v = rng() % int(1e6) + 1;
      FT.update(i, v);
    } else {
      int l = rng() % N, r = rng() % N;
      if (l > r) swap(l, r);
      ans.push_back(FT.query(l, r));
    }
  }
  const auto end_time = chrono::system_clock::now();
  double sec = ((end_time - start_time).count() / double(chrono::system_clock::period::den));
  cout << "Subtest 3 (ND) Passed" << endl;
  cout << "  N: " << N << endl;
  cout << "  Q: " << Q << endl;
  cout << "  Time: " << fixed << setprecision(3) << sec << "s" << endl;
  long long checkSum = 0;
  for (auto &&a : ans) checkSum = 31 * checkSum + a;
  cout << "  Checksum: " << checkSum << endl;
}

void test4() {
  const auto start_time = chrono::system_clock::now();
  mt19937_64 rng(0);
  int N = 1e7;
  vector<long long> A(N);
  for (auto &&ai : A) ai = rng() % int(1e6) + 1;
  FenwickTree1D<long long> FT(A);
  int Q = 1e7;
  vector<long long> ans;
  for (int i = 0; i < Q; i++) {
    int t = rng() % 4;
    if (t == 0) {
      int i = rng() % N;
      long long v = rng() % int(1e6) + 1;
      FT.update(i, v);
    } else if (t == 1) {
      int l = rng() % N, r = rng() % N;
      if (l > r) swap(l, r);
      ans.push_back(FT.query(l, r));
    } else if (t == 2) {
      long long TOT = FT.query(N - 1) * 2 + 1;
      long long v = rng() % TOT;
      ans.push_back(FT.bsearch(v, less<long long>()));
    } else {
      long long TOT = FT.query(N - 1) * 2 + 1;
      long long v = rng() % TOT;
      ans.push_back(FT.bsearch(v, less_equal<long long>()));
    }
  }
  const auto end_time = chrono::system_clock::now();
  double sec = ((end_time - start_time).count() / double(chrono::system_clock::period::den));
  cout << "Subtest 4 (lower_bound, upper_bound) Passed" << endl;
  cout << "  N: " << N << endl;
  cout << "  Q: " << Q << endl;
  cout << "  Time: " << fixed << setprecision(3) << sec << "s" << endl;
  long long checkSum = 0;
  for (auto &&a : ans) checkSum = 31 * checkSum + a;
  cout << "  Checksum: " << checkSum << endl;
}

void test5() {
  const auto start_time = chrono::system_clock::now();
  mt19937_64 rng(0);
  int N = 1e7;
  vector<long long> A(N);
  for (auto &&ai : A) ai = rng() % int(1e6) + 1;
  FenwickTree1D<long long> FT(A);
  int Q = 1e7;
  vector<long long> ans;
  for (int i = 0; i < Q; i++) {
    int t = rng() % 4;
    if (t == 0) {
      int i = rng() % N;
      long long v = rng() % int(1e6) + 1;
      FT.update(i, v);
    } else if (t == 1) {
      int l = rng() % N, r = rng() % N;
      if (l > r) swap(l, r);
      ans.push_back(FT.query(l, r));
    } else if (t == 2) {
      long long TOT = FT.query(N - 1) * 2 + 1;
      long long v = rng() % TOT;
      ans.push_back(bsearch<FIRST>(0, N, [&] (int k) {
        return FT.query(k) >= v;
      }));
    } else {
      long long TOT = FT.query(N - 1) * 2 + 1;
      long long v = rng() % TOT;
      ans.push_back(bsearch<FIRST>(0, N, [&] (int k) {
        return FT.query(k) > v;
      }));
    }
  }
  const auto end_time = chrono::system_clock::now();
  double sec = ((end_time - start_time).count() / double(chrono::system_clock::period::den));
  cout << "Subtest 5 (bsearch with query) Passed" << endl;
  cout << "  N: " << N << endl;
  cout << "  Q: " << Q << endl;
  cout << "  Time: " << fixed << setprecision(3) << sec << "s" << endl;
  long long checkSum = 0;
  for (auto &&a : ans) checkSum = 31 * checkSum + a;
  cout << "  Checksum: " << checkSum << endl;
}

int main() {
  test1();
  test2();
  test3();
  test4();
  test5();
  cout << "Test Passed" << endl;
  return 0;
}


text:
#include <bits/stdc++.h>
#include "../../../../../Content/C++/datastructures/trees/fenwicktrees/FenwickTree.h"
using namespace std;

void test1() {
  const auto start_time = chrono::system_clock::now();
  mt19937_64 rng(0);
  const int TESTCASES = 1e5;
  long long checkSum = 0;
  for (int ti = 0; ti < TESTCASES; ti++) {
    int N = rng() % 11;
    int M = rng() % 21;
    vector<vector<long long>> A(N, vector<long long>(M));
    for (auto &&ai : A) for (auto &&aij : ai) aij = rng() % int(1e9) + 1;
    FenwickTree<2, long long> FT(N, M);
    for (int i = 0; i < N; i++) for (int j = 0; j < M; j++) FT.update(i, j, A[i][j]);
    int Q = N == 0 || M == 0 ? 0 : 100 - rng() % 5;
    vector<long long> ans0, ans1;
    for (int i = 0; i < Q; i++) {
      int t = rng() % 2;
      if (t == 0) {
        int i = rng() % N, j = rng() % M;
        long long v = rng() % int(1e9) + 1;
        A[i][j] += v;
        FT.update(i, j, v);
      } else {
        int u = rng() % N, d = rng() % N, l = rng() % M, r = rng() % M;
        if (u > d) swap(u, d);
        if (l > r) swap(l, r);
        long long sm = 0;
        for (int j = u; j <= d; j++) for (int k = l; k <= r; k++) sm += A[j][k];
        ans0.push_back(sm);
        ans1.push_back(FT.query(u, d, l, r));
      }
    }
    assert(ans0 == ans1);
    for (auto &&a : ans0) checkSum = 31 * checkSum + a;
  }
  const auto end_time = chrono::system_clock::now();
  double sec = ((end_time - start_time).count() / double(chrono::system_clock::period::den));
  cout << "Subtest 1 Passed" << endl;
  cout << "  Time: " << fixed << setprecision(3) << sec << "s" << endl;
  cout << "  Checksum: " << checkSum << endl;
}

int main() {
  test1();
  cout << "Test Passed" << endl;
  return 0;
}

repo_name: wesley-a-leung/Resources

mime_type: text/x-c++


text:
#include <bits/stdc++.h>
#include "../../../../../Content/C++/datastructures/trees/fenwicktrees/FenwickTreeRangePoint.h"
using namespace std;

void test1() {
  const auto start_time = chrono::system_clock::now();
  mt19937_64 rng(0);
  int N = 4000;
  int M = 6000;
  vector<vector<long long>> A(N, vector<long long>(M));
  for (auto &&ai : A) for (auto &&aij : ai) aij = rng() % int(1e9) + 1;
  FenwickTreeRangePoint<2, long long> FT(N, M);
  for (int i = 0; i < N; i++) for (int j = 0; j < M; j++) FT.update(A[i][j], i, i, j, j);
  int Q = 1e6;
  vector<long long> ans;
  for (int i = 0; i < Q; i++) {
    int t = rng() % 2;
    if (t == 0) {
      int u = rng() % N, d = rng() % N, l = rng() % M, r = rng() % M;
      if (u > d) swap(u, d);
      if (l > r) swap(l, r);
      long long v = rng() % int(1e9) + 1;
      FT.update(v, u, d, l, r);
    } else {
      int i = rng() % N, j = rng() % M;
      ans.push_back(FT.get(i, j));
    }
  }
  const auto end_time = chrono::system_clock::now();
  double sec = ((end_time - start_time).count() / double(chrono::system_clock::period::den));
  cout << "Subtest 1 Passed" << endl;
  cout << "  N: " << N << endl;
  cout << "  M: " << M << endl;
  cout << "  Q: " << Q << endl;
  cout << "  Time: " << fixed << setprecision(3) << sec << "s" << endl;
  long long checkSum = 0;
  for (auto &&a : ans) checkSum = 31 * checkSum + a;
  cout << "  Checksum: " << checkSum << endl;
}

int main() {
  test1();
  cout << "Test Passed" << endl;
  return 0;
}


text:
#include <bits/stdc++.h>
#include "../../../../../Content/C++/datastructures/trees/fenwicktrees/FenwickTreeRange1D.h"
using namespace std;

void test1() {
  const auto start_time = chrono::system_clock::now();
  mt19937_64 rng(0);
  const int TESTCASES = 1e5;
  long long checkSum = 0;
  for (int ti = 0; ti < TESTCASES; ti++) {
    int N = rng() % 101;
    vector<long long> A(N);
    for (auto &&ai : A) ai = rng() % int(1e9) + 1;
    FenwickTreeRange1D<long long> FT(A);
    int Q = N == 0 ? 0 : 100 - rng() % 5;
    vector<long long> ans0, ans1;
    for (int i = 0; i < Q; i++) {
      int t = rng() % 2;
      if (t == 0) {
        int l = rng() % N, r = rng() % N;
        if (l > r) swap(l, r);
        long long v = rng() % int(1e9) + 1;
        for (int j = l; j <= r; j++) A[j] += v;
        FT.update(l, r, v);
      } else {
        int l = rng() % N, r = rng() % N;
        if (l > r) swap(l, r);
        long long sm = 0;
        for (int j = l; j <= r; j++) sm += A[j];
        ans0.push_back(sm);
        ans1.push_back(FT.query(l, r));
      }
    }
    assert(ans0 == ans1);
    for (auto &&a : ans0) checkSum = 31 * checkSum + a;
  }
  const auto end_time = chrono::system_clock::now();
  double sec = ((end_time - start_time).count() / double(chrono::system_clock::period::den));
  cout << "Subtest 1 Passed" << endl;
  cout << "  Time: " << fixed << setprecision(3) << sec << "s" << endl;
  cout << "  Checksum: " << checkSum << endl;
}

int main() {
  test1();
  cout << "Test Passed" << endl;
  return 0;
}


text:
#!/usr/bin/env bash
# first command line argument is the compile command, following arguments are
# the files to compile
numArgs=$#
declare -i total=0
declare -i pass=0
declare -i fail=0
for i in $(seq 2 $numArgs); do
  test=${!i}
  total+=1
  echo ""
  echo "$test:"
  $1 $test
  retCode=$?
  if (($retCode == 0)); then
    echo "Compiled Successfully"
    pass+=1
  else
    fail+=1
  fi
done

echo ""
echo "$total files(s) compiled"
echo "$pass successful"
echo "$fail failed"

if (($fail != 0)); then
  exit 1
fi


text:
import argparse
import pathlib

parser = argparse.ArgumentParser(
  description="Script to remove pramga once at the top of C++ header files "
              "and replace them with header guards",
)
parser.add_argument("filenames", metavar="file", type=str, nargs="+")
filenames = parser.parse_args().filenames

replaced = 0
for filename in filenames:
  print()
  print(filename + ":")
  headerguard = pathlib.Path(filename).stem.upper() + "_H"
  output = []
  with open(filename, "r") as file:
    replacePragmaOnce = False
    for line in file.read().splitlines():
      if line == "#pragma once":
        output.append("#ifndef " + headerguard)
        output.append("#define " + headerguard)
        output.append("")
        replacePragmaOnce = True
      else:
        output.append(line)
    if replacePragmaOnce:
      output.append("")
      output.append("#endif")
      print("Replaced #pragma once")
      replaced += 1
    else:
      print("Nothing to replace")
  with open(filename, "w") as file:
    file.write("\n".join(output) + "\n")

print()
print(len(filenames), "file(s) checked")
print(replaced, "file(s) modified")


text:
#!/usr/bin/env bash
# first command line argument is the compile command, following arguments are
# the files to compile
ulimit -s 524288
numArgs=$#
declare -i total=0
declare -i pass=0
declare -i fail=0
for i in $(seq 2 $numArgs); do
  test=${!i}
  total+=1
  echo ""
  echo "$test:"
  $1 $test -o $test.o
  retCode=$?
  if (($retCode == 0)); then
    timeout 5m ./$test.o
    retCode=$?
    if (($retCode == 0)); then
      pass+=1
    else
      fail+=1
      echo "Test Failed"
      echo "  Exit Code: $retCode"
    fi
    rm $test.o
  else
    fail+=1
    echo "Failed to compile"
  fi
done

echo ""
echo "$total test(s) ran"
echo "$pass passed"
echo "$fail failed"

if (($fail != 0)); then
  exit 1
fi


text:
import argparse
import sys

parser = argparse.ArgumentParser(
  description="Script to check that every line ends in \\n and not \\r\\n, "
              "does not contain \\t, "
              "and that every line does not exceed 79 characters",
)
parser.add_argument("filenames", metavar="file", type=str, nargs="+")
filenames = parser.parse_args().filenames

good = 0
bad = 0
for filename in filenames:
  print()
  print(filename + ":")
  with open(filename, "rb") as file:
    ok = True
    for curLine, line in enumerate(file, 1):
      if not line.endswith(b"\n"):
        print(f"line {curLine} does not end in \\n")
        ok = False
      if line.endswith(b"\r\n"):
        print(f"line {curLine} ends in \\r\\n")
        ok = False
      if line.find(b"\t") != -1:
        print(f"line {curLine} contains \\t")
        ok = False
      if line.find(b"http") == -1 and len(line.rstrip(b"\r\n")) > 79:
        print(f"line {curLine} exceeds maximum line length of 79")
        ok = False
    if ok:
      print("All lines good")
      good += 1
    else:
      bad += 1

print()
print(len(filenames), "file(s) checked")
print(good, "good")
print(bad, "with errors")
sys.exit(bad != 0)


text:
import java.io.*
import java.math.*
import java.util.*

class Reader {
  private val In: BufferedReader
  private var st: StringTokenizer? = null
  constructor(inputStream: InputStream) {
    In = BufferedReader(InputStreamReader(inputStream))
  }
  constructor(fileName: String) {
    In = BufferedReader(FileReader(fileName))
  }
  fun next(): String {
    while (st == null || !st!!.hasMoreTokens())
      st = StringTokenizer(In.readLine().trim())
    return st!!.nextToken()
  }
  fun nextLine(): String {
    st = null
    return In.readLine()
  }
  fun nextChar(): Char = next()[0]
  fun nextDouble(): Double = next().toDouble()
  fun nextInt(): Int = next().toInt()
  fun nextLong(): Long = next().toLong()
  fun close(): Unit = In.close()
}


text:
import java.io.*
import java.math.*
import java.util.*

class FastReader {
  private val BUFFER_SIZE: Int = 1 shl 12
  private var LENGTH: Int = -1
  private val din: DataInputStream
  private val buffer: ByteArray = ByteArray(BUFFER_SIZE)
  private var bufferPointer: Int = 0
  private var bytesRead: Int = 0
  private var buf: CharArray = CharArray(0)
  constructor(inputStream: InputStream) {
    din = DataInputStream(inputStream)
  }
  constructor(fileName: String) {
    din = DataInputStream(FileInputStream(fileName))
  }
  fun nextInt(): Int {
    var ret: Int = 0
    var c: Byte
    do {
      c = read()
    } while (c <= 32)
    var neg: Boolean = c == 45.toByte()
    if (neg) c = read()
    do {
      ret = ret * 10 + c - 48
      c = read()
    } while (c >= 48)
    if (neg) return -ret
    return ret
  }
  fun nextLong(): Long {
    var ret: Long = 0L
    var c: Byte
    do {
      c = read()
    } while (c <= 32)
    var neg: Boolean = c == 45.toByte()
    if (neg) c = read()
    do {
      ret = ret * 10 + c - 48
      c = read()
    } while (c >= 48)
    if (neg) return -ret
    return ret
  }
  fun nextDouble(): Double {
    var ret: Double = 0.0
    var div: Double = 1.0
    var c: Byte
    do {
      c = read()
    } while (c <= 32)
    var neg: Boolean = c == 45.toByte()
    if (neg) c = read()
    do {
      ret = ret * 10 + c - 48
      c = read()
    } while (c >= 48)
    if (c == 46.toByte()) {
      c = read()
      while (c >= 48) {
        div *= 10
        ret += (c - 48) / div
        c = read()
      }
    }
    if (neg) return -ret
    return ret
  }
  fun nextChar(): Char {
    var c: Byte
    do {
      c = read()
    } while (c <= 32)
    return c.toChar()
  }
  fun next(): String {
    var c: Byte
    var cnt: Int = 0
    do {
      c = read()
    } while (c <= 32)
    do {
      buf[cnt++] = c.toChar()
      c = read()
    } while (c > 32)
    return String(buf, 0, cnt)
  }
  fun nextLine(): String {
    while (bufferPointer > 0 && buffer[bufferPointer - 1] == 13.toByte())
      read()
    var c: Byte
    var cnt: Int = 0
    c = read()
    while (c != 10.toByte() && c != 0.toByte()) {
      if (c != 13.toByte()) buf[cnt++] = c.toChar()
      c = read()
    }
    return String(buf, 0, cnt)
  }
  fun setLength(len: Int) {
    LENGTH = len
    buf = CharArray(len)
  }
  fun hasNext(): Boolean {
    while (peek() > -1 && peek() <= 32) read()
    return peek() > -1
  }
  fun hasNextLine(): Boolean {
    while (peek() == 13.toByte()) read()
    return peek() > -1
  }
  private fun fillBuffer() {
    bufferPointer = 0
    bytesRead = din.read(buffer, bufferPointer, BUFFER_SIZE)
    if (bytesRead == -1) buffer[0] = -1
  }
  private fun read(): Byte {
    if (bufferPointer == bytesRead) fillBuffer()
    return buffer[bufferPointer++]
  }
  private fun peek(): Byte {
    if (bufferPointer == bytesRead) fillBuffer()
    return buffer[bufferPointer]
  }
  fun close() {
    din.close()
  }
}


text:
import java.io.*;
import java.math.*;
import java.util.*;
import java.awt.geom.*;

// Functions with polygons using the Area object
// Functions:
//   makeArea(P): creates an Area from an array of points P
//   union(areas): returns an Area that is the union of the Area
//     objects in areas
//   intersection(areas): returns an Area that is the intersection of the Area
//     objects in areas
//   getArea(area): returns the area of the Area object area
// Tested:
//   https://open.kattis.com/problems/abstractart
//   https://codeforces.com/problemset/gymProblem/100952/J
public class Polygon {
  public static Area makeArea(double[][] P) {
    Path2D.Double path = new Path2D.Double();
    path.moveTo(P[0][0], P[0][1]);
    for (int i = 1; i < P.length; i++) path.lineTo(P[i][0], P[i][1]);
    path.closePath();
    return new Area(path);
  }
  public static Area union(Area[] areas) {
    Area ret = new Area(); for (Area area : areas) ret.add(area);
    return ret;
  }
  public static Area intersection(Area[] areas) {
    if (areas.length == 0) return new Area();
    Area ret = areas[0];
    for (int i = 1; i < areas.length; i++) ret.intersect(areas[i]);
    return ret;
  }
  public static double getArea(Area area) {
    PathIterator iter = area.getPathIterator(null); double ret = 0;
    double[] buf = new double[6]; ArrayList<double[]> P = new ArrayList<>();
    for (; !iter.isDone(); iter.next()) {
      switch (iter.currentSegment(buf)) {
        case PathIterator.SEG_MOVETO:
        case PathIterator.SEG_LINETO:
          P.add(new double[]{buf[0], buf[1]}); break;
        case PathIterator.SEG_CLOSE:
          for (int i = 0; i < P.size(); i++) {
            double[] a = P.get(i), b = P.get((i + 1) % P.size());
            ret -= a[0] * b[1] - a[1] * b[0];
          }
          P = new ArrayList<>(); break;
      }
    }
    return ret / 2;
  }
}


text:
import java.io.*;
import java.math.*;
import java.util.*;

public class FastReader {
  private final int BUFFER_SIZE = 1 << 16; private int LENGTH = -1;
  private DataInputStream din; private byte[] buffer, buf;
  private int bufferPointer, bytesRead;
  public FastReader(InputStream inputStream) {
    din = new DataInputStream(inputStream); buffer = new byte[BUFFER_SIZE];
    bufferPointer = bytesRead = 0;
  }
  public FastReader(String fileName) throws IOException {
    din = new DataInputStream(new FileInputStream(fileName));
    buffer = new byte[BUFFER_SIZE]; bufferPointer = bytesRead = 0;
  }
  public int nextInt() throws IOException {
    int ret = 0; byte c; do {
      c = read();
    } while (c <= ' ');
    boolean neg = (c == '-'); if (neg) c = read();
    do {
      ret = ret * 10 + c - '0';
    } while ((c = read()) >= '0');
    if (neg) return -ret;
    return ret;
  }
  public long nextLong() throws IOException {
    long ret = 0; byte c; do {
        c = read();
    } while (c <= ' ');
    boolean neg = (c == '-');
    if (neg) c = read();
    do {
        ret = ret * 10 + c - '0';
    } while ((c = read()) >= '0');
    if (neg) return -ret;
    return ret;
  }
  public double nextDouble() throws IOException {
    double ret = 0, div = 1; byte c; do {
      c = read();
    } while (c <= ' ');
    boolean neg = (c == '-');
    if (neg) c = read();
    do {
        ret = ret * 10 + c - '0';
    } while ((c = read()) >= '0');
    if (c == '.') while ((c = read()) >= '0') ret += (c - '0') / (div *= 10);
    if (neg) return -ret;
    return ret;
  }
  public char nextChar() throws IOException {
    byte c; do {
      c = read();
    } while (c <= ' ');
    return (char) c;
  }
  public String next() throws IOException {
    int cnt = 0; byte c; do {
      c = read();
    } while (c <= ' ');
    do {
      buf[cnt++] = c;
    } while ((c = read()) > ' ');
    return new String(buf, 0, cnt);
  }
  public String nextLine() throws IOException {
    while (bufferPointer > 0 && buffer[bufferPointer - 1] == '\r') read();
    int cnt = 0; byte c;
    while ((c = read()) != '\n' && c != '\0') if (c != '\r') buf[cnt++] = c;
    return new String(buf, 0, cnt);
  }
  public void setLength(int len) { LENGTH = len; buf = new byte[len]; }
  public boolean hasNext() throws IOException {
    while (peek() != -1 && peek() <= ' ') read();
    return peek() != -1;
  }
  public boolean hasNextLine() throws IOException {
    while (peek() == '\r') read();
    return peek() != -1;
  }
  private void fillBuffer() throws IOException {
    bytesRead = din.read(buffer, bufferPointer = 0, BUFFER_SIZE);
    if (bytesRead == -1) buffer[0] = -1;
  }
  private byte read() throws IOException {
    if (bufferPointer == bytesRead) fillBuffer();
    return buffer[bufferPointer++];
  }
  private byte peek() throws IOException {
    if (bufferPointer == bytesRead) fillBuffer();
    return buffer[bufferPointer];
  }
  public void close() throws IOException {
    if (din == null) return;
    din.close();
  }
}


text:
public class Pair<Item1 extends Comparable<Item1>,
                  Item2 extends Comparable<Item2>>
    implements Comparable<Pair<Item1, Item2>> {
  public Item1 first;
  public Item2 second;
  public Pair(Item1 first, Item2 second) {
    this.first = first;
    this.second = second;
  }
  @Override
  public int hashCode() {
    return 31 * first.hashCode() + second.hashCode();
  }
  @Override
  public boolean equals(Object o) {
    if (o == this) return true;
    if (!(o instanceof Pair)) return false;
    Pair p = (Pair) o;
    return p.first.equals(first) && p.second.equals(second);
  }
  @Override
  public int compareTo(Pair<Item1, Item2> p) {
    int o = first.compareTo(p.first);
    return o == 0 ? second.compareTo(p.second) : o;
  }
}


text:
import java.io.*;
import java.math.*;
import java.util.*;

public class Reader {
  private BufferedReader in; private StringTokenizer st;
  public Reader(InputStream inputStream) {
    in = new BufferedReader(new InputStreamReader(inputStream));
  }
  public Reader(String fileName) throws FileNotFoundException {
    in = new BufferedReader(new FileReader(fileName));
  }
  public String next() throws IOException {
    while (st == null || !st.hasMoreTokens())
      st = new StringTokenizer(in.readLine().trim());
    return st.nextToken();
  }
  public String nextLine() throws IOException {
    st = null; return in.readLine();
  }
  public BigInteger nextBigInteger() throws IOException {
    return new BigInteger(next());
  }
  public byte nextByte() throws IOException {
    return Byte.parseByte(next());
  }
  public char nextChar() throws IOException { return next().charAt(0); }
  public double nextDouble() throws IOException {
    return Double.parseDouble(next());
  }
  public float nextFloat() throws IOException {
    return Float.parseFloat(next());
  }
  public int nextInt() throws IOException {
    return Integer.parseInt(next());
  }
  public long nextLong() throws IOException {
    return Long.parseLong(next());
  }
  public short nextShort() throws IOException {
    return Short.parseShort(next());
  }
  public void close() throws IOException { in.close(); }
}


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

// Cuthill-McKee reordering of a graph to reduce the bandwith
// 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
//   rev: a boolean indicating whether the order should be reversed or not,
//     performs better with Gaussian Elimination
// Fields:
//   mapping: a mapping from the original vertex to the reordered vertex
//   revMapping: a mapping from the reordered vertex to the original vertex
// In practice, has a small constant
// Time Complexity:
//   constructor: O(V log V + E log E)
// Memory Complexity: O(V)
// Tested:
//   https://dmoj.ca/problem/ddrp3
struct CuthillMcKee {
  int V; vector<int> mapping, revMapping, deg;
  template <class Graph> CuthillMcKee(const Graph &G, bool rev = false)
      : V(G.size()), mapping(V, 0), revMapping(V), deg(V, 0) {
    for (int v = 0; v < V; v++) for (int w : G[v]) { deg[v]++; deg[w]++; }
    auto cmpDeg = [&] (int v, int w) {
      return make_pair(deg[v], v) < make_pair(deg[w], w);
    };
    int front = 0, back = 0; vector<int> P(V); iota(P.begin(), P.end(), 0);
    sort(P.begin(), P.end(), cmpDeg); for (int i = 0; i < V; i++) {
      int s = P[i]; if (mapping[s]) continue;
      mapping[revMapping[back++] = s] = 1; while (front < back) {
        int v = revMapping[front++]; vector<int> adj; adj.reserve(deg[v]);
        for (int w : G[v]) adj.push_back(w);
        sort(adj.begin(), adj.end(), cmpDeg);
        for (int w : adj) if (!mapping[w]) mapping[revMapping[back++] = w] = 1;
      }
    }
    if (rev) reverse(revMapping.begin(), revMapping.end());
    for (int i = 0; i < V; i++) mapping[revMapping[i]] = i;
  }
};


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

// Breadth First Traversal of a graph (weighted or unweighted)
// Vertices are 0-indexed
// Template Arguments:
//   T: the type of the weight of the edges in the graph
// Constructor Arguments:
//   G: a generic 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
// 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 breadth first search 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://www.spoj.com/problems/BITMAP/
//   https://dmoj.ca/problem/ddrp3
template <class T = int> struct BFS {
  using Edge = tuple<int, int, T>; vector<T> dist; vector<int> par; T INF;
  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 Graph> BFS(const Graph &G, const vector<int> &srcs,
                             T INF = numeric_limits<T>::max())
      : dist(G.size(), INF), par(G.size(), -1), INF(INF) {
    vector<int> q(G.size()); int front = 0, back = 0;
    for (int s : srcs) dist[q[back++] = s] = T();
    while (front < back) {
      int v = q[front++]; for (auto &&e : G[v]) {
        int w = getTo(e); if (dist[w] >= INF)
          dist[q[back++] = w] = dist[par[w] = v] + getWeight(e);
      }
    }
  }
  template <class Graph> BFS(const Graph &G, int s,
                             T INF = numeric_limits<T>::max())
      : BFS(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;

// Multidirectional breadth first search of an undirected, unweighted graph to
//   query for the closest distance between two distinct vertices in a set
//   of vertices
// Vertices are 0-indexed
// Template Arguments:
//   Graph: 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
// Constructor Arguments:
//   G: a instance of Graph
// Functions:
//   closest(src): returns the closest distance between any two distinct
//     vertices in src
// In practice, has a moderate constant
// Time Complexity:
//   constructor: O(V + E)
//   closest: O(V + E) worst case, much better in practice if two source
//     vertices are connected and the graph is random and sparse
// Memory Complexity: O(V)
// Tested:
//   https://dmoj.ca/problem/acc1p2
//   https://dmoj.ca/problem/wac1p6
template <class Graph> struct MultidirectionalBFS {
  Graph G; vector<int> d, from, vis, q; int stamp;
  MultidirectionalBFS(const Graph &G)
      : G(G), d(G.size(), INT_MAX), from(G.size(), -1), vis(G.size(), 0),
        q(G.size()), stamp(0) {}
  int closest(const vector<int> &srcs) {
    assert(int(srcs.size()) >= 2);
    int front = 0, back = 0; ++stamp; for (int s : srcs) {
      if (vis[s] == stamp) return 0;
      vis[q[back++] = from[s] = s] = stamp; d[s] = 0;
    }
    int mn = INT_MAX; while (front < back) {
      int v = q[front++]; for (int w : G[v]) {
        if (vis[v] != vis[w]) {
          d[q[back++] = w] = d[v] + 1; vis[w] = vis[v]; from[w] = from[v];
        } else if (from[v] != from[w]) {
          int sm = d[v] + d[w] + 1; if (sm > mn) return mn;
          mn = sm;
        }
      }
    }
    return mn;
  }
};


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

// Computes the Transitive Closure in a graph using Floyd Warshall
//   with bitset optimizations
// Template Arguments:
//   MAXV: the maximum number of vertices in the graph
// Constructor Arguments:
//   matrix: a V x V matrix (represented by V bitsets of size MAXV) such that
//     matrix[v][w] is 1 if there is a directed edge from v to w and
//       0 otherwise
// Functions:
//   reachable(v, w): returns true if w is reachable from v and false otherwise
// In practice, has a very small constant, slower than the SCC variant
// Time Complexity:
//   constructor: O(V^3 / 64)
//   reachable: O(1)
// Memory Complexity: O(MAXV V / 64)
// Tested:
//   Stress Tested
//   https://dmoj.ca/problem/acc2p2
//   https://open.kattis.com/problems/watchyourstep
template <const int MAXV> struct TransitiveClosureFloydWarshall {
  vector<bitset<MAXV>> dp;
  TransitiveClosureFloydWarshall(vector<bitset<MAXV>> matrix)
      : dp(move(matrix)) {
    int V = dp.size(); for (int v = 0; v < V; v++) dp[v][v] = 1;
    for (int u = 0; u < V; u++) for (int v = 0; v < V; v++)
      if (dp[v][u]) dp[v] |= dp[u];
  }
  bool reachable(int v, int w) { return dp[v][w]; }
};


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

// Maintains the Transitive Closure in a graph after edges are added
// Constructor Arguments:
//   V: the number of vertices in the graph
// Functions:
//   reachable(v, w): returns true if w is reachable from v and false otherwise
//   addEdge(v, w): adds an edge between vertices v and w, returns true if the
//     resulting edge does not create a new cycle and false otherwise
// In practice, has a small constant
// Time Complexity:
//   constructor: O(V^2)
//   reachable: O(1)
//   addEdge: O(V) amortized
// Memory Complexity: O(V^2)
// Tested:
//   https://www.spoj.com/problems/GHOSTS/
struct IncrementalTransitiveClosure {
  int V; vector<vector<int>> par; vector<vector<vector<int>>> ch;
  IncrementalTransitiveClosure(int V)
      : V(V), par(V, vector<int>(V, -1)), ch(V, vector<vector<int>>(V)) {}
  bool reachable(int v, int w) { return v == w || par[v][w] >= 0; }
  void meld(int root, int sub, int v, int w) {
    par[root][w] = v; ch[root][v].push_back(w);
    for (int c : ch[sub][w]) if (!reachable(root, c)) meld(root, sub, w, c);
  }
  bool addEdge(int v, int w) {
    if (reachable(w, v)) return false;
    if (reachable(v, w)) return true;
    for (int u = 0; u < V; u++) if (reachable(u, v) && !reachable(u, w))
      meld(u, w, v, w);
    return true;
  }
};


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

// Computes the Topological Order of a directed acyclic graph
// Vertices are 0-indexed
// Constructor Arguments:
//   G: a generic directed acyclic graph structure
//       which can be weighted or unweighted, though weights do not change
//       the topological order
//     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
// Fields:
//   ind: a vector of the topological order index for each vertex
//   ord: a vector of the vertices sorted by topological order index
// In practice, has a moderate constant, faster than DepthFirstOrder
// Time Complexity:
//   constructor: O(V + E)
// Memory Complexity: O(V)
// Tested:
//   https://atcoder.jp/contests/dp/tasks/dp_g
struct TopologicalOrder {
  int V; vector<int> ind, ord, inDeg;
  int getTo(int e) { return e; }
  template <class T> int getTo(const pair<int, T> &e) { return e.first; }
  template <class DAG> TopologicalOrder(const DAG &G)
      : V(G.size()), ind(V, -1), ord(V), inDeg(V, 0) {
    int front = 0, back = 0;
    for (int v = 0; v < V; v++) for (auto &&e : G[v]) inDeg[getTo(e)]++;
    for (int v = 0; v < V; v++) if (inDeg[v] == 0) ord[back++] = v;
    while (front < back) {
      int v = ord[front]; ind[v] = front++;
      for (auto &&e : G[v]) if (--inDeg[getTo(e)] == 0) ord[back++] = getTo(e);
    }
  }
};


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

// Computes the Transitive Closure in a graph using strongly connected
//   components and dynamic programming with bitset optimizations
// Template Arguments:
//   MAXV: the maximum number of vertices in the graph
// Constructor Arguments:
//   G: a generic 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
// Functions:
//   reachable(v, w): returns true if w is reachable from v and false otherwise
// In practice, has a very small constant, faster than the
//   Floyd Warshall variant
// Time Complexity:
//   constructor: O(V + E + MAXV E / 64)
//   reachable: O(1)
// Memory Complexity: O(V + E + MAXV V / 64)
// Tested:
//   Stress Tested
//   https://dmoj.ca/problem/acc2p2
//   https://open.kattis.com/problems/watchyourstep
template <const int MAXV> struct TransitiveClosureSCC {
  vector<pair<int, int>> DAG; SCC scc; vector<bitset<MAXV>> dp;
  template <class Graph> TransitiveClosureSCC(const Graph &G)
      : DAG(), scc(G, DAG), dp(scc.components.size()) {
    for (int i = 0; i < int(dp.size()); i++) dp[i][i] = 1;
    for (auto &&e : DAG) dp[e.first] |= dp[e.second];
  }
  bool reachable(int v, int w) { return dp[scc.id[v]][scc.id[w]]; }
};


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

// Finds the centroids for each connected component in a forest
// Vertices are 0-indexed
// Function Arguments:
//   G: a generic forest 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 forest
// Return Value: a vector of vectors of integers containing the centroids for
//   each connected component with each connected component having either one
//   or two centroids; connected components are ordered by their minimum vertex
// In practice, has a moderate constant
// Time Complexity: O(V)
// Memory Complexity: O(V)
// Tested:
//   https://codeforces.com/contest/1406/problem/C
template <class Forest> vector<vector<int>> getCentroids(const Forest &G) {
  int V = G.size(); vector<int> size(V, 0); vector<vector<int>> centroids;
  function<int(int, int)> getSize = [&] (int v, int prev) {
    size[v] = 1; for (int w : G[v]) if (w != prev) size[v] += getSize(w, v);
    return size[v];
  };
  function<void(int, int, int)> dfs = [&] (int v, int prev, int compSize) {
    bool b = true; for (int w : G[v])
      if (w != prev) { dfs(w, v, compSize); b &= size[w] <= compSize / 2; }
    if (b && compSize - size[v] <= compSize / 2) centroids.back().push_back(v);
  };
  for (int v = 0; v < V; v++)
    if (size[v] == 0) { centroids.emplace_back(); dfs(v, -1, getSize(v, -1)); }
  return centroids;
}


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

// Computes the diameter of a undirected tree (weighted or unweighted)
// Vertices are 0-indexed
// Template Arguments:
//   T: the type of the weight of the edges in the tree
// Constructor Arguments:
//   G: a generic undirected tree 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 tree, or a list of
//         pair<int, T> for a weighted tree with weights of type T)
//       size() const: returns the number of vertices in the tree
//   INF: a value for infinity
// Fields:
//   endpoints: a pair containing the vertices on the endpoints of the diameter
//   diameter: the length of the diameter
// Functions:
//   getPath(): returns the list of edges on the diameter
// In practice, has a moderate constant
// Time Complexity:
//   constructor: O(V)
//   getPath: O(V)
// Memory Complexity: O(V)
// Tested:
//   https://judge.yosupo.jp/problem/tree_diameter
template <class T = int> struct TreeDiameter {
  BFS<T> bfs; pair<int, int> endpoints; T diameter;
  template <class Tree>
  TreeDiameter(const Tree &G, T INF = numeric_limits<T>::max())
      : bfs(G, 0, INF) {
    endpoints.first = max_element(bfs.dist.begin(), bfs.dist.end())
        - bfs.dist.begin();
    move(bfs); bfs = BFS<T>(G, endpoints.first);
    endpoints.second = max_element(bfs.dist.begin(), bfs.dist.end())
        - bfs.dist.begin();
    diameter = bfs.dist[endpoints.second];
  }
  vector<typename BFS<T>::Edge> getPath() {
    return bfs.getPath(endpoints.second);
  }
};


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

// Bidirectional breadth first search of an undirected, unweighted graph to
//   query for the distance between two vertices
// Vertices are 0-indexed
// Template Arguments:
//   Graph: 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
// Constructor Arguments:
//   G: a instance of Graph
// Functions:
//   dist(s, t): returns the distance between vertices s and t
// In practice, has a moderate constant
// Time Complexity:
//   constructor: O(V + E)
//   dist: O(V + E) worst case, much better in practice if s and t are
//     connected and the graph is random and sparse
// Memory Complexity: O(V)
// Tested:
//   https://dmoj.ca/problem/acc1p2
template <class Graph> struct BidirectionalBFS {
  Graph G; vector<int> d, vis, q; int stamp;
  BidirectionalBFS(const Graph &G)
      : G(G), d(G.size(), INT_MAX), vis(G.size(), 0), q(G.size()), stamp(0) {}
  int dist(int s, int t) {
    if (s == t) return 0;
    int front = 0, back = 0; d[s] = d[t] = 0;
    vis[q[back++] = s] = ++stamp; vis[q[back++] = t] = -stamp;
    while (front < back) {
      int v = q[front++]; for (int w : G[v]) {
        if (vis[v] == -vis[w]) return d[v] + d[w] + 1;
        else if (vis[v] != vis[w]) {
          d[q[back++] = w] = d[v] + 1; vis[w] = vis[v];
        }
      }
    }
    return INT_MAX;
  }
};


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

// Breadth First Traversal of a graph where all edge weights are 0 or 1
// 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
// Fields:
//   dist: vector of the shortest distance from the closest source vertex to
//     each vertex, or INT_MAX if unreachable
//   par: the parent vertex for each vertex in the breadth first search 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://www.codechef.com/problems/REVERSE
//   https://atcoder.jp/contests/abc176/tasks/abc176_d
struct ZeroOneBFS {
  using Edge = tuple<int, int, int>; vector<int> dist, par;
  template <class WeightedGraph>
  ZeroOneBFS(const WeightedGraph &G, const vector<int> &srcs)
      : dist(G.size(), INT_MAX), par(G.size(), -1) {
    vector<int> q(G.size()), stk(G.size()); int front = 0, back = 0, top = 0;
    for (int s : srcs) dist[stk[top++] = s] = 0;
    while (top > 0 || front < back) {
      int v = top > 0 ? stk[--top] : q[front++]; for (auto &&e : G[v]) {
        int w = e.first; if (dist[w] > dist[v] + e.second) {
          (e.second == 0 ? stk[top++] : q[back++]) = w;
          dist[w] = dist[par[w] = v] + e.second;
        }
      }
    }
  }
  template <class WeightedGraph> ZeroOneBFS(const WeightedGraph &G, int s)
      : ZeroOneBFS(G, vector<int>{s}) {}
  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 Depth First Orders of a graph (pre order, post order,
//   topological/reverse post order)
// Vertices are 0-indexed
// Constructor Arguments:
//   G: a generic 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
//   rt: a single root vertex
//   roots: a vector of root vertices for each connected component
// Fields:
//   preInd: a vector of the pre order index for each vertex
//   postInd: a vector of the post order index for each vertex
//   revPostInd: a vector of the topological/reverse post order index
//     for each vertex
//   preVert: a vector of the vertices sorted by pre order index
//   postVert: a vector of the vertices sorted by post order index
//   revPostVert: a vector of the vertices sorted by topological/reverse post
//     order index
// In practice, has a moderate constant, slower than TopologicalOrder
// Time Complexity:
//   constructor: O(V + E)
// Memory Complexity: O(V)
// Tested:
//   https://atcoder.jp/contests/nikkei2019-qual/tasks/nikkei2019_qual_d
//   https://codeforces.com/contest/24/problem/A
struct DFSOrder {
  int V, curPre, curPost, curRevPost;
  vector<int> preInd, postInd, revPostInd, preVert, postVert, revPostVert;
  template <class Graph> void dfs(const Graph &G, int v) {
    preVert[preInd[v] = curPre++] = v;
    for (int w : G[v]) if (preInd[w] == -1) dfs(G, w);
    postVert[postInd[v] = curPost++] = v;
    revPostVert[revPostInd[v] = curRevPost--] = v;
  }
  template <class Graph>
  DFSOrder(const Graph &G, const vector<int> &roots = vector<int>())
      : V(G.size()), curPre(0), curPost(0), curRevPost(V - 1), preInd(V, -1),
        postInd(V), revPostInd(V), preVert(V), postVert(V), revPostVert(V) {
    if (roots.empty()) {
      for (int v = 0; v < V; v++) if (preInd[v] == -1) dfs(G, v);
    } else for (int v : roots) if (preInd[v] == -1) dfs(G, v);
  }
  template <class Graph> DFSOrder(const Graph &G, int rt)
      : DFSOrder(G, vector<int>{rt}) {}
};


text:
#pragma once
#include <bits/stdc++.h>
#include <ext/pb_ds/priority_queue.hpp>
using namespace std;
using namespace __gnu_pbds;

// Computes the maximum flow using a path with the minimum cost by finding
//   Shortest Augmenting Paths
// Can handle negative edges but not negative cost cycles
//   (assertion failure if there is a negative cost cycle)
// Vertices are 0-indexed
// Template Arguments:
//   FlowUnit: the type of the flow
//   CostUnit: the type of the cost, must be integral
// Constructor Arguments:
//   V: the number of vertices in the flow network
//   FLOW_EPS: a value for the flow epsilon
//   COST_INF: a value for the cost infinity
// Fields:
//   V: the number of vertices in the flow network
//   FLOW_EPS: a value for the flow epsilon
//   COST_INF: a value for the cost infinity
//   G: an adjacency list of all edges and reverse edges in the flow network
// Functions:
//   addEdge(v, w, vwCap, vwCost): adds an edge from v to w with capacity
//     vwCap and cost vwCost
//   getFlowMinCost(s, t): returns the maximum flow from s to t with
//     the minimum cost
// In practice, has a very small constant
// Time Complexity:
//   constructor: O(V)
//   addEdge: O(1)
//   getFlowMinCost: O(E^2 V log V), much faster in practice
// Memory Complexity: O(V + E)
// Tested:
//   https://open.kattis.com/problems/mincostmaxflow
//   https://loj.ac/p/102
//   https://www.spoj.com/problems/SCITIES/
//   https://dmoj.ca/problem/tle17c7p5
template <class FlowUnit, class CostUnit> struct SAPMinCostMaxFlow {
  struct Edge {
    int to, rev; FlowUnit cap, resCap; CostUnit cost;
    Edge(int to, int rev, FlowUnit cap, CostUnit cost)
        : to(to), rev(rev), cap(cap), resCap(cap), cost(cost) {}
  };
  struct Node {
    CostUnit d; int v; Node(CostUnit d, int v) : d(d), v(v) {}
    bool operator < (const Node &o) const { return d > o.d; }
  };
  int V; FlowUnit FLOW_EPS; CostUnit COST_INF; vector<vector<Edge>> G;
  bool hasNegativeCost; vector<CostUnit> phi, dist;
  vector<Edge*> to; vector<int> par;
  SAPMinCostMaxFlow(int V, FlowUnit FLOW_EPS = FlowUnit(1e-9),
                    CostUnit COST_INF = numeric_limits<CostUnit>::max())
      : V(V), FLOW_EPS(FLOW_EPS), COST_INF(COST_INF), G(V),
        hasNegativeCost(false), phi(V), dist(V), to(V), par(V) {}
  void addEdge(int v, int w, FlowUnit vwCap, CostUnit vwCost) {
    if (v == w) return;
    if (vwCost < CostUnit()) hasNegativeCost = true;
    G[v].emplace_back(w, int(G[w].size()), vwCap, vwCost);
    G[w].emplace_back(v, int(G[v].size()) - 1, FlowUnit(), -vwCost);
  }
  void bellmanFord(int s) {
    fill(phi.begin(), phi.end(), COST_INF); phi[s] = CostUnit();
    for (int j = 0; j < V - 1; j++) for (int v = 0; v < V; v++)
      for (auto &&e : G[v]) if (e.resCap > FLOW_EPS && phi[v] < COST_INF)
        phi[e.to] = min(phi[e.to], phi[v] + e.cost);
    for (int v = 0; v < V; v++) for (auto &&e : G[v])
      if (e.resCap > FLOW_EPS && phi[v] < COST_INF)
        assert(phi[e.to] <= phi[v] + e.cost);
  }
  bool dijkstra(int s, int t) {
    fill(dist.begin(), dist.end(), COST_INF); fill(par.begin(), par.end(), -1);
    vector<bool> seen(V, false); __gnu_pbds::priority_queue<Node> PQ;
    vector<typename decltype(PQ)::point_iterator> ptr(V, PQ.end());
    ptr[s] = PQ.push(Node(dist[s] = CostUnit(), s)); while (!PQ.empty()) {
      int v = PQ.top().v; PQ.pop(); ptr[v] = PQ.end(); seen[v] = true;
      for (auto &&e : G[v]) {
        int w = e.to; if (seen[w] || e.resCap <= FLOW_EPS) continue;
        CostUnit d = dist[v] + e.cost + phi[v] - phi[w];
        if (dist[w] <= d) continue;
        par[w] = v; to[w] = &e;
        if (ptr[w] == PQ.end()) ptr[w] = PQ.push(Node(dist[w] = d, w));
        else PQ.modify(ptr[w], Node(dist[w] = d, w));
      }
    }
    return dist[t] < COST_INF;
  }
  pair<FlowUnit, CostUnit> getFlowMinCost(int s, int t) {
    pair<FlowUnit, CostUnit> ret = make_pair(FlowUnit(), CostUnit());
    fill(phi.begin(), phi.end(), CostUnit()); if (s == t) return ret;
    if (hasNegativeCost) bellmanFord(s);
    while (dijkstra(s, t)) {
      FlowUnit aug = FlowUnit(); for (int v = t; par[v] != -1; v = par[v])
        if (v == t || aug > to[v]->resCap) aug = to[v]->resCap;
      ret.first += aug; for (int v = t; par[v] != -1; v = par[v]) {
        to[v]->resCap -= aug; G[to[v]->to][to[v]->rev].resCap += aug;
        ret.second += aug * to[v]->cost;
      }
      for (int v = 0; v < V; v++) if (dist[v] < COST_INF) phi[v] += dist[v];
    }
    return ret;
  }
};


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

// Given a list of weighted edges representing an undirected flow graph,
//   compute a tree such that the max flow/min cut between any pair of vertices
//   is given by the minimum weight edge between the two vertices
// Template Arguments:
//   T: the type of the weight of the edges
// Constructor Arguments:
//   V: number of vertices in the graph
//   edges: a vector of tuples in the form (v, w, weight) representing
//     an undirected edge in the graph between vertices v and w with
//     weight of weight
// Fields:
//   treeEdges: a vector of tuples of the edges in the Gomory-Hu Tree
// In practice, has a very small constant
// Time Complexity:
//   constructor: O(V^3 sqrt E), much faster in practice
// Memory Complexity: O(V + E)
// Tested:
//   https://codeforces.com/gym/101480/problem/J
template <class T> struct GomoryHu {
  using Edge = tuple<int, int, T>; vector<Edge> treeEdges;
  GomoryHu(int V, const vector<Edge> &edges, T EPS = T(1e-9)) {
    PushRelabelMaxFlow<FlowEdge<T>> mf(V, EPS); for (auto &&e : edges)
      mf.addEdge(get<0>(e), get<1>(e), get<2>(e), get<2>(e));
    vector<int> par(V, 0); for (int i = 1; i < V; i++) {
      for (int v = 0; v < V; v++) for (auto &&e : mf.G[v]) e.resCap = e.cap;
      treeEdges.emplace_back(i, par[i], mf.getFlow(i, par[i]));
      for (int j = i + 1; j < V; j++)
        if (par[j] == par[i] && mf.cut[j]) par[j] = i;
    }
  }
};


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

// A sample edge struct for maximum flow
// Flow of edge can be found with cap - resCap
// Template Arguments:
//   _FlowUnit: the type of the flow
// Constructor Arguments:
//   to: the vertex that this directed edge ends at
//   rev: the index in the adjacency list of vertex to of the reverse edge
//   cap: the initial capacity of this edge
// Fields:
//   FlowUnit: the type of the flow
//   to: the vertex that this directed edge ends at
//   rev: the index in the adjacency list of vertex to of the reverse edge
//   cap: the initial capacity of this edge
//   resCap: the residual capacity of this edge
// Tested:
//   https://open.kattis.com/problems/maxflow
//   https://open.kattis.com/problems/mincut
//   https://www.spoj.com/problems/FASTFLOW/
//   https://vn.spoj.com/problems/FFLOW/
//   https://loj.ac/p/127
template <class _FlowUnit> struct FlowEdge {
  using FlowUnit = _FlowUnit; int to, rev; FlowUnit cap, resCap;
  FlowEdge(int to, int rev, FlowUnit cap)
      : to(to), rev(rev), cap(cap), resCap(cap) {}
};

// Computes the maximum flow and minimum cut in a flow network using the
//   Push Relabel algorithm with the highest label selection rule and
//   gap relabelling heuristic
// Vertices are 0-indexed
// Template Arguments:
//   Edge: a generic edge class (such as FlowEdge)
//     Required Fields:
//       FlowUnit: the type of the flow
//       to: the vertex that this directed edge ends at
//       rev: the index in the adjacency list of vertex to of the reverse edge
//       resCap: the residual capacity of this edge
//     Required Functions:
//       constructor(to, rev, cap): initializes the edge to vertex to with
//         the reverse index rev and capacity cap
// Constructor Arguments:
//   V: the number of vertices in the flow network
//   FLOW_EPS: a value for the flow epsilon
// Fields:
//   V: the number of vertices in the flow network
//   FLOW_EPS: a value for the flow epsilon
//   G: an adjacency list of all edges and reverse edges in the flow network
//   cut: a vector of booleans indicating which side of a minimum s-t cut it
//     is on (true if on the same side as s, false if on the same side as t)
// Functions:
//   addEdge(v, w, vwCap, wvCap): adds an edge from v to w with capacity
//     vwCap and a reverse capacity of wvCap
//   getFlow(s, t): returns the maximum flow from s to t as well as the minimum
//     s-t cut
// In practice, has a very small constant
// Time Complexity:
//   constructor: O(V)
//   addEdge: O(1)
//   getFlow: O(V^2 sqrt E), much faster in practice
// Memory Complexity: O(V + E)
// Tested:
//   https://open.kattis.com/problems/maxflow
//   https://open.kattis.com/problems/mincut
//   https://www.spoj.com/problems/FASTFLOW/
//   https://vn.spoj.com/problems/FFLOW/
//   https://loj.ac/p/127
template <class Edge> struct PushRelabelMaxFlow {
  using FlowUnit = typename Edge::FlowUnit;
  int V; FlowUnit FLOW_EPS; vector<vector<Edge>> G; vector<bool> cut;
  PushRelabelMaxFlow(int V, FlowUnit FLOW_EPS = FlowUnit(1e-9))
      : V(V), FLOW_EPS(FLOW_EPS), G(V) {}
  void addEdge(int v, int w, FlowUnit vwCap, FlowUnit wvCap = FlowUnit()) {
    if (v == w) return;
    G[v].emplace_back(w, int(G[w].size()), vwCap);
    G[w].emplace_back(v, int(G[v].size()) - 1, wvCap);
  }
  FlowUnit getFlow(int s, int t) {
    cut.assign(V, false); if (s == t) return FlowUnit();
    vector<FlowUnit> ex(V, FlowUnit()); vector<int> h(V, 0), cnt(V * 2, 0);
    vector<vector<int>> hs(V * 2);
    vector<typename vector<Edge>::iterator> cur(V);
    auto push = [&] (int v, Edge &e, FlowUnit df) {
      int w = e.to;
      if (abs(ex[w]) <= FLOW_EPS && df > FLOW_EPS) hs[h[w]].push_back(w);
      e.resCap -= df; G[w][e.rev].resCap += df; ex[v] -= df; ex[w] += df;
    };
    h[s] = V; ex[t] = FlowUnit(1); cnt[0] = V - 1;
    for (int v = 0; v < V; v++) cur[v] = G[v].begin();
    for (auto &&e : G[s]) push(s, e, e.resCap);
    if (!hs[0].empty()) for (int hi = 0; hi >= 0;) {
      int v = hs[hi].back(); hs[hi].pop_back(); while (ex[v] > FLOW_EPS) {
        if (cur[v] == G[v].end()) {
          h[v] = INT_MAX; for (auto e = G[v].begin(); e != G[v].end(); e++)
            if (e->resCap > FLOW_EPS && h[v] > h[e->to] + 1)
              h[v] = h[(cur[v] = e)->to] + 1;
          cnt[h[v]]++;
          if (--cnt[hi] == 0 && hi < V) for (int w = 0; w < V; w++)
            if (hi < h[w] && h[w] < V) { cnt[h[w]]--; h[w] = V + 1; }
          hi = h[v];
        } else if (cur[v]->resCap > FLOW_EPS && h[v] == h[cur[v]->to] + 1) {
          push(v, *cur[v], min(ex[v], cur[v]->resCap));
        } else cur[v]++;
      }
      while (hi >= 0 && hs[hi].empty()) hi--;
    }
    for (int v = 0; v < V; v++) cut[v] = h[v] >= V;
    return -ex[s];
  }
};


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

// A sample edge struct for flow with demands
// Flow of edge can be found with cap - resCap
// Template Arguments:
//   _FlowUnit: the type of the flow
// Constructor Arguments:
//   to: the vertex that this directed edge ends at
//   rev: the index in the adjacency list of vertex to of the reverse edge
//   dem: the initial demand of this edge
//   cap: the initial capacity of this edge
// Fields:
//   FlowUnit: the type of the flow
//   to: the vertex that this directed edge ends at
//   rev: the index in the adjacency list of vertex to of the reverse edge
//   dem: the initial demand of this edge
//   cap: the initial capacity of this edge
//   resCap: the residual capacity of this edge
// Tested:
//   https://loj.ac/p/117
//   https://dmoj.ca/problem/wac4p6
//   https://codeforces.com/gym/100199/problem/B
template <class _FlowUnit> struct FlowEdgeDemands {
  using FlowUnit = _FlowUnit; int to, rev; FlowUnit dem, cap, resCap;
  FlowEdgeDemands(int to, int rev, FlowUnit dem, FlowUnit cap)
      : to(to), rev(rev), dem(dem), cap(cap), resCap(cap) {}
};

// Computes the minimum/feasible/maximum flow in a flow network with
//   edge demands using the Push Relabel algorithm with the highest label
//   selection rule and gap relabelling heuristic
// Vertices are 0-indexed
// Template Arguments:
//   Edge: a generic edge class (such as FlowEdgeDemands)
//     Required Fields:
//       FlowUnit: the type of the flow
//       to: the vertex that this directed edge ends at
//       rev: the index in the adjacency list of vertex to of the reverse edge
//       dem: the initial demand of this edge
//       cap: the initial capacity of this edge
//       resCap: the residual capacity of this edge
//     Required Functions:
//       constructor(to, rev, dem, cap): initializes the edge to vertex to with
//         the reverse index rev, demand dem, and capacity cap
// Constructor Arguments:
//   V: the number of vertices in the flow network
//   FLOW_INF: a value for the flow infinity
//   FLOW_EPS: a value for the flow epsilon
// Fields:
//   V: the number of vertices in the flow network
//   FLOW_INF: a value for the flow infinity
//   FLOW_EPS: a value for the flow epsilon
//   G: an adjacency list of all edges and reverse edges in the flow network
// Functions:
//   addEdge(v, w, vwDem, wvCap): adds an edge from v to w with capacity
//     vwCap and a demand of vwDem
//   getFeasibleFlow(s, t, flowType): computes a feasible flow (or circulation
//     if s and t are -1), maximizes the flow value if flowType is positive,
//     minimizes it if flowType is negative, and finds any flow if 0; returns
//     a pair of a boolean and FlowUnit indicating whether a flow is feasible
//     and the flow value
//   getMinFlow(s, t): returns a pair of a boolean and FlowUnit indicating
//     whether a flow is feasible and the minimum flow value
//   getMaxFlow(s, t): returns a pair of a boolean and FlowUnit indicating
//     whether a flow is feasible and the maximum flow value
// In practice, has a very small constant
// Time Complexity:
//   constructor: O(V)
//   addEdge: O(1)
//   getFeasibleFlow, getMinFlow, getMaxFlow: O(V^2 sqrt E),
//     much faster in practice
// Memory Complexity: O(V + E)
// Tested:
//   https://loj.ac/p/117
//   https://dmoj.ca/problem/wac4p6
//   https://codeforces.com/gym/100199/problem/B
template <class Edge>
struct PushRelabelFlowDemands : public PushRelabelMaxFlow<Edge> {
  using MF = PushRelabelMaxFlow<Edge>;
  using FlowUnit = typename Edge::FlowUnit; using MF::V; using MF::G;
  FlowUnit FLOW_INF; vector<FlowUnit> outDem, inDem;
  PushRelabelFlowDemands(
      int V, FlowUnit FLOW_INF = numeric_limits<FlowUnit>::max(),
      FlowUnit FLOW_EPS = FlowUnit(1e-9))
      : MF(V, FLOW_EPS), FLOW_INF(FLOW_INF),
        outDem(V, FlowUnit()), inDem(V, FlowUnit()) {}
  void addEdge(int v, int w, FlowUnit vwDem, FlowUnit vwCap, int type = 1) {
    if (v == w) return;
    G[v].emplace_back(w, int(G[w].size()), vwDem, vwCap);
    G[w].emplace_back(v, int(G[v].size()) - 1, -vwDem, -vwDem);
    if (type == 1) { outDem[v] += vwDem; inDem[w] += vwDem; }
  }
  pair<bool, FlowUnit> getFeasibleFlow(int s = -1, int t = -1,
                                       int flowType = 0) {
    int ss = V, tt = V + 1; G.emplace_back(); G.emplace_back();
    FlowUnit bnd = FLOW_INF, sm = 0;
    pair<bool, FlowUnit> ret(true, flowType < 0 ? FLOW_INF : FlowUnit());
    for (int v = 0; v < V; v++) {
      for (auto &&e : G[v]) e.cap -= e.dem;
      addEdge(ss, v, FlowUnit(), inDem[v], 2);
      addEdge(v, tt, FlowUnit(), outDem[v], 2); sm += inDem[v];
    }
    for (int h = 0; h < 2; h++) {
      for (int v = 0; v < V + 2; v++) for (auto &&e : G[v]) e.resCap = e.cap;
      if (s != -1 && t != -1) addEdge(t, s, FlowUnit(), bnd, 2);
      V += 2;
      if (sm - (bnd = MF::getFlow(ss, tt)) > MF::FLOW_EPS) ret.first = false;
      V -= 2; if (s != -1 && t != -1) {
        G[s].pop_back();
        if (ret.first) ret.second = G[t].back().cap - G[t].back().resCap;
        G[t].pop_back();
      }
      if (flowType >= 0 || !ret.first || h > 0) break;
      for (int v = 0; v < V + 2; v++) for (auto &&e : G[v]) e.resCap = e.cap;
      V += 2; bnd -= MF::getFlow(ss, tt); V -= 2;
    }
    G.pop_back(); G.pop_back();
    for (int v = 0; v < V; v++) {
      G[v].pop_back(); G[v].pop_back(); for (auto &&e : G[v]) e.cap += e.dem;
    }
    if (flowType > 0 && ret.first) ret.second += MF::getFlow(s, t);
    return ret;
  }
  pair<bool, FlowUnit> getMinFlow(int s, int t) {
    return getFeasibleFlow(s, t, -1);
  }
  pair<bool, FlowUnit> getMaxFlow(int s, int t) {
    return getFeasibleFlow(s, t, 1);
  }
};


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

// A sample edge struct for minimum cost maximum flow
// Flow of edge can be found with cap - resCap
// Template Arguments:
//   _FlowUnit: the type of the flow
//   _CostUnit: the type of the cost
// Constructor Arguments:
//   to: the vertex that this directed edge ends at
//   rev: the index in the adjacency list of vertex to of the reverse edge
//   cap: the initial capacity of this edge
//   cost: the cost of this edge
// Fields:
//   FlowUnit: the type of the flow
//   CostUnit: the type of the cost
//   to: the vertex that this directed edge ends at
//   rev: the index in the adjacency list of vertex to of the reverse edge
//   cap: the initial capacity of this edge
//   resCap: the residual capacity of this edge
//   cost: the cost of this edge
// Tested:
//   https://open.kattis.com/problems/mincostmaxflow
//   https://loj.ac/p/102
//   https://uoj.ac/problem/487
//   https://www.spoj.com/problems/SCITIES/
//   https://open.kattis.com/problems/workers
template <class _FlowUnit, class _CostUnit> struct FlowCostEdge {
  using FlowUnit = _FlowUnit; using CostUnit = _CostUnit;
  int to, rev; FlowUnit cap, resCap; CostUnit cost;
  FlowCostEdge(int to, int rev, FlowUnit cap, CostUnit cost)
      : to(to), rev(rev), cap(cap), resCap(cap), cost(cost) {}
};

// Computes the maximum flow with the minimum cost in a flow network using the
//   Push Relabel algorithm with the highest label selection rule and
//   gap relabelling heuristic
// Costs are scaled by a factor that is O(V)
// CostUnit must be integral
// Flow circulations are added to negative cost cycles
// Vertices are 0-indexed
// Template Arguments:
//   Edge: a generic edge class (such as FlowCostEdge)
//     Required Fields:
//       FlowUnit: the type of the flow
//       CostUnit: the type of the cost, must be integral
//       to: the vertex that this directed edge ends at
//       rev: the index in the adjacency list of vertex to of the reverse edge
//       resCap: the residual capacity of this edge
//       cost: the cost of this edge
//     Required Functions:
//       constructor(to, rev, cap, cost): initializes the edge to vertex to
//         with the reverse index rev, capacity cap, and cost cost
// Constructor Arguments:
//   V: the number of vertices in the flow network
//   FLOW_EPS: a value for the flow epsilon
//   COST_INF: a value for the cost infinity
//   COST_EPS: a value for the cost epsilon
// Fields:
//   V: the number of vertices in the flow network
//   FLOW_EPS: a value for the flow epsilon
//   COST_INF: a value for the cost infinity
//   COST_EPS: a value for the cost epsilon
//   G: an adjacency list of all edges and reverse edges in the flow network
// Functions:
//   addEdge(v, w, vwCap, vwCost): adds an edge from v to w with capacity
//     vwCap and cost vwCost
//   getFlowMinCost(s, t): returns the maximum flow (or circulation if s and t
//     are -1) from s to t with the minimum cost
// In practice, has a very small constant
// Time Complexity:
//   constructor: O(V)
//   addEdge: O(1)
//   getFlowMinCost: O(E V^2 log (VC)) where C is the maximum edge cost,
//     much faster in practice
// Memory Complexity: O(V + E)
// Tested:
//   https://open.kattis.com/problems/mincostmaxflow
//   https://loj.ac/p/102
//   https://uoj.ac/problem/487
//   https://www.spoj.com/problems/SCITIES/
//   https://open.kattis.com/problems/workers
template <class Edge>
struct PushRelabelMinCostMaxFlow : public PushRelabelMaxFlow<Edge> {
  using MF = PushRelabelMaxFlow<Edge>;
  using FlowUnit = typename Edge::FlowUnit;
  using CostUnit = typename Edge::CostUnit;
  static_assert(is_integral<CostUnit>::value, "CostUnit must be integral");
  using MF::V; using MF::G;
  using MF::FLOW_EPS; CostUnit COST_INF, COST_EPS, negCost;
  PushRelabelMinCostMaxFlow(
      int V, FlowUnit FLOW_EPS = FlowUnit(1e-9),
      CostUnit COST_INF = numeric_limits<CostUnit>::max(),
      CostUnit COST_EPS = CostUnit(1e-9))
      : MF(V, FLOW_EPS), COST_INF(COST_INF), COST_EPS(COST_EPS),
        negCost(CostUnit()) {}
  void addEdge(int v, int w, FlowUnit vwCap, CostUnit vwCost) {
    if (v == w) {
      if (vwCost < CostUnit()) negCost += vwCap * vwCost;
      return;
    }
    G[v].emplace_back(w, int(G[w].size()), vwCap, vwCost);
    G[w].emplace_back(v, int(G[v].size()) - 1, FlowUnit(), -vwCost);
  }
  pair<FlowUnit, CostUnit> getFlowMinCost(int s = -1, int t = -1) {
    CostUnit minCost = CostUnit(), mul = CostUnit(2) << __lg(V);
    CostUnit bnd = CostUnit(); vector<int> stk(V), infs(V, 0); int top = 0;
    vector<CostUnit> phi(V, CostUnit()); vector<FlowUnit> ex(V, FlowUnit());
    auto costP = [&] (int v, const Edge &e) {
      int netInfs = infs[v] - infs[e.to]; if (netInfs > 0) return COST_INF;
      else if (netInfs < 0) return -COST_INF;
      else return e.cost + phi[v] - phi[e.to];
    };
    auto push = [&] (int v, Edge &e, FlowUnit df, bool pushToStack) {
      if (e.resCap < df) df = e.resCap;
      int w = e.to; e.resCap -= df; G[w][e.rev].resCap += df;
      ex[v] -= df; ex[w] += df;
      if (pushToStack && FLOW_EPS < ex[e.to] && ex[e.to] <= df + FLOW_EPS)
        stk[top++] = e.to;
    };
    auto relabel = [&] (int v, CostUnit delta) {
      if (delta < COST_INF) phi[v] -= delta + bnd;
      else { infs[v]--; phi[v] -= bnd; }
    };
    auto lookAhead = [&] (int v) {
      if (abs(ex[v]) > FLOW_EPS) return false;
      CostUnit delta = COST_INF; for (auto &&e : G[v]) {
        if (e.resCap <= FLOW_EPS) continue;
        CostUnit c = costP(v, e); if (c < -COST_EPS) return false;
        else delta = min(delta, c);
      }
      relabel(v, delta); return true;
    };
    auto discharge = [&] (int v) {
      CostUnit delta = COST_INF; for (int i = 0; i < int(G[v].size()); i++) {
        Edge &e = G[v][i]; if (e.resCap <= FLOW_EPS) continue;
        if (costP(v, e) < -COST_EPS) {
          if (lookAhead(e.to)) { i--; continue; }
          push(v, e, ex[v], true); if (abs(ex[v]) <= FLOW_EPS) return;
        } else delta = min(delta, costP(v, e));
      }
      relabel(v, delta); stk[top++] = v;
    };
    for (int v = 0; v < V; v++) for (auto &&e : G[v]) {
      minCost += e.cost * e.resCap; e.cost *= mul; bnd = max(bnd, e.cost);
    }
    FlowUnit maxFlow = (s == -1 || t == -1) ? FlowUnit() : MF::getFlow(s, t);
    while (bnd > 1) {
      bnd = max(bnd / CostUnit(8), CostUnit(1)); top = 0;
      for (int v = 0; v < V; v++) for (auto &&e : G[v])
        if (costP(v, e) < -COST_EPS && e.resCap > FLOW_EPS)
          push(v, e, e.resCap, false);
      for (int v = 0; v < V; v++) if (ex[v] > FLOW_EPS) stk[top++] = v;
      while (top > 0) discharge(stk[--top]);
    }
    for (int v = 0; v < V; v++)
      for (auto &&e : G[v]) { e.cost /= mul; minCost -= e.cost * e.resCap; }
    return make_pair(maxFlow, (minCost /= 2) += negCost);
  }
};


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

// Assigns colors to the edges of an undirected simple graph such that
//   no edges that share an endpoint have the same color
// If D is the maximum degree of any vertex, then at most D + 1 colors
//   will be used
// Vertices and colors are 0-indexed
// Constructor Arguments:
//   V: the number of vertices in the simple graph
//   edges: a vector of pairs in the form (v, w) representing
//     an undirected edge in the simple graph (no self loops or parallel edges)
//     between vertices v and w
// Fields:
//   color: a vector of integers in the range [0, D] that has a length equal
//     to the length of edges representing the coloring of the edges
// In practice, has a very small constant
// Time Complexity: O(VE)
// Memory Complexity: O(VD) where D is the maximum degree of any vertex
// Tested:
//   https://open.kattis.com/problems/gamescheduling
struct EdgeColoring {
  vector<int> color;
  EdgeColoring(int V, const vector<pair<int, int>> &edges)
      : color(edges.size(), 0) {
    vector<int> deg(V, 0), fan(V, 0), avail(V, 0);
    for (auto &&e : edges) { deg[e.first]++; deg[e.second]++; }
    int D = *max_element(deg.begin(), deg.end());
    vector<vector<int>> adj(V, vector<int>(D + 1, -1));
    vector<int> loc(D + 1, 0); for (auto &&e : edges) {
      int v, w; tie(v, w) = e; fan[0] = w; fill(loc.begin(), loc.end(), 0);
      int at = v, en = v, d = avail[w], c = avail[v], ind = 0, i = 0;
      for (; !loc[d] && (w = adj[v][d]) != -1; d = avail[w]) {
        deg[loc[d] = ++ind] = d; fan[ind] = w;
      }
      for (int cd = d; at != -1; cd ^= c ^ d, at = adj[at][cd])
        swap(adj[at][cd], adj[en = at][cd ^ c ^ d]);
      deg[loc[d]] = c; while (adj[fan[i]][d] != -1) {
        int l = fan[i], r = fan[++i], e = deg[i];
        adj[v][e] = l; adj[l][e] = v; adj[r][e] = -1; avail[r] = e;
      }
      adj[v][d] = fan[i]; adj[fan[i]][d] = v; for (int y : {fan[0], v, en})
        for (int &z = avail[y] = 0; adj[y][z] != -1; z++);
    }
    for (int i = 0; i < int(edges.size()); i++)
      while (adj[edges[i].first][color[i]] != edges[i].second) color[i]++;
  }
};


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

// Computes the minimum number of colors required to color a simple graph such
//   that no two adjacent vertices share the same color
// Vertices are 0-indexed
// Function Arguments:
//   V: the number of vertices in the simple graph
//   edges: a vector of pairs in the form (v, w) representing
//     an undirected edge in the simple graph (no self loops or parallel edges)
//     between vertices v and w
// Return Value: the chromatic number (minimum number of colors)
// In practice, has a small constant
// Time Complexity: O(V 2^V)
// Memory Complexity: O(2^V)
// Tested:
//   https://judge.yosupo.jp/problem/chromatic_number
int chromaticNumber(int V, const vector<pair<int, int>> &edges) {
  int ans = V, PV = 1 << V; vector<int> mask(V, 0), ind(PV), aux(PV);
  for (auto &&e : edges) {
    mask[e.first] |= 1 << e.second; mask[e.second] |= 1 << e.first;
  }
  for (int d : {7, 9, 21, 33, 87, 93, 97}) {
    long long mod = 1e9 + d; fill(ind.begin(), ind.end(), 0); ind[0] = 1;
    fill(aux.begin(), aux.end(), 1); for (int i = 0; i < PV; i++) {
      int v = __builtin_ctz(i);
      ind[i] = ind[i ^ (1 << v)] + ind[(i ^ (1 << v)) & ~mask[v]];
    }
    for (int k = 1; k < ans; k++) {
      long long chi = 0; for (int i = 0; i < PV; i++) {
        int j = i ^ (i >> 1); aux[j] = (long long) aux[j] * ind[j] % mod;
        chi += (i & 1) ? aux[j] : -aux[j];
      }
      if (chi % mod != 0) ans = k;
    }
  }
  return ans;
}


text:
#pragma once
#include <bits/stdc++.h>
#include "../../datastructures/unionfind/UnionFind.h"
#include "../matching/HopcroftKarpMaxMatch.h"
using namespace std;

// Assigns colors to the edges of a bipartite graph such that
//   no edges that share an endpoint have the same color
// If D is the maximum degree of any vertex, then exactly D colors
//   will be used
// Vertices and colors are 0-indexed
// Constructor Arguments:
//   V: the number of vertices in the simple graph
//   edges: a vector of pairs in the form (v, w) representing
//     an undirected edge in the graph; side[v] != side[w] must hold
//   side: the side of the bipartition each vertex is part of
// Fields:
//   color: a vector of integers in the range [0, D) that has a length equal
//     to the length of edges representing the coloring of the edges
// In practice, has a small constant
// Time Complexity: O(V log V + E log D sqrt (E / D)) where D is the maximum
//   degree of any vertex
// Memory Complexity: O(V + E)
// Tested:
//   https://judge.yosupo.jp/problem/bipartite_edge_coloring
struct BipartiteEdgeColoring {
  struct Pair {
    int s, v; Pair(int s, int v) : s(s), v(v) {}
    bool operator < (const Pair &o) const { return s > o.s; }
  };
  int makeDRegular(int &V, vector<pair<int, int>> &edges, vector<bool> &side) {
    vector<int> deg(V, 0); for (auto &&e : edges) {
      if (side[e.first]) swap(e.first, e.second);
      deg[e.first]++; deg[e.second]++;
    }
    int D = *max_element(deg.begin(), deg.end()); UnionFind uf(V);
    for (int s = 0; s < 2; s++) {
      std::priority_queue<Pair> PQ;
      for (int v = 0; v < V; v++) if (side[v] == s) PQ.emplace(deg[v], v);
      while (int(PQ.size()) >= 2) {
        Pair a = PQ.top(); PQ.pop(); Pair b = PQ.top(); PQ.pop();
        if (a.s + b.s <= D) { uf.join(a.v, b.v); PQ.emplace(a.s + b.s, a.v); }
      }
    }
    vector<int> cnt(2, 0), id(V, -1); int curId = 0;
    for (int s = 0; s < 2; s++) for (int v = 0; v < V; v++)
      if (uf.find(v) == v && side[v] == s) { id[v] = curId++; cnt[s]++; }
    deg.assign(V = max(cnt[0], cnt[1]) * 2, 0); edges.reserve(V * D / 2);
    side.reserve(V); side.assign(cnt[0] + cnt[1], true);
    fill(side.begin(), side.begin() + cnt[0], false);
    for (int s = 0; s < 2; s++) for (; cnt[s] * 2 < V; cnt[s]++)
      side.push_back(s);
    for (auto &&e : edges) {
      deg[e.first = id[uf.find(e.first)]]++;
      deg[e.second = id[uf.find(e.second)]]++;
    }
    for (int v = 0, w = 0; v < V; v++) while (!side[v] && deg[v] < D) {
      while (!side[w] || deg[w] == D) w++;
      edges.emplace_back(v, w); deg[v]++; deg[w]++;
    }
    return D;
  }
  vector<int> eulerianCircuit(
      int V, const vector<pair<int, int>> &edges, const vector<int> &inds) {
    vector<vector<pair<int, int>>> G(V); vector<int> circuit;
    for (int i = 0; i < int(inds.size()); i++) {
      int v, w; tie(v, w) = edges[inds[i]];
      G[v].emplace_back(w, i); G[w].emplace_back(v, i);
    }
    vector<bool> vis1(V, false), vis2(inds.size(), false);
    vector<pair<int, int>> stk; for (int s = 0; s < V; s++) if (!vis1[s]) {
      stk.clear(); stk.emplace_back(s, -1); while (!stk.empty()) {
        int v, w, e; tie(v, e) = stk.back(); vis1[v] = true;
        if (G[v].empty()) { circuit.emplace_back(e); stk.pop_back(); }
        else {
          tie(w, e) = G[v].back(); G[v].pop_back();
          if (!vis2[e]) { vis2[e] = true; stk.emplace_back(w, e); }
        }
      }
      circuit.pop_back();
    }
    for (auto &&e : circuit) e = inds[e];
    return circuit;
  }
  vector<int> color;
  BipartiteEdgeColoring(int V, vector<pair<int, int>> edges, vector<bool> side)
      : color(edges.size(), 0) {
    for (auto &&e : edges) assert(side[e.first] != side[e.second]);
    int D = makeDRegular(V, edges, side), curCol = 0;
    function<void(int, const vector<int> &)> rec = [&] (
        int d, const vector<int> &inds) {
      if (d == 0) return;
      else if (d == 1) {
        for (int e : inds) if (e < int(color.size())) color[e] = curCol;
        curCol++;
      } else if (d % 2 == 0) {
        vector<int> circuit = eulerianCircuit(V, edges, inds), half1, half2;
        half1.reserve(circuit.size() / 2); half2.reserve(circuit.size() / 2);
        for (int i = 0; i < int(circuit.size()); i += 2) {
          half1.push_back(circuit[i]); half2.push_back(circuit[i + 1]);
        }
        rec(d / 2, half1); rec(d / 2, half2);
      } else {
        vector<vector<int>> G(V); for (int e : inds) {
          int v, w; tie(v, w) = edges[e]; G[v].push_back(w); G[w].push_back(v);
        }
        vector<int> unmatched; HopcroftKarpMaxMatch mm(G, side);
        for (int e : inds) {
          int v, w; tie(v, w) = edges[e]; if (mm.mate[v] == w) {
            mm.mate[v] = -1; if (e < int(color.size())) color[e] = curCol;
          } else unmatched.push_back(e);
        }
        curCol++; rec(d - 1, unmatched);
      }
    };
    vector<int> inds(edges.size()); iota(inds.begin(), inds.end(), 0);
    rec(D, inds);
  }
};


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

// Link Cut Tree supporting path operations on a dynamic tree,
//   backed by a splay tree
// Vertices are 0-indexed
// Template Arguments:
//   Node: a generic node class (sample structs are in BSTNode)
//     Required Fields:
//       Data: the data type
//       Lazy: the lazy type
//       static const RANGE_UPDATES: a boolean indicating whether
//         range updates are permitted
//       static const RANGE_QUERIES: a boolean indicating whether
//         range queries are permitted
//       static const RANGE_REVERSALS: a boolean indicating whether
//         range reversals are permitted
//       sz: the number of nodes in the subtree
//       l: a pointer to the left child
//       r: a pointer to the right child
//       p: a pointer to the parent
//       val: only required if getFirst is called, the value being stored
//       sbtr: only required if RANGE_QUERIES is true to support path queries,
//         the aggregate value of type Data for the subtree
//     Required Functions:
//       constructor(v): initializes a node with the value v
//       update(): updates the current node's information based on its children
//       propagate(): propagates the current node's lazy information (including
//         rev) to its children
//       apply(v): applies the lazy value v to the node
//       reverse(): only required if RANGE_REVERSALS is true to support
//         rerooting, reverse this node's subtree (aggregate data
//         and any lazy flags should be reversed)
//       static qdef(): returns the query default value
// Constructor Arguments:
//   A: a vector of type Node::Data
// Functions:
//   makeRoot(x): only valid if Node::RANGE_REVERSALS is true, makes x the
//     root of its connected component
//   lca(x, y): returns the lowest common ancestor of x and y in
//     the current forest, returns -1 if not connected
//   connected(x, y): returns whether x and y are connected
//   link(x, y): only valid if Node::RANGE_REVERSALS is true, links the nodes
//     x and y, assumes x and y are not connected, reroots the tree at node y
//   safeLink(x, y): only valid if Node::RANGE_REVERSALS is true, links
//     the nodes x and y, returns false if already connected, true otherwise,
//     reroots the tree at node y
//   linkParent(par, ch): makes par the parent of the node ch, returns false if
//     par and ch are already connected, true otherwise
//   cut(x, y): only valid if Node::RANGE_REVERSALS is true, cuts the edge
//     between nodes x and y, returns false if this edge doesn't exist,
//     true otherwise, reroots the forest at node x
//   cutParent(x): cuts the edge between node x and its parent, returns false
//     if no parent exists (x is a root), true otherwise
//   findParent(x): returns the parent of node x, -1 if it doesn't exist
//   findRoot(x): returns the root of the forest containing node x
//   depth(x): returns the depth of node x, where the depth of the root is 0
//   kthParent(x): returns the kth parent of node x (0th parent is x,
//     1st parent is the parent of x), -1 if it doesn't exist
//   updateVertex(x, v): updates the node x with the lazy value v
//   updatePathFromRoot(to, v): only valid if Node::RANGE_UPDATES is true,
//     updates the path from the root of the forest containing
//     node to, to node to, with the lazy value v
//   updatePath(from, to, v): only valid if Node::RANGE_UPDATES and
//     Node::RANGE_REVERSALS are true, updates the path from node
//     from to node to, reroots the forest at node from, with the lazy value v,
//     reroots the forest at node from
//   queryVertex(x): returns the value of node x
//   queryPathFromRoot(to): only valid if Node::RANGE_QUERIES is true,
//     returns the aggregate value of the path from the root of the forest
//     containing node to, to node to
//   queryPath(from, to): only valid if Node::RANGE_QUERIES and
//     Node::RANGE_REVERSALS are true, returns the aggregate value of the path
//     from node from to node to, reroots the forest at node from, reroots the
//     forest at node from
// In practice, has a moderate constant
// Time Complexity:
//   constructor: O(N)
//   makeRoot, lca, connected, link, safeLink, linkParent, cut, cutParent,
//     findParent, findRoot, depth, kthParent, updateVertex,
//     updatePathFromRoot, updatePath, queryVertex, queryPathFromRoot,
//     queryPath: O(log N) amortized
// Memory Complexity: O(N)
// Tested:
//   https://dmoj.ca/problem/coi08p2
//   https://codeforces.com/contest/13/problem/E
//   https://dmoj.ca/problem/ds5easy
//   https://www.spoj.com/problems/QTREE2/
//   https://judge.yosupo.jp/problem/dynamic_tree_vertex_set_path_composite
//   https://oj.uz/problem/view/JOI13_synchronization
template <class Node> struct LCT : public Splay<Node, vector<Node>> {
  using Tree = Splay<Node, vector<Node>>;
  using Tree::TR; using Tree::makeNode; using Tree::splay; using Tree::select;
  using Data = typename Node::Data; using Lazy = typename Node::Lazy;
  int vert(Node *x) { return x - TR.data(); }
  Node *access(Node *x) {
    Node *last = nullptr;
    for (Node *y = x; y; y = y->p) { splay(y); y->r = last; last = y; }
    splay(x); return last;
  }
  template <const int _ = Node::RANGE_REVERSALS>
  typename enable_if<_>::type makeRoot(Node *x) { access(x); x->reverse(); }
  Node *findMin(Node *x) {
    for (x->propagate(); x->l; (x = x->l)->propagate());
    splay(x); return x;
  }
  Node *findMax(Node *x) {
    for (x->propagate(); x->r; (x = x->r)->propagate());
    splay(x); return x;
  }
  template <const int _ = Node::RANGE_REVERSALS>
  typename enable_if<_>::type makeRoot(int x) { makeRoot(&TR[x]); }
  int lca(int x, int y) {
    if (x == y) return x;
    access(&TR[x]); Node *ny = access(&TR[y]); return TR[x].p ? vert(ny) : -1;
  }
  bool connected(int x, int y) { return lca(x, y) != -1; }
  template <const int _ = Node::RANGE_REVERSALS>
  typename enable_if<_>::type link(int x, int y) {
    makeRoot(y); TR[y].p = &TR[x];
  }
  template <const int _ = Node::RANGE_REVERSALS>
  typename enable_if<_, bool>::type safeLink(int x, int y) {
    if (connected(x, y)) return false;
    link(x, y); return true;
  }
  bool linkParent(int par, int ch) {
    access(&TR[ch]); if (TR[ch].l) return false;
    TR[ch].p = &TR[par]; return true;
  }
  template <const int _ = Node::RANGE_REVERSALS>
  typename enable_if<_, bool>::type cut(int x, int y) {
    makeRoot(x); access(&TR[y]);
    if (&TR[x] != TR[y].l || TR[x].r) return false;
    TR[y].l->p = nullptr; TR[y].l = nullptr; return true;
  }
  bool cutParent(int x) {
    access(&TR[x]); if (!TR[x].l) return false;
    TR[x].l->p = nullptr; TR[x].l = nullptr; return true;
  }
  int findParent(int x) {
    access(&TR[x]); return TR[x].l ? vert(findMax(TR[x].l)) : -1;
  }
  int findRoot(int x) { access(&TR[x]); return vert(findMin(&TR[x])); }
  int depth(int x) { access(&TR[x]); return TR[x].l ? TR[x].l->sz : 0; }
  int kthParent(int x, int k) {
    int d = depth(x); Node *nx = &TR[x];
    return k <= d ? vert(select(nx, d - k)) : -1;
  }
  void updateVertex(int x, const Lazy &v) {
    access(&TR[x]); Node *l = TR[x].l; TR[x].l = nullptr;
    TR[x].apply(v); TR[x].propagate(); TR[x].update(); TR[x].l = l;
  }
  template <const int _ = Node::RANGE_UPDATES>
  typename enable_if<_>::type updatePathFromRoot(int to, const Lazy &v) {
    access(&TR[to]); TR[to].apply(v);
  }
  template <const int _ = Node::RANGE_UPDATES && Node::RANGE_REVERSALS>
  typename enable_if<_, bool>::type updatePath(
      int from, int to, const Lazy &v) {
    makeRoot(from); access(&TR[to]);
    if (from != to && !TR[from].p) return false;
    TR[to].apply(v); return true;
  }
  Data queryVertex(int x) { access(&TR[x]); return TR[x].val; }
  template <const int _ = Node::RANGE_QUERIES>
  typename enable_if<_, Data>::type queryPathFromRoot(int to) {
    access(&TR[to]); return TR[to].sbtr;
  }
  template <const int _ = Node::RANGE_QUERIES && Node::RANGE_REVERSALS>
  typename enable_if<_, Data>::type queryPath(int from, int to) {
    makeRoot(from); access(&TR[to]);
    return from == to || TR[from].p ? TR[to].sbtr : Node::qdef();
  }
  LCT(const vector<Data> &A) {
    TR.reserve(A.size()); for (auto &&a : A) makeNode(a);
  }
};


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

// Top Tree supporting path and subtree operations on a dynamic tree
// Vertices are 0-indexed
// Template Arguments:
//   C: struct to combine data and lazy values
//     Required Fields:
//       Data: the data type
//       Lazy: the lazy type
//     Required Functions:
//       static qdef(): returns the query default value of type Data
//       static merge(l, r): returns the values l of type Data merged with
//         r of type Data, must be associative
//       static applyLazy(l, r, k): returns the value r of type Lazy applied to
//         l of type Data over a segment of length k, must be associative
//       static mergeLazy(l, r): returns the values l of type Lazy merged with
//         r of type Lazy, must be associative
//       static revData(v): reverses the value v of type Data
//     Sample Struct: supporting range assignments and range sum queries
//       struct C {
//         using Data = int;
//         using Lazy = int;
//         static Data qdef() { return 0; }
//         static Lazy ldef() { return numeric_limits<int>::min(); }
//         static Data merge(const Data &l, const Data &r) { return l + r; }
//         static Data applyLazy(const Data &l, const Lazy &r, int k) {
//           return r * k;
//         }
//         static Lazy mergeLazy(const Lazy &l, const Lazy &r) { return r; }
//         static void revData(Data &v) {}
//       };
// Constructor Arguments:
//   A: a vector of type C::Data
// Functions:
//   makeRoot(x): makes x the root of its connected component
//   lca(x, y): returns the lowest common ancestor of x and y in
//     the current forest, returns -1 if not connected
//   connected(x, y): returns whether x and y are connected
//   link(par, ch): makes par the parent of the node ch, assumes par and ch
//     are not connected, reroots the tree at node par
//   safeLink(par, ch): makes par the parent of the node ch, returns false if
//     already connected, true otherwise, reroots the tree at node par
//   cutParent(x): cuts the edge between node x and its parent, returns false
//     if no parent exists (x is a root), true otherwise
//   findParent(x): returns the parent of node x, -1 if it doesn't exist
//   findRoot(x): returns the root of the forest containing node x
//   depth(x): returns the depth of node x, where the depth of the root is 0
//   kthParent(x): returns the kth parent of node x (0th parent is x,
//     1st parent is the parent of x), -1 if it doesn't exist
//   updateVertex(x, v): updates the node x with the lazy value v
//   updatePathFromRoot(to, v): updates the path from the root of the forest
//     containing node to, to node to, with the lazy value v
//   updatePath(from, to, v): updates the path from node
//     from to node to, reroots the forest at node from, with the lazy value v,
//     reroots the forest at node from
//   updateSubtree(x): updates the subtree of node x with the lazy value v
//   queryVertex(x): returns the value of node x
//   queryPathFromRoot(to): returns the aggregate value of the path from
//     the root of the forest containing node to, to node to
//   queryPath(from, to): returns the aggregate value of the path
//     from node from to node to, reroots the forest at node from, reroots the
//     forest at node from
//   querySubtree(x): returns the aggregate value of the subtree of node x
// In practice, has a large constant
// Time Complexity:
//   constructor: O(N)
//   makeRoot, lca, connected, link, safeLink, cutParent,
//     findParent, findRoot, depth, kthParent, updateVertex,
//     updatePathFromRoot, updatePath, updateSubtree, queryVertex,
//     queryPathFromRoot, queryPath, querySubtree: O(log N) amortized
// Memory Complexity: O(N)
// Tested:
//   https://www.spoj.com/problems/QTREE2/
//   https://judge.yosupo.jp/problem/dynamic_tree_vertex_set_path_composite
//   https://judge.yosupo.jp/problem/dynamic_tree_subtree_add_subtree_sum
//   https://dmoj.ca/problem/ds5
template <class C> struct TopTree {
  using Data = typename C::Data; using Lazy = typename C::Lazy;
  struct Node {
    bool rev, aux; int szpath, szvtr; array<Node *, 4> ch; Node *p;
    Lazy lzpath, lzvtr; Data val, path, vtr;
    Node(bool aux, const Data &v)
        : rev(false), aux(aux), szpath(aux ? 0 : 1), szvtr(0),
          p(nullptr), lzpath(C::ldef()), lzvtr(C::ldef()), val(v),
          path(aux ? C::qdef() : v), vtr(C::qdef()) {
      ch.fill(nullptr);
    }
    void applyVal(const Lazy &v) { val = C::applyLazy(val, v, 1); }
    void applyPath(const Lazy &v) {
      applyVal(v); lzpath = C::mergeLazy(lzpath, v);
      if (szpath > 0) path = C::applyLazy(path, v, szpath);
    }
    void applyVtr(const Lazy &v, bool ap = true) {
      lzvtr = C::mergeLazy(lzvtr, v);
      if (szvtr > 0) vtr = C::applyLazy(vtr, v, szvtr);
      if (!aux && ap) applyPath(v);
    }
    void update() {
      szpath = aux ? 0 : 1; path = aux ? C::qdef() : val;
      szvtr = 0; vtr = C::qdef(); for (int i = 0; i < 4; i++) if (ch[i]) {
        vtr = C::merge(vtr, ch[i]->vtr); szvtr += ch[i]->szvtr; if (i < 2) {
          path = i ? C::merge(path, ch[i]->path) : C::merge(ch[i]->path, path);
          szpath += ch[i]->szpath;
        } else { vtr = C::merge(vtr, ch[i]->path); szvtr += ch[i]->szpath; }
      }
    }
    void propagate() {
      if (rev) {
        for (int i = 0; i < 2; i++) if (ch[i]) ch[i]->reverse();
        rev = false;
      }
      if (lzpath != C::ldef() && !aux) {
        for (int i = 0; i < 2; i++) if (ch[i]) ch[i]->applyPath(lzpath);
        lzpath = C::ldef();
      }
      if (lzvtr != C::ldef()) {
        for (int i = 0; i < 4; i++) if (ch[i]) ch[i]->applyVtr(lzvtr, i >= 2);
        lzvtr = C::ldef();
      }
    }
    void reverse() { rev = !rev; swap(ch[0], ch[1]); C::revData(path); }
  };
  vector<Node> TR; vector<Node *> deleted, stk;
  Node *makeNode(bool aux, const Data &v) {
    if (deleted.empty()) { TR.emplace_back(aux, v); return &TR.back(); }
    Node *x = deleted.back(); deleted.pop_back(); *x = Node(aux, v); return x;
  }
  int vert(Node *x) { return x - TR.data(); }
  bool isRoot(Node *x, bool t) {
    if (t) return !x->p || !x->aux || !x->p->aux;
    else return !x->p || (x != x->p->ch[0] && x != x->p->ch[1]);
  }
  void connect(Node *x, Node *p, int i) {
    if (x) x->p = p;
    if (i != -1) p->ch[i] = x;
  }
  int findInd(Node *x) {
    for (int i = 0; i < 4; i++) if (x->p->ch[i] == x) return i;
    return -1;
  }
  void rotate(Node *x, int t) {
    Node *p = x->p, *g = p->p; bool isL = x == p->ch[t];
    if (g) connect(x, g, findInd(p));
    else x->p = nullptr;
    connect(x->ch[t ^ isL], p, t ^ !isL); connect(p, x, t ^ isL);
    p->update();
  }
  void splay(Node *x, int t) {
    while (!isRoot(x, t)) {
      Node *p = x->p, *g = p->p;
      if (!isRoot(p, t)) rotate((x == p->ch[t]) == (p == g->ch[t]) ? p : x, t);
      rotate(x, t);
    }
    x->update();
  }
  Node *select(Node *&root, int k) {
    Node *last = nullptr; while (root) {
      (last = root)->propagate();
      int t = root->ch[0] ? root->ch[0]->szpath : 0;
      if (t > k) root = root->ch[0];
      else if (t < k) { root = root->ch[1]; k -= t + 1; }
      else break;
    }
    if (last) splay(root = last, 0);
    return root;
  }
  void add(Node *x, Node *y) {
    Node *z = makeNode(true, C::qdef());
    connect(y->ch[2], z, 2); connect(x, z, 3); connect(z, y, 2); z->update();
  }
  void rem(Node *x) {
    Node *p = x->p, *g = p->p; connect(p->ch[findInd(x) ^ 1], g, findInd(p));
    splay(g, 2); deleted.push_back(x->p); x->p = nullptr;
  }
  void touch(Node *x) {
    for (Node *y = x; y->p; y = y->p) stk.push_back(y->p);
    for (; !stk.empty(); stk.pop_back()) stk.back()->propagate();
    x->propagate();
  }
  Node *nextChain(Node *x) {
    splay(x, 0); if (!x->p) return nullptr;
    Node *p = x->p; splay(p, 2); return p->p;
  }
  Node *access(Node *x) {
    touch(x); Node *last = nullptr; for (Node *y = x; y; y = nextChain(y)) {
      splay(y, 0); if (last) rem(last);
      if (y->ch[1]) add(y->ch[1], y);
      connect(last, y, 1); last = y;
    }
    splay(x, 0); return last;
  }
  void makeRoot(Node *x) { access(x); x->reverse(); }
  Node *findMin(Node *x) {
    for (x->propagate(); x->ch[0]; (x = x->ch[0])->propagate());
    splay(x, 0); return x;
  }
  Node *findMax(Node *x) {
    for (x->propagate(); x->ch[1]; (x = x->ch[1])->propagate());
    splay(x, 0); return x;
  }
  void makeRoot(int x) { makeRoot(&TR[x]); }
  int lca(int x, int y) {
    if (x == y) return x;
    access(&TR[x]); Node *ny = access(&TR[y]); return TR[x].p ? vert(ny) : -1;
  }
  bool connected(int x, int y) { return lca(x, y) != -1; }
  void link(int par, int ch) {
    makeRoot(par); access(&TR[ch]); add(&TR[par], &TR[ch]);
  }
  bool safeLink(int par, int ch) {
    if (connected(par, ch)) return false;
    link(par, ch); return true;
  }
  bool cutParent(int x) {
    access(&TR[x]); if (!TR[x].ch[0]) return false;
    TR[x].ch[0]->p = nullptr; TR[x].ch[0] = nullptr; return true;
  }
  int findParent(int x) {
    access(&TR[x]); return TR[x].ch[0] ? vert(findMax(TR[x].ch[0])) : -1;
  }
  int findRoot(int x) { access(&TR[x]); return vert(findMin(&TR[x])); }
  int depth(int x) {
    access(&TR[x]); return TR[x].ch[0] ? TR[x].ch[0]->szpath : 0;
  }
  int kthParent(int x, int k) {
    int d = depth(x); Node *nx = &TR[x];
    return k <= d ? vert(select(nx, d - k)) : -1;
  }
  void updateVertex(int x, const Lazy &v) {
    access(&TR[x]); TR[x].applyVal(v); TR[x].update();
  }
  void updatePathFromRoot(int to, const Lazy &v) {
    access(&TR[to]); TR[to].applyPath(v);
  }
  bool updatePath(int from, int to, const Lazy &v) {
    makeRoot(from); access(&TR[to]);
    if (from != to && !TR[from].p) return false;
    TR[to].applyPath(v); return true;
  }
  void updateSubtree(int x, const Lazy &v) {
    access(&TR[x]); TR[x].applyVal(v);
    for (int i = 2; i < 4; i++) if (TR[x].ch[i]) TR[x].ch[i]->applyVtr(v);
  }
  Data queryVertex(int x) { access(&TR[x]); return TR[x].val; }
  Data queryPathFromRoot(int to) { access(&TR[to]); return TR[to].sbtr; }
  Data queryPath(int from, int to) {
    makeRoot(from); access(&TR[to]);
    return from == to || TR[from].p ? TR[to].path : C::qdef();
  }
  Data querySubtree(int x) {
    access(&TR[x]); Data ret = TR[x].val;
    for (int i = 2; i < 4; i++) if (TR[x].ch[i]) {
      ret = C::merge(C::merge(ret, TR[x].ch[i]->vtr), TR[x].ch[i]->path);
    }
    return ret;
  }
  TopTree(const vector<Data> &A) {
    TR.reserve(A.size() * 2); for (auto &&a : A) makeNode(false, a);
  }
};