// Mide las tres formas de calcular todos los pares. Con -O2, una corrida, en la máquina del momento.
#include <chrono>
#include <cstdio>
#include <functional>
#include <random>
#include <vector>

#include "all_pairs.h"

namespace {

double seconds_of(const std::function<void()>& work) {
    const auto start = std::chrono::steady_clock::now();
    work();
    return std::chrono::duration<double>(std::chrono::steady_clock::now() - start).count();
}

std::vector<Edge> random_sparse(int n, int per_node, std::mt19937& gen) {
    std::vector<Edge> edges;
    for (int i = 0; i < n * per_node; i++) {
        const int a = static_cast<int>(gen() % static_cast<unsigned>(n));
        const int b = static_cast<int>(gen() % static_cast<unsigned>(n));
        if (a != b) edges.push_back({a, b, static_cast<long long>(gen() % 1000) + 1});
    }
    return edges;
}

std::vector<Edge> complete(int n, std::mt19937& gen) {
    std::vector<Edge> edges;
    for (int a = 0; a < n; a++) {
        for (int b = 0; b < n; b++) {
            if (a != b) edges.push_back({a, b, static_cast<long long>(gen() % 1000) + 1});
        }
    }
    return edges;
}

// Una cadena con sus aristas al revés: el peor orden para una sola fila.
std::vector<Edge> reversed_chain(int n) {
    std::vector<Edge> edges;
    for (int a = n - 2; a >= 0; a--) edges.push_back({a, a + 1, 1});
    return edges;
}

void measure(const char* family, int n, const std::vector<Edge>& edges) {
    Matrix a;
    Matrix b;
    FloydResult c;
    const double t_dijkstra = seconds_of([&] { a = all_pairs_dijkstra(n, edges); });
    const double t_bellman = seconds_of([&] { b = all_pairs_bellman_ford(n, edges); });
    const double t_floyd = seconds_of([&] { c = floyd_warshall(n, edges); });
    std::printf("%-16s %5d %8zu %10.3f %10.3f %10.3f  %s\n", family, n, edges.size(), t_dijkstra, t_bellman, t_floyd,
                a == b && b == c.dist ? "iguales" : "DISTINTOS");
}

}  // namespace

int main() {
    std::mt19937 gen(2026);
    std::printf("%-16s %5s %8s %10s %10s %10s\n", "grafo", "V", "E", "V Dijkstra", "V B-F", "F-W");
    for (int n : {100, 200, 400}) measure("disperso, E=4V", n, random_sparse(n, 4, gen));
    for (int n : {100, 200, 400}) measure("completo", n, complete(n, gen));
    for (int n : {200, 400, 800}) measure("cadena al revés", n, reversed_chain(n));
}

// Caminos mínimos entre todos los pares: Floyd-Warshall bien escrito, con el ciclo de k por dentro,
// y correr un algoritmo de un origen V veces.
#pragma once

#include <cstddef>
#include <vector>

#include "shortest_paths.h"

using Matrix = std::vector<std::vector<long long>>;

// El estado inicial, d_0: sin nodos intermedios.
inline Matrix initial_matrix(int n, const std::vector<Edge>& edges) {
    Matrix d(idx(n), std::vector<long long>(idx(n), INF));
    for (int v = 0; v < n; v++) d[idx(v)][idx(v)] = 0;
    for (const Edge& e : edges) {
        if (e.weight < d[idx(e.from)][idx(e.to)]) d[idx(e.from)][idx(e.to)] = e.weight;
    }
    return d;
}

struct FloydResult {
    Matrix dist;
    bool negative_cycle = false;  // si es verdadero, dist quedó a medias y no significa nada
};

// Una vuelta: permite k como nodo intermedio. Convierte d_k en d_{k + 1}, en el mismo lugar.
inline void floyd_round(Matrix& d, int n, int k) {
    for (int i = 0; i < n; i++) {
        if (d[idx(i)][idx(k)] == INF) continue;
        for (int j = 0; j < n; j++) {
            if (d[idx(k)][idx(j)] == INF) continue;
            const long long through = d[idx(i)][idx(k)] + d[idx(k)][idx(j)];
            if (through < d[idx(i)][idx(j)]) d[idx(i)][idx(j)] = through;
        }
    }
}

// Antes de cada vuelta mira la diagonal: si d[k][k] < 0, hay un ciclo negativo y se detiene.
// Por eso los valores no crecen sin control.
inline FloydResult floyd_warshall(int n, const std::vector<Edge>& edges) {
    FloydResult result;
    result.dist = initial_matrix(n, edges);
    for (int k = 0; k < n; k++) {
        if (result.dist[idx(k)][idx(k)] < 0) {
            result.negative_cycle = true;
            return result;
        }
        floyd_round(result.dist, n, k);
    }
    return result;
}

// El error de esta semana: compila, termina y da números creíbles.
inline Matrix floyd_warshall_k_inside(int n, const std::vector<Edge>& edges) {
    Matrix d = initial_matrix(n, edges);
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            for (int k = 0; k < n; k++) {
                if (d[idx(i)][idx(k)] == INF || d[idx(k)][idx(j)] == INF) continue;
                const long long through = d[idx(i)][idx(k)] + d[idx(k)][idx(j)];
                if (through < d[idx(i)][idx(j)]) d[idx(i)][idx(j)] = through;
            }
        }
    }
    return d;
}

// Todos los pares con V corridas de un algoritmo de un origen. El de Bellman-Ford supone que
// no hay ciclos negativos, o devuelve la fila tal como quedó.
inline Matrix all_pairs_dijkstra(int n, const std::vector<Edge>& edges) {
    Matrix d;
    for (int s = 0; s < n; s++) d.push_back(dijkstra(n, edges, s));
    return d;
}

inline Matrix all_pairs_bellman_ford(int n, const std::vector<Edge>& edges) {
    Matrix d;
    for (int s = 0; s < n; s++) d.push_back(bellman_ford(n, edges, s).dist);
    return d;
}

// Bellman-Ford: el grafo donde Dijkstra falla, la tabla d_i(v), la versión de una sola fila y cuánto
// importa el orden de las aristas.
#include <iomanip>
#include <iostream>
#include <string>
#include <vector>

#include "shortest_paths.h"

namespace {

std::string show(long long x) { return x >= INF ? "inf" : std::to_string(x); }

void print_row(const std::vector<long long>& row) {
    for (long long x : row) std::cout << std::setw(5) << show(x);
    std::cout << '\n';
}

// Una cadena 0 -> 1 -> ... -> n - 1, con sus aristas en el orden dado.
std::vector<Edge> chain(int n, bool forward) {
    std::vector<Edge> edges;
    for (int i = 0; i + 1 < n; i++) {
        const int a = forward ? i : n - 2 - i;
        edges.push_back({a, a + 1, 1});
    }
    return edges;
}

}  // namespace

int main() {
    const std::vector<Edge> edges = {{0, 1, 1}, {0, 2, 4}, {2, 1, -5}, {1, 3, 2}};
    std::cout << "Aristas: 0->1 (1), 0->2 (4), 2->1 (-5), 1->3 (2)\n";
    std::cout << "Dijkstra:    ";
    print_row(dijkstra(4, edges, 0));

    std::cout << "\nLa tabla d_i(v):\n  i    0    1    2    3\n";
    const auto table = bellman_ford_table(4, edges, 0, 4);
    for (std::size_t i = 0; i < table.size(); i++) {
        std::cout << std::setw(3) << i;
        for (long long x : table[i]) std::cout << std::setw(5) << show(x);
        std::cout << '\n';
    }

    std::cout << "\nUna sola fila, con las aristas en el orden de arriba:\n";
    std::cout << "vuelta  0    1    2    3\n";
    bellman_ford(4, edges, 0, [](int round, const BellmanFordResult& r) {
        std::cout << std::setw(3) << round << ' ';
        print_row(r.dist);
    });

    const std::vector<Edge> shuffled = {{1, 3, 2}, {2, 1, -5}, {0, 2, 4}, {0, 1, 1}};
    std::cout << "\nUna sola fila, con las aristas al revés:\n";
    std::cout << "vuelta  0    1    2    3\n";
    bellman_ford(4, shuffled, 0, [](int round, const BellmanFordResult& r) {
        std::cout << std::setw(3) << round << ' ';
        print_row(r.dist);
    });

    std::cout << "\nCuántas vueltas hacen falta, hasta una que no cambia nada, en una cadena de n nodos:\n";
    std::cout << "    n  aristas en orden  aristas al revés\n";
    for (int n : {10, 100, 1000}) {
        std::cout << std::setw(5) << n << std::setw(18) << bellman_ford(n, chain(n, true), 0).rounds << std::setw(18)
                  << bellman_ford(n, chain(n, false), 0).rounds << '\n';
    }
}

// Detectar un ciclo negativo y recuperarlo: el nodo relajado en la vuelta V puede no estar en el ciclo.
#include <iostream>
#include <vector>

#include "shortest_paths.h"

namespace {

void print_values(const char* name, const std::vector<long long>& values) {
    std::cout << name;
    for (long long x : values) std::cout << ' ' << (x >= INF ? "inf" : std::to_string(x));
    std::cout << '\n';
}

}  // namespace

int main() {
    // 0 -> 1, el ciclo 1 -> 2 -> 3 -> 1 de peso -1, y la cola 3 -> 4 -> 5.
    const int n = 6;
    const std::vector<Edge> edges = {{0, 1, 2}, {1, 2, 1}, {2, 3, 1}, {3, 1, -3}, {3, 4, 1}, {4, 5, 1}};
    std::cout << "Aristas: 0->1 (2), 1->2 (1), 2->3 (1), 3->1 (-3), 3->4 (1), 4->5 (1)\n";
    std::cout << "vuelta: d[0..5]\n";
    const BellmanFordResult r = bellman_ford(n, edges, 0, [](int round, const BellmanFordResult& partial) {
        std::cout << round << ":";
        print_values("", partial.dist);
    });
    std::cout << "nodo relajado al final de la vuelta " << n << ": " << r.witness << '\n';
    std::cout << "pred:";
    for (int p : r.pred) std::cout << ' ' << p;
    std::cout << '\n';

    std::cout << "\nLeer directamente desde ese nodo, hasta volver a él (se corta a los 12 pasos):\n ";
    int x = r.witness;
    for (int step = 0; step < 12; step++) {
        std::cout << ' ' << x;
        x = r.pred[idx(x)];
        if (x == r.witness) break;
    }
    std::cout << " ...\n";

    std::cout << "\nRetroceder " << n << " veces y recién entonces leer:\n ";
    x = r.witness;
    for (int step = 0; step < n; step++) {
        std::cout << ' ' << x;
        x = r.pred[idx(x)];
    }
    std::cout << " -> " << x << '\n';
    const std::vector<int> cycle = recover_cycle(n, r.pred, r.witness);
    std::cout << "ciclo:";
    for (int v : cycle) std::cout << ' ' << v;
    std::cout << " -> " << cycle[0] << ", peso " << cycle_weight(edges, cycle) << '\n';
}

// Cada fila del cuadro de decisión con el mismo grafo de prueba: el algoritmo rápido frente a
// Bellman-Ford, que sirve en todas pero cuesta más.
#include <iostream>
#include <vector>

#include "all_pairs.h"

namespace {

void print_row(const char* name, const std::vector<long long>& row) {
    std::cout << "  " << name;
    for (long long x : row) std::cout << ' ' << (x >= INF ? "inf" : std::to_string(x));
    std::cout << '\n';
}

}  // namespace

int main() {
    std::cout << "Sin pesos: 0->1, 0->2, 1->3, 2->3, 3->4\n";
    std::vector<Edge> unweighted = {{0, 1, 1}, {0, 2, 1}, {1, 3, 1}, {2, 3, 1}, {3, 4, 1}};
    print_row("BFS:         ", bfs_distances(5, unweighted, 0));
    print_row("Bellman-Ford:", bellman_ford(5, unweighted, 0).dist);

    std::cout << "\nDAG con un peso negativo: 0->1 (1), 0->2 (4), 2->1 (-5), 1->3 (2)\n";
    std::vector<Edge> dag = {{0, 1, 1}, {0, 2, 4}, {2, 1, -5}, {1, 3, 2}};
    print_row("orden topológico:", dag_distances(4, dag, 0));
    print_row("Bellman-Ford:    ", bellman_ford(4, dag, 0).dist);
    print_row("Dijkstra:        ", dijkstra(4, dag, 0));

    std::cout << "\nPesos no negativos: 0->1 (7), 0->2 (1), 2->1 (2), 1->3 (1), 2->3 (6)\n";
    std::vector<Edge> roads = {{0, 1, 7}, {0, 2, 1}, {2, 1, 2}, {1, 3, 1}, {2, 3, 6}};
    print_row("Dijkstra:    ", dijkstra(4, roads, 0));
    print_row("Bellman-Ford:", bellman_ford(4, roads, 0).dist);

    std::cout << "\nCon un ciclo: el orden topológico no se aplica. Aristas 0->1 (1), 1->2 (1), 2->1 (1)\n";
    const std::vector<Edge> cyclic = {{0, 1, 1}, {1, 2, 1}, {2, 1, 1}};
    std::cout << "  dag_distances devuelve " << dag_distances(3, cyclic, 0).size() << " valores\n";

    std::cout << "\nTodos los pares con un negativo: 0->1 (3), 1->2 (-1), 2->0 (2), 2->1 (4)\n";
    const std::vector<Edge> small = {{0, 1, 3}, {1, 2, -1}, {2, 0, 2}, {2, 1, 4}};
    const Matrix floyd = floyd_warshall(3, small).dist;
    const Matrix repeated = all_pairs_bellman_ford(3, small);
    for (int s = 0; s < 3; s++) {
        print_row("Floyd-Warshall:", floyd[idx(s)]);
        print_row("V Bellman-Ford:", repeated[idx(s)]);
    }
}

// Floyd-Warshall: la matriz vuelta por vuelta, el ciclo de k por dentro, la diagonal y el crecimiento
// de los valores cuando hay un ciclo negativo y no se detiene.
#include <cstdio>
#include <limits>
#include <random>
#include <vector>

#include "all_pairs.h"

namespace {

void print_matrix(const Matrix& d) {
    for (const auto& row : d) {
        for (long long x : row) {
            if (x >= INF) {
                std::printf("%5s", "inf");
            } else {
                std::printf("%5lld", x);
            }
        }
        std::printf("\n");
    }
}

// Sin guardia, en double para ver el tamaño sin que haya overflow de enteros.
void growth_without_guard(int n) {
    std::vector<std::vector<double>> d(idx(n), std::vector<double>(idx(n), -1));
    for (int i = 0; i < n; i++) d[idx(i)][idx(i)] = 0;
    const double limit = static_cast<double>(std::numeric_limits<long long>::max());
    bool passed = false;
    for (int k = 0; k < n; k++) {
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++) {
                const double through = d[idx(i)][idx(k)] + d[idx(k)][idx(j)];
                if (through < d[idx(i)][idx(j)]) d[idx(i)][idx(j)] = through;
            }
        }
        double lowest = 0;
        for (const auto& row : d) {
            for (double x : row) {
                if (x < lowest) lowest = x;
            }
        }
        if (k < 6 || (!passed && -lowest > limit)) {
            std::printf("  tras la vuelta k = %2d: menor valor %.3e%s\n", k, lowest,
                        -lowest > limit ? "  (ya no cabe en long long)" : "");
        }
        if (-lowest > limit) passed = true;
    }
}

}  // namespace

int main() {
    const int n = 4;
    const std::vector<Edge> edges = {{0, 1, 4}, {0, 3, 9}, {1, 2, -2}, {2, 3, 1}, {3, 1, 5}, {2, 0, 3}};
    std::printf("Aristas: 0->1 (4), 0->3 (9), 1->2 (-2), 2->3 (1), 3->1 (5), 2->0 (3)\n");
    Matrix d = initial_matrix(n, edges);
    std::printf("d_0, sin nodos intermedios:\n");
    print_matrix(d);
    for (int k = 0; k < n; k++) {
        floyd_round(d, n, k);
        std::printf("d_%d, tras permitir el nodo %d:\n", k + 1, k);
        print_matrix(d);
    }

    const std::vector<Edge> trap = {{1, 3, 5}, {3, 2, 5}, {2, 0, 9}, {1, 0, 30}};
    std::printf("\nEl ciclo de k por dentro. Aristas: 1->3 (5), 3->2 (5), 2->0 (9), 1->0 (30)\n");
    const Matrix right = floyd_warshall(n, trap).dist;
    const Matrix wrong = floyd_warshall_k_inside(n, trap);
    std::printf("k por fuera:\n");
    print_matrix(right);
    std::printf("k por dentro:\n");
    print_matrix(wrong);

    std::printf("\nGrafos al azar en los que k por dentro se equivoca (20000 por fila, 35 %% de aristas, pesos 1 a 9):\n");
    std::mt19937 gen(10);
    for (int size : {3, 4, 5, 8}) {
        int bad = 0;
        for (int t = 0; t < 20000; t++) {
            std::vector<Edge> random_edges;
            for (int a = 0; a < size; a++) {
                for (int b = 0; b < size; b++) {
                    if (a != b && gen() % 100 < 35) random_edges.push_back({a, b, static_cast<long long>(gen() % 9) + 1});
                }
            }
            if (floyd_warshall(size, random_edges).dist != floyd_warshall_k_inside(size, random_edges)) bad++;
        }
        std::printf("  n = %d: %d\n", size, bad);
    }

    std::printf("\nLa diagonal. Aristas: 0->1 (1), 1->2 (-3), 2->1 (2), 2->3 (1)\n");
    const std::vector<Edge> negative = {{0, 1, 1}, {1, 2, -3}, {2, 1, 2}, {2, 3, 1}};
    Matrix m = initial_matrix(n, negative);
    for (int k = 0; k < n; k++) {
        std::printf("  antes de k = %d, diagonal:", k);
        for (int v = 0; v < n; v++) std::printf(" %lld", m[idx(v)][idx(v)]);
        std::printf("\n");
        if (m[idx(k)][idx(k)] < 0) {
            std::printf("  d[%d][%d] < 0: hay un ciclo negativo, y se detiene\n", k, k);
            break;
        }
        floyd_round(m, n, k);
    }

    std::printf("\nSin detenerse: 40 nodos, todas las aristas de peso -1.\n");
    growth_without_guard(40);
}

// Mide las permutaciones frente a Held-Karp. Con -O2, una corrida, en la máquina del momento.
#include <chrono>
#include <cstdio>
#include <random>
#include <vector>

#include "held_karp.h"

namespace {

template <typename Work>
double seconds_of(Work work) {
    const auto start = std::chrono::steady_clock::now();
    work();
    return std::chrono::duration<double>(std::chrono::steady_clock::now() - start).count();
}

}  // namespace

int main() {
    std::mt19937 gen(20);
    std::printf("%3s %14s %14s %16s\n", "n", "permutaciones", "Held-Karp", "memoria de tabla");
    for (int n = 8; n <= 20; n += (n < 12 ? 1 : 2)) {
        tsp::Costs cost(tsp::idx(n), std::vector<long long>(tsp::idx(n), 0));
        for (auto& row : cost) {
            for (long long& c : row) c = static_cast<long long>(gen() % 100) + 1;
        }
        long long a = 0;
        long long b = 0;
        const double t_table = seconds_of([&] { b = tsp::held_karp(cost); });
        const double megabytes = static_cast<double>((1ULL << (n - 1)) * static_cast<unsigned long long>(n) * 8) / 1e6;
        if (n <= 12) {
            const double t_brute = seconds_of([&] { a = tsp::brute_force(cost); });
            std::printf("%3d %12.3f s %12.3f s %13.1f MB  %s\n", n, t_brute, t_table, megabytes, a == b ? "iguales" : "DISTINTOS");
        } else {
            std::printf("%3d %14s %12.3f s %13.1f MB\n", n, "-", t_table, megabytes);
        }
    }
}

// Held-Karp en cuatro ciudades, casilla por casilla, y cuánto trabajo ahorra frente a las permutaciones.
#include <iomanip>
#include <iostream>
#include <random>
#include <string>
#include <vector>

#include "held_karp.h"

namespace {

std::string subset_name(std::size_t mask, int n) {
    std::string text = "{";
    for (int v = 1; v < n; v++) {
        if ((mask & (std::size_t{1} << (v - 1))) != 0) {
            if (text.size() > 1) text += ",";
            text += std::to_string(v);
        }
    }
    return text + "}";
}

tsp::Costs random_costs(int n, std::mt19937& gen) {
    tsp::Costs cost(tsp::idx(n), std::vector<long long>(tsp::idx(n), 0));
    for (auto& row : cost) {
        for (long long& c : row) c = static_cast<long long>(gen() % 100) + 1;
    }
    return cost;
}

}  // namespace

int main() {
    const tsp::Costs cost = {{0, 10, 15, 20}, {5, 0, 9, 10}, {6, 13, 0, 12}, {8, 8, 9, 0}};
    std::cout << "Costos:\n";
    for (const auto& row : cost) {
        for (long long c : row) std::cout << std::setw(4) << c;
        std::cout << '\n';
    }
    const auto table = tsp::held_karp_table(cost);
    std::cout << "\nmáscara  subconjunto  termina en 1  en 2  en 3\n";
    for (std::size_t mask = 1; mask < table.size(); mask++) {
        std::cout << std::setw(7) << mask << "  " << std::setw(11) << subset_name(mask, 4);
        for (int v = 1; v <= 3; v++) {
            const long long value = table[mask][tsp::idx(v)];
            std::cout << std::setw(v == 1 ? 14 : 6) << (value == tsp::NONE ? std::string("-") : std::to_string(value));
        }
        std::cout << '\n';
    }
    std::cout << "\nmejor recorrido, Held-Karp: " << tsp::held_karp(cost) << '\n';
    std::cout << "mejor recorrido, permutaciones: " << tsp::brute_force(cost) << '\n';

    std::cout << "\n n  recorridos (n-1)!  transiciones  n^2 2^n     mismo costo\n";
    std::mt19937 gen(15);
    for (int n = 4; n <= 11; n++) {
        const tsp::Costs random = random_costs(n, gen);
        long long tours = 0;
        long long transitions = 0;
        const long long by_brute = tsp::brute_force(random, &tours);
        const long long by_table = tsp::held_karp(random, &transitions);
        std::cout << std::setw(2) << n << std::setw(19) << tours << std::setw(14) << transitions << std::setw(9)
                  << n * n * (1LL << n) << std::setw(16) << (by_brute == by_table ? "sí" : "NO") << '\n';
    }
}

// El problema del viajero: fuerza bruta sobre las permutaciones, y Held-Karp sobre subconjuntos.
// Un subconjunto de nodos es un entero: el bit v vale 1 si v está en el subconjunto.
#pragma once

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

namespace tsp {

using Costs = std::vector<std::vector<long long>>;

// Un costo que ningún recorrido alcanza. Marca que no hay arista.
inline constexpr long long NONE = std::numeric_limits<long long>::max() / 4;

inline std::size_t idx(int v) { return static_cast<std::size_t>(v); }

// El costo mínimo de un recorrido que sale de 0, pasa una vez por cada nodo y vuelve a 0.
// Prueba las (n - 1)! órdenes. Devuelve NONE si no hay recorrido.
inline long long brute_force(const Costs& cost, long long* tours = nullptr) {
    const int n = static_cast<int>(cost.size());
    if (n <= 1) return 0;
    std::vector<int> middle(idx(n - 1));
    std::iota(middle.begin(), middle.end(), 1);
    long long best = NONE;
    long long count = 0;
    do {
        count++;
        long long total = cost[0][idx(middle[0])];
        for (std::size_t i = 0; i + 1 < middle.size() && total < NONE; i++) {
            const long long step = cost[idx(middle[i])][idx(middle[i + 1])];
            total = step == NONE ? NONE : total + step;
        }
        const long long back = cost[idx(middle.back())][0];
        if (total < NONE && back < NONE) best = std::min(best, total + back);
    } while (std::next_permutation(middle.begin(), middle.end()));
    if (tours != nullptr) *tours = count;
    return best;
}

// Held-Karp. table[mask][v] es el costo mínimo de un camino que sale de 0, visita exactamente los
// nodos de mask y termina en v, con v en mask. Los bits 1, ..., n - 1 son los nodos distintos de
// 0; el nodo 0 no se guarda en mask. Por eso el bit v - 1 representa al nodo v.
inline std::vector<std::vector<long long>> held_karp_table(const Costs& cost, long long* transitions = nullptr) {
    const int n = static_cast<int>(cost.size());
    const std::size_t masks = std::size_t{1} << (n - 1);
    std::vector<std::vector<long long>> table(masks, std::vector<long long>(idx(n), NONE));
    long long count = 0;
    for (int v = 1; v < n; v++) table[std::size_t{1} << (v - 1)][idx(v)] = cost[0][idx(v)];
    // Quitar un bit da un entero menor: al llegar a mask, sus subconjuntos ya están listos.
    for (std::size_t mask = 1; mask < masks; mask++) {
        for (int v = 1; v < n; v++) {
            const std::size_t bit_v = std::size_t{1} << (v - 1);
            if ((mask & bit_v) == 0 || mask == bit_v) continue;
            const std::size_t without_v = mask ^ bit_v;
            for (int u = 1; u < n; u++) {
                if ((without_v & (std::size_t{1} << (u - 1))) == 0) continue;
                count++;
                if (table[without_v][idx(u)] == NONE || cost[idx(u)][idx(v)] == NONE) continue;
                table[mask][idx(v)] = std::min(table[mask][idx(v)], table[without_v][idx(u)] + cost[idx(u)][idx(v)]);
            }
        }
    }
    if (transitions != nullptr) *transitions = count;
    return table;
}

// El costo del mejor recorrido: se elige el último nodo antes de volver a 0.
inline long long held_karp(const Costs& cost, long long* transitions = nullptr) {
    const int n = static_cast<int>(cost.size());
    if (n <= 1) return 0;
    const auto table = held_karp_table(cost, transitions);
    const std::size_t masks = table.size();
    long long best = NONE;
    for (int v = 1; v < n; v++) {
        if (table[masks - 1][idx(v)] != NONE && cost[idx(v)][0] != NONE) {
            best = std::min(best, table[masks - 1][idx(v)] + cost[idx(v)][0]);
        }
    }
    return best;
}

}  // namespace tsp

// Los dos errores de partir de un "infinito" y sumarle algo a un nodo no alcanzado.
// El primero es un overflow de int: es comportamiento indefinido, y aquí lo señala el sanitizer.
#include <climits>
#include <iostream>
#include <vector>

#include "shortest_paths.h"

namespace {

void print_row(const char* name, const std::vector<long long>& row) {
    std::cout << name;
    for (long long x : row) std::cout << ' ' << (x >= INF ? "inf" : std::to_string(x));
    std::cout << '\n';
}

// Primer error: INT_MAX como infinito, sin mirar si el nodo fue alcanzado.
std::vector<int> relax_with_int_max(int n, const std::vector<Edge>& edges, int source) {
    std::vector<int> d(idx(n), INT_MAX);
    d[idx(source)] = 0;
    for (int round = 1; round < n; round++) {
        for (const Edge& e : edges) {
            if (d[idx(e.from)] + static_cast<int>(e.weight) < d[idx(e.to)]) {  // overflow si d[e.from] es INT_MAX
                d[idx(e.to)] = d[idx(e.from)] + static_cast<int>(e.weight);
            }
        }
    }
    return d;
}

// Segundo error: un infinito sin overflow, pero tampoco se trata como especial.
std::vector<long long> relax_without_check(int n, const std::vector<Edge>& edges, int source, bool& changed_last) {
    const long long big = 1'000'000'000'000;
    std::vector<long long> d(idx(n), big);
    d[idx(source)] = 0;
    changed_last = false;
    for (int round = 1; round <= n; round++) {
        changed_last = false;
        for (const Edge& e : edges) {
            if (d[idx(e.from)] + e.weight < d[idx(e.to)]) {
                d[idx(e.to)] = d[idx(e.from)] + e.weight;
                changed_last = true;
            }
        }
    }
    return d;
}

}  // namespace

int main() {
    // El origen es 0, que solo llega a 1. Los nodos 2 y 3 no se alcanzan.
    const std::vector<Edge> positive = {{2, 3, 5}, {0, 1, 7}};
    std::cout << "Primer error: INT_MAX como infinito. Aristas 2->3 (5), 0->1 (7).\n";
    std::vector<int> wrong = relax_with_int_max(4, positive, 0);
    std::cout << "d con INT_MAX:";
    for (int x : wrong) std::cout << ' ' << x;
    std::cout << '\n';
    print_row("correcto:     ", bellman_ford(4, positive, 0).dist);

    std::cout << "\nSegundo error: infinito + un negativo. Aristas 0->1 (7), 2->3 (-5).\n";
    const std::vector<Edge> negative = {{0, 1, 7}, {2, 3, -5}};
    bool changed = false;
    print_row("sin revisar:  ", relax_without_check(4, negative, 0, changed));
    print_row("revisando:    ", bellman_ford(4, negative, 0).dist);

    std::cout << "\nY un falso ciclo: 0->1 (7), y el ciclo 2->3 (-5), 3->2 (1), que el origen no alcanza.\n";
    const std::vector<Edge> unreachable = {{0, 1, 7}, {2, 3, -5}, {3, 2, 1}};
    relax_without_check(4, unreachable, 0, changed);
    std::cout << "sin revisar: la vuelta 4 cambió algo: " << (changed ? "sí" : "no") << '\n';
    const BellmanFordResult r = bellman_ford(4, unreachable, 0);
    std::cout << "revisando:   la vuelta 4 cambió algo: " << (r.witness != -1 ? "sí" : "no") << '\n';
}

// De productos a sumas, y lo que el redondeo le hace a esa cuenta.
#include <cmath>
#include <cstdio>
#include <vector>

#include "rates.h"

int main() {
    std::printf("Producto real y suma de -ln(tasa) en double:\n");
    std::printf("%-22s %10s %12s\n", "tasas", "producto", "suma");
    const std::vector<std::pair<const char*, std::vector<double>>> chains = {
        {"2 y 0,5", {2, 0.5}},
        {"1,25 y 0,8", {1.25, 0.8}},
        {"1,6 y 0,625", {1.6, 0.625}},
        {"2,5 y 0,4", {2.5, 0.4}},
        {"1,6, 1,25 y 0,5", {1.6, 1.25, 0.5}},
        {"1,25, 1,25 y 0,64", {1.25, 1.25, 0.64}},
    };
    for (const auto& [name, rates] : chains) std::printf("%-22s %10s %12.3e\n", name, "1 exacto", log_sum(rates));
    std::printf("%-22s %10s %12.3e\n", "1,1, 1,1 y 1/1,21 (*)", "menor a 1", log_sum({1.1, 1.1, 0.826446280991735}));
    std::printf("(*) 0,826446280991735 es 1/1,21 truncado a 15 decimales.\n");

    std::printf("\nCon pesos -ln(tasa) + epsilon, ¿la vuelta n relaja algo?\n");
    std::printf("%10s %28s %28s\n", "epsilon", "tabla exacta (sin arbitraje)", "ciclo de 400, 1,000001");
    for (double epsilon : {0.0, 1e-15, 1e-12, 1e-9, 1e-8, 1e-6}) {
        const bool table = round_n_relaxes(5, consistent_table(), epsilon);
        const bool cycle = round_n_relaxes(400, long_cycle(400), epsilon);
        std::printf("%10.0e %28s %28s\n", epsilon, table ? "sí (falso arbitraje)" : "no", cycle ? "sí" : "no (se le escapa)");
    }
}

// Tipos de cambio: un producto mayor que 1 es una suma negativa de -ln(tasa). Estas funciones solo
// responden si, desde la moneda 0, se detecta un ciclo. Con ellas se mide qué hace el redondeo.
#pragma once

#include <cmath>
#include <cstddef>
#include <limits>
#include <vector>

struct Rate {
    int from;
    int to;
    double rate;  // una unidad de from compra rate unidades de to
};

// La suma de -ln(tasa) a lo largo de una cadena de tasas.
inline double log_sum(const std::vector<double>& rates) {
    double total = 0;
    for (double r : rates) total += -std::log(r);
    return total;
}

// Bellman-Ford desde la moneda 0, con pesos -ln(tasa) + epsilon. Dice si la vuelta n relajó algo.
inline bool round_n_relaxes(int n, const std::vector<Rate>& rates, double epsilon) {
    const double unreached = std::numeric_limits<double>::infinity();
    std::vector<double> d(static_cast<std::size_t>(n), unreached);
    d[0] = 0;
    bool relaxed = false;
    for (int round = 1; round <= n; round++) {
        relaxed = false;
        for (const Rate& e : rates) {
            const auto u = static_cast<std::size_t>(e.from);
            const auto v = static_cast<std::size_t>(e.to);
            if (d[u] == unreached) continue;
            const double through = d[u] + (-std::log(e.rate) + epsilon);
            if (through < d[v]) {
                d[v] = through;
                relaxed = true;
            }
        }
        if (!relaxed) return false;
    }
    return relaxed;
}

// Cinco monedas con valores 1, 1,25, 1,6, 2,5 y 0,4, y la tasa exacta entre cada par. Todos los
// productos alrededor de un ciclo valen exactamente 1, así que no hay arbitraje.
inline std::vector<Rate> consistent_table() {
    return {{0, 1, 1.25},  {0, 2, 1.6},  {0, 3, 2.5},    {0, 4, 0.4},  {1, 0, 0.8},  {1, 2, 1.28}, {1, 3, 2},
            {1, 4, 0.32},  {2, 0, 0.625}, {2, 1, 0.78125}, {2, 3, 1.5625}, {2, 4, 0.25}, {3, 0, 0.4},  {3, 1, 0.5},
            {3, 2, 0.64},  {3, 4, 0.16}, {4, 0, 2.5},    {4, 1, 3.125}, {4, 2, 4},    {4, 3, 6.25}};
}

// Un ciclo de largo L, todas las tasas 1 menos una: su producto es exactamente 1,000001.
inline std::vector<Rate> long_cycle(int length) {
    std::vector<Rate> rates;
    for (int i = 0; i < length; i++) rates.push_back({i, (i + 1) % length, i == 0 ? 1.000001 : 1.0});
    return rates;
}

// Caminos mínimos desde un origen: Bellman-Ford en sus dos formas, la recuperación de un ciclo
// negativo, y los algoritmos con los que se compara (BFS, DAG y Dijkstra).
#pragma once

#include <cstddef>
#include <functional>
#include <limits>
#include <queue>
#include <utility>
#include <vector>

struct Edge {
    int from;
    int to;
    long long weight;
};

// Un valor que ningún camino alcanza, y lejos del límite de long long.
inline constexpr long long INF = std::numeric_limits<long long>::max() / 4;

inline std::size_t idx(int v) { return static_cast<std::size_t>(v); }

// La tabla completa: table[i][v] es d_i(v), el menor largo de un camino de source a v con a lo
// sumo i aristas. Tiene las filas 0, ..., rows.
inline std::vector<std::vector<long long>> bellman_ford_table(int n, const std::vector<Edge>& edges, int source,
                                                              int rows) {
    std::vector<std::vector<long long>> table(idx(rows) + 1, std::vector<long long>(idx(n), INF));
    table[0][idx(source)] = 0;
    for (int i = 1; i <= rows; i++) {
        table[idx(i)] = table[idx(i - 1)];  // el caso de a lo sumo i - 1 aristas
        for (const Edge& e : edges) {       // el caso de exactamente i aristas, con e al final
            if (table[idx(i - 1)][idx(e.from)] == INF) continue;
            const long long through = table[idx(i - 1)][idx(e.from)] + e.weight;
            if (through < table[idx(i)][idx(e.to)]) table[idx(i)][idx(e.to)] = through;
        }
    }
    return table;
}

struct BellmanFordResult {
    std::vector<long long> dist;
    std::vector<int> pred;  // pred[v] es el nodo desde el que se relajó v por última vez, o -1
    int rounds = 0;         // vueltas que se hicieron
    int witness = -1;       // un nodo relajado en la vuelta n, o -1 si no hubo ninguno
};

// Una sola fila, hasta n vueltas, y se detiene si una vuelta no cambia nada.
// Invariante: tras la vuelta i, d[v] es el largo de algún camino de source a v (o INF), y d[v] <= d_i(v).
// after_round, si se pasa, se llama al final de cada vuelta.
inline BellmanFordResult bellman_ford(int n, const std::vector<Edge>& edges, int source,
                                      const std::function<void(int, const BellmanFordResult&)>& after_round = {}) {
    BellmanFordResult result;
    result.dist.assign(idx(n), INF);
    result.pred.assign(idx(n), -1);
    result.dist[idx(source)] = 0;
    for (int round = 1; round <= n; round++) {
        result.rounds = round;
        int last_relaxed = -1;
        for (const Edge& e : edges) {
            if (result.dist[idx(e.from)] == INF) continue;  // no alcanzado: no hay nada que propagar
            if (result.dist[idx(e.from)] + e.weight < result.dist[idx(e.to)]) {
                result.dist[idx(e.to)] = result.dist[idx(e.from)] + e.weight;
                result.pred[idx(e.to)] = e.from;
                last_relaxed = e.to;
            }
        }
        if (after_round) after_round(round, result);
        if (last_relaxed == -1) break;
        if (round == n) result.witness = last_relaxed;
    }
    return result;
}

// Dado un nodo relajado en la vuelta n, devuelve un ciclo negativo: c[0], c[1], ..., con
// c[i] -> c[i + 1] aristas del grafo y c.back() -> c[0] la que lo cierra.
inline std::vector<int> recover_cycle(int n, const std::vector<int>& pred, int witness) {
    int x = witness;
    for (int step = 0; step < n; step++) x = pred[idx(x)];  // ahora x está en el ciclo
    std::vector<int> cycle;
    int y = x;
    do {
        cycle.push_back(y);
        y = pred[idx(y)];
    } while (y != x);
    // Se leyó siguiendo a los predecesores, así que está al revés.
    std::vector<int> forward(cycle.rbegin(), cycle.rend());
    return forward;
}

// El menor peso de una arista de a hacia b. Sirve para sumar un ciclo o un camino ya recuperado.
inline long long edge_weight(const std::vector<Edge>& edges, int a, int b) {
    long long best = INF;
    for (const Edge& e : edges) {
        if (e.from == a && e.to == b && e.weight < best) best = e.weight;
    }
    return best;
}

inline long long cycle_weight(const std::vector<Edge>& edges, const std::vector<int>& cycle) {
    long long total = 0;
    for (std::size_t i = 0; i < cycle.size(); i++) {
        total += edge_weight(edges, cycle[i], cycle[(i + 1) % cycle.size()]);
    }
    return total;
}

inline std::vector<std::vector<std::pair<int, long long>>> adjacency(int n, const std::vector<Edge>& edges) {
    std::vector<std::vector<std::pair<int, long long>>> adj(idx(n));
    for (const Edge& e : edges) adj[idx(e.from)].push_back({e.to, e.weight});
    return adj;
}

// BFS: sirve si todas las aristas pesan lo mismo. Cuenta aristas.
inline std::vector<long long> bfs_distances(int n, const std::vector<Edge>& edges, int source) {
    const auto adj = adjacency(n, edges);
    std::vector<long long> dist(idx(n), INF);
    std::queue<int> pending;
    dist[idx(source)] = 0;
    pending.push(source);
    while (!pending.empty()) {
        const int u = pending.front();
        pending.pop();
        for (const auto& [v, weight] : adj[idx(u)]) {
            (void)weight;
            if (dist[idx(v)] == INF) {
                dist[idx(v)] = dist[idx(u)] + 1;
                pending.push(v);
            }
        }
    }
    return dist;
}

// DAG: procesa los nodos en orden topológico, con pesos de cualquier signo. Devuelve una lista
// vacía si el grafo tiene un ciclo.
inline std::vector<long long> dag_distances(int n, const std::vector<Edge>& edges, int source) {
    const auto adj = adjacency(n, edges);
    std::vector<int> in_degree(idx(n), 0);
    for (const Edge& e : edges) in_degree[idx(e.to)]++;
    std::vector<int> order;
    for (int v = 0; v < n; v++) {
        if (in_degree[idx(v)] == 0) order.push_back(v);
    }
    std::vector<long long> dist(idx(n), INF);
    dist[idx(source)] = 0;
    for (std::size_t next = 0; next < order.size(); next++) {
        const int u = order[next];
        for (const auto& [v, weight] : adj[idx(u)]) {
            if (dist[idx(u)] != INF && dist[idx(u)] + weight < dist[idx(v)]) dist[idx(v)] = dist[idx(u)] + weight;
            if (--in_degree[idx(v)] == 0) order.push_back(v);
        }
    }
    if (order.size() != idx(n)) return {};
    return dist;
}

// Dijkstra con heap y borrado perezoso, de la semana 5. Solo es correcto con pesos no negativos: si
// hay negativos, un nodo fijado puede bajar después, pero nadie vuelve a leer sus aristas.
inline std::vector<long long> dijkstra(int n, const std::vector<Edge>& edges, int source) {
    const auto adj = adjacency(n, edges);
    std::vector<long long> dist(idx(n), INF);
    using Entry = std::pair<long long, int>;
    std::priority_queue<Entry, std::vector<Entry>, std::greater<Entry>> heap;
    std::vector<bool> settled(idx(n), false);
    dist[idx(source)] = 0;
    heap.push({0, source});
    while (!heap.empty()) {
        const int u = heap.top().second;
        heap.pop();
        if (settled[idx(u)]) continue;
        settled[idx(u)] = true;
        for (const auto& [v, weight] : adj[idx(u)]) {
            if (dist[idx(u)] + weight < dist[idx(v)]) {
                dist[idx(v)] = dist[idx(u)] + weight;
                heap.push({dist[idx(v)], v});
            }
        }
    }
    return dist;
}