// Flujo máximo. Las aristas se guardan en pares: la arista original tiene un número par e y su
// inversa es e ^ 1. Lo único que se guarda de cada una es su capacidad residual.
#pragma once

#include <algorithm>
#include <cstddef>
#include <limits>
#include <queue>
#include <vector>

class FlowNetwork {
public:
    explicit FlowNetwork(int nodes) : arcs_from_(static_cast<std::size_t>(nodes)) {}

    // Agrega from → to y, junto con ella, su inversa con capacidad 0. Devuelve el número de la
    // arista, que siempre es par. Dos aristas antiparalelas son dos pares independientes.
    int add_edge(int from, int to, long long capacity) {
        const int edge = static_cast<int>(to_.size());
        to_.push_back(to);
        residual_.push_back(capacity);
        arcs_from_[static_cast<std::size_t>(from)].push_back(edge);
        to_.push_back(from);
        residual_.push_back(0);
        arcs_from_[static_cast<std::size_t>(to)].push_back(edge + 1);
        return edge;
    }

    int node_count() const { return static_cast<int>(arcs_from_.size()); }
    int edge_count() const { return static_cast<int>(to_.size()) / 2; }
    int head(int arc) const { return to_[static_cast<std::size_t>(arc)]; }
    int tail(int arc) const { return to_[static_cast<std::size_t>(arc ^ 1)]; }
    long long residual(int arc) const { return residual_[static_cast<std::size_t>(arc)]; }
    const std::vector<int>& arcs_from(int node) const { return arcs_from_[static_cast<std::size_t>(node)]; }

    // De una arista original: lo que pasa por ella es lo que su inversa puede devolver.
    long long flow_on(int edge) const { return residual(edge ^ 1); }
    long long capacity(int edge) const { return residual(edge) + residual(edge ^ 1); }

    // Valor del flujo: lo que sale de source menos lo que entra.
    long long flow_value(int source) const {
        long long net = 0;
        for (int arc : arcs_from(source)) {
            if (arc % 2 == 0) net += flow_on(arc);
            else net -= flow_on(arc ^ 1);
        }
        return net;
    }

    // Empuja amount por cada arco del camino. Empujar por un arco es quitarle a él y darle a su par.
    void push(const std::vector<int>& path, long long amount) {
        for (int arc : path) {
            residual_[static_cast<std::size_t>(arc)] -= amount;
            residual_[static_cast<std::size_t>(arc ^ 1)] += amount;
        }
    }

    long long bottleneck(const std::vector<int>& path) const {
        long long smallest = std::numeric_limits<long long>::max();
        for (int arc : path) smallest = std::min(smallest, residual(arc));
        return smallest;
    }

private:
    std::vector<int> to_;
    std::vector<long long> residual_;
    std::vector<std::vector<int>> arcs_from_;
};

// Distancias en aristas desde source en la red residual. -1 si el nodo no se alcanza.
inline std::vector<int> residual_distances(const FlowNetwork& net, int source) {
    std::vector<int> dist(static_cast<std::size_t>(net.node_count()), -1);
    std::queue<int> pending;
    dist[static_cast<std::size_t>(source)] = 0;
    pending.push(source);
    while (!pending.empty()) {
        const int u = pending.front();
        pending.pop();
        for (int arc : net.arcs_from(u)) {
            const int v = net.head(arc);
            if (net.residual(arc) > 0 && dist[static_cast<std::size_t>(v)] < 0) {
                dist[static_cast<std::size_t>(v)] = dist[static_cast<std::size_t>(u)] + 1;
                pending.push(v);
            }
        }
    }
    return dist;
}

// Regla de Edmonds-Karp: el camino de aumento con menos aristas, con BFS. Vacío si no hay.
inline std::vector<int> shortest_path_rule(const FlowNetwork& net, int source, int sink) {
    std::vector<int> arc_into(static_cast<std::size_t>(net.node_count()), -1);
    std::queue<int> pending;
    pending.push(source);
    while (!pending.empty() && arc_into[static_cast<std::size_t>(sink)] < 0) {
        const int u = pending.front();
        pending.pop();
        for (int arc : net.arcs_from(u)) {
            const int v = net.head(arc);
            if (v != source && net.residual(arc) > 0 && arc_into[static_cast<std::size_t>(v)] < 0) {
                arc_into[static_cast<std::size_t>(v)] = arc;
                pending.push(v);
            }
        }
    }
    std::vector<int> path;
    if (arc_into[static_cast<std::size_t>(sink)] < 0) return path;
    for (int v = sink; v != source; v = net.tail(arc_into[static_cast<std::size_t>(v)])) {
        path.push_back(arc_into[static_cast<std::size_t>(v)]);
    }
    std::reverse(path.begin(), path.end());
    return path;
}

// Regla de DFS: el primer camino que encuentra una búsqueda en profundidad. Vacío si no hay.
inline std::vector<int> depth_first_rule(const FlowNetwork& net, int source, int sink) {
    std::vector<char> seen(static_cast<std::size_t>(net.node_count()), 0);
    std::vector<int> path;
    // Pila explícita: por cada nodo del camino, cuántos de sus arcos ya se probaron.
    std::vector<std::size_t> next_arc = {0};
    int u = source;
    seen[static_cast<std::size_t>(source)] = 1;
    while (u != sink) {
        const std::vector<int>& arcs = net.arcs_from(u);
        if (next_arc.back() == arcs.size()) {
            next_arc.pop_back();
            if (path.empty()) return {};
            u = net.tail(path.back());
            path.pop_back();
            continue;
        }
        const int arc = arcs[next_arc.back()++];
        const int v = net.head(arc);
        if (net.residual(arc) > 0 && !seen[static_cast<std::size_t>(v)]) {
            seen[static_cast<std::size_t>(v)] = 1;
            path.push_back(arc);
            next_arc.push_back(0);
            u = v;
        }
    }
    return path;
}

struct Augmentation {
    std::vector<int> path;
    long long amount;
};

struct FlowRun {
    long long value;
    long long iterations;
};

// Ford-Fulkerson: mientras la regla encuentre un camino de aumento, empujar su cuello de botella.
// Si log no es nulo, guarda cada aumento.
template <typename Rule>
FlowRun ford_fulkerson(FlowNetwork& net, int source, int sink, Rule rule, std::vector<Augmentation>* log = nullptr) {
    long long iterations = 0;
    while (true) {
        const std::vector<int> path = rule(net, source, sink);
        if (path.empty()) break;
        const long long amount = net.bottleneck(path);
        net.push(path, amount);
        iterations++;
        if (log != nullptr) log->push_back({path, amount});
    }
    return {net.flow_value(source), iterations};
}

// Edmonds-Karp es Ford-Fulkerson con la regla del camino más corto.
inline FlowRun edmonds_karp(FlowNetwork& net, int source, int sink, std::vector<Augmentation>* log = nullptr) {
    return ford_fulkerson(net, source, sink, shortest_path_rule, log);
}

// Al terminar el algoritmo, estos nodos son el lado de la fuente de un corte mínimo.
inline std::vector<char> residual_side(const FlowNetwork& net, int source) {
    const std::vector<int> dist = residual_distances(net, source);
    std::vector<char> side(dist.size(), 0);
    for (std::size_t v = 0; v < dist.size(); v++) side[v] = dist[v] >= 0;
    return side;
}