// Medidor de tiempos del curso. Compilar con -O2 y sin sanitizers.
#ifndef BENCH_H
#define BENCH_H

#include <algorithm>
#include <chrono>
#include <vector>

// Corre fn `reps` veces y devuelve la mediana en segundos.
template <typename Fn>
double median_seconds(int reps, Fn fn) {
    std::vector<double> times;
    times.reserve(reps);
    for (int i = 0; i < reps; i++) {
        auto start = std::chrono::steady_clock::now();
        fn();
        auto stop = std::chrono::steady_clock::now();
        times.push_back(std::chrono::duration<double>(stop - start).count());
    }
    std::sort(times.begin(), times.end());
    return times[times.size() / 2];
}

// Evita que el compilador borre un cálculo cuyo resultado no se usa.
template <typename T>
void keep(const T& value) {
    volatile T sink = value;
    (void)sink;
}

#endif

// Mide la multiplicación de la escuela contra Karatsuba, con -O2 y sin sanitizers.
// Los tiempos dependen de la máquina y de esta implementación. Cada dato es la mediana de varias corridas,
// y cada corrida repite la multiplicación hasta sumar al menos unos milisegundos.
//
//   build/karatsuba_bench            hasta n = 65536
//   build/karatsuba_bench 262144     hasta el n que se indique, potencia de 2; tarda varios minutos
#include <cstdio>
#include <cstdlib>
#include <random>
#include <vector>

#include "bench.h"
#include "karatsuba.h"

static Digits random_number(int n, std::mt19937& rng) {
    Digits a(n);
    for (int& digit : a) {
        digit = static_cast<int>(rng() % 10);
    }
    a[n - 1] = 1 + static_cast<int>(rng() % 9);
    return a;
}

// Microsegundos por multiplicación: 7 corridas, o 3 si cada una tarda más de un segundo.
template <typename Fn>
double micros_per_call(Fn fn) {
    int repeats = 1;
    double first = 0;
    while (true) {
        first = median_seconds(1, [&] {
            for (int i = 0; i < repeats; i++) {
                fn();
            }
        });
        if (first >= 0.005) {
            break;
        }
        repeats *= 2;
    }
    int runs = first > 1.0 ? 3 : 7;
    return 1e6 * median_seconds(runs, [&] {
               for (int i = 0; i < repeats; i++) {
                   fn();
               }
           }) / repeats;
}

int main(int argc, char** argv) {
    int largest = argc > 1 ? std::atoi(argv[1]) : 65536;
    std::mt19937 rng(2026);
    const std::size_t hybrid = 128;
    std::printf("microsegundos por multiplicacion\n");
    std::printf("%8s %12s %12s %12s %9s %9s\n", "n", "escuela", "umbral 3", "umbral 128", "esc/u3", "esc/u128");
    for (int n = 16; n <= largest; n *= 2) {
        Digits x = random_number(n, rng);
        Digits y = random_number(n, rng);
        double school = micros_per_call([&] {
            long long ignored = 0;
            keep(school_multiply(x, y, ignored).size());
        });
        double pure = micros_per_call([&] {
            long long ignored = 0;
            keep(karatsuba(x, y, 3, ignored).size());
        });
        double mixed = micros_per_call([&] {
            long long ignored = 0;
            keep(karatsuba(x, y, hybrid, ignored).size());
        });
        std::printf("%8d %12.1f %12.1f %12.1f %9.2f %9.2f\n", n, school, pure, mixed, school / pure, school / mixed);
    }

    const int n = 16384;
    Digits x = random_number(n, rng);
    Digits y = random_number(n, rng);
    std::printf("\nn = %d con distintos umbrales\n", n);
    std::printf("%8s %12s\n", "umbral", "microseg");
    for (std::size_t threshold : {3, 8, 16, 32, 64, 128, 256, 512, 1024}) {
        double kar = micros_per_call([&] {
            long long ignored = 0;
            keep(karatsuba(x, y, threshold, ignored).size());
        });
        std::printf("%8zu %12.1f\n", threshold, kar);
    }
}

// Cuenta productos de un dígito por un dígito. Es un conteo exacto: no depende de la máquina.
#include <cstdio>
#include <random>

#include "karatsuba.h"

static Digits random_number(int n, std::mt19937& rng) {
    Digits a(n);
    for (int& digit : a) {
        digit = static_cast<int>(rng() % 10);
    }
    a[n - 1] = 1 + static_cast<int>(rng() % 9);
    return a;
}

int main() {
    std::mt19937 rng(2026);
    std::printf("%8s %14s %8s %14s %8s\n", "n", "escuela", "razon", "karatsuba", "razon");
    long long previous_school = 0;
    long long previous_kar = 0;
    for (int n = 16; n <= 65536; n *= 2) {
        Digits x = random_number(n, rng);
        Digits y = random_number(n, rng);
        long long school = 0;
        long long kar = 0;
        Digits slow = n <= 8192 ? school_multiply(x, y, school) : Digits{};
        if (n > 8192) {
            school = static_cast<long long>(n) * n;
        }
        Digits fast = karatsuba(x, y, 3, kar);
        if (n <= 8192) {
            assert(slow == fast);
        }
        std::printf("%8d %14lld %8.2f %14lld %8.2f\n", n, school,
                    previous_school == 0 ? 0.0 : static_cast<double>(school) / static_cast<double>(previous_school),
                    kar,
                    previous_kar == 0 ? 0.0 : static_cast<double>(kar) / static_cast<double>(previous_kar));
        previous_school = school;
        previous_kar = kar;
    }
}

// Sigue el primer nivel de Karatsuba con 1234 x 5678.
#include <iostream>

#include "karatsuba.h"

int main() {
    long long products = 0;
    Digits x = from_string("1234");
    Digits y = from_string("5678");
    Digits a = from_string("12");
    Digits b = from_string("34");
    Digits c = from_string("56");
    Digits d = from_string("78");
    Digits ac = karatsuba(a, c, 3, products);
    Digits bd = karatsuba(b, d, 3, products);
    Digits sum_ab = add(a, b);
    Digits sum_cd = add(c, d);
    Digits sums = karatsuba(sum_ab, sum_cd, 3, products);
    Digits middle = subtract(subtract(sums, ac), bd);
    std::cout << "a = 12, b = 34, c = 56, d = 78, m = 2\n";
    std::cout << "ac = " << to_string(ac) << "\n";
    std::cout << "bd = " << to_string(bd) << "\n";
    std::cout << "(a + b)(c + d) = " << to_string(sum_ab) << " x " << to_string(sum_cd) << " = " << to_string(sums) << "\n";
    std::cout << "ad + bc = " << to_string(sums) << " - " << to_string(ac) << " - " << to_string(bd) << " = " << to_string(middle) << "\n";
    long long other = 0;
    std::cout << "resultado = " << to_string(karatsuba(x, y, 3, other)) << "\n";
    std::cout << "la escuela da = " << to_string(school_multiply(x, y, other)) << "\n";
}

// Multiplicación de enteros no negativos de muchos dígitos.
// Un número es un vector de dígitos en base 10, del menos significativo al más significativo,
// sin ceros a la izquierda. El cero es el vector vacío.
#ifndef KARATSUBA_H
#define KARATSUBA_H

#include <algorithm>
#include <cassert>
#include <cstddef>
#include <string>
#include <vector>

using Digits = std::vector<int>;

inline void trim(Digits& a) {
    while (!a.empty() && a.back() == 0) {
        a.pop_back();
    }
}

inline Digits from_string(const std::string& text) {
    Digits a;
    for (auto it = text.rbegin(); it != text.rend(); ++it) {
        a.push_back(*it - '0');
    }
    trim(a);
    return a;
}

inline std::string to_string(const Digits& a) {
    if (a.empty()) {
        return "0";
    }
    std::string text;
    for (auto it = a.rbegin(); it != a.rend(); ++it) {
        text.push_back(static_cast<char>('0' + *it));
    }
    return text;
}

inline Digits add(const Digits& a, const Digits& b) {
    Digits sum;
    sum.reserve(std::max(a.size(), b.size()) + 1);
    int carry = 0;
    for (std::size_t i = 0; i < a.size() || i < b.size() || carry != 0; i++) {
        int value = carry;
        if (i < a.size()) {
            value += a[i];
        }
        if (i < b.size()) {
            value += b[i];
        }
        sum.push_back(value % 10);
        carry = value / 10;
    }
    return sum;
}

// Resta a - b. Requiere a >= b.
inline Digits subtract(const Digits& a, const Digits& b) {
    Digits result = a;
    int borrow = 0;
    for (std::size_t i = 0; i < a.size(); i++) {
        int value = a[i] - borrow - (i < b.size() ? b[i] : 0);
        borrow = value < 0 ? 1 : 0;
        result[i] = value + 10 * borrow;
    }
    assert(borrow == 0);
    trim(result);
    return result;
}

// Multiplica por 10^positions.
inline Digits shifted(const Digits& a, std::size_t positions) {
    if (a.empty()) {
        return a;
    }
    Digits result(positions, 0);
    result.insert(result.end(), a.begin(), a.end());
    return result;
}

// Los dígitos desde `start` hasta antes de `stop`, sin ceros a la izquierda.
inline Digits slice(const Digits& a, std::size_t start, std::size_t stop) {
    start = std::min(start, a.size());
    stop = std::min(stop, a.size());
    Digits part(a.begin() + static_cast<long>(start), a.begin() + static_cast<long>(stop));
    trim(part);
    return part;
}

// La multiplicación de la escuela: cada dígito de a con cada dígito de b.
// `digit_products` cuenta los productos de un dígito por un dígito.
inline Digits school_multiply(const Digits& a, const Digits& b, long long& digit_products) {
    if (a.empty() || b.empty()) {
        return {};
    }
    std::vector<long long> cells(a.size() + b.size(), 0);
    for (std::size_t i = 0; i < a.size(); i++) {
        for (std::size_t j = 0; j < b.size(); j++) {
            cells[i + j] += a[i] * b[j];
        }
    }
    digit_products += static_cast<long long>(a.size()) * static_cast<long long>(b.size());
    Digits result(cells.size());
    long long carry = 0;
    for (std::size_t i = 0; i < cells.size(); i++) {
        long long value = cells[i] + carry;
        result[i] = static_cast<int>(value % 10);
        carry = value / 10;
    }
    trim(result);
    return result;
}

// Karatsuba: x = a·10^m + b, y = c·10^m + d, y tres productos en vez de cuatro.
// Con n <= threshold usa la multiplicación de la escuela. threshold debe ser al menos 3: a + b
// puede tener un dígito más que sus partes, y solo con n >= 4 los subproblemas son más cortos que n.
inline Digits karatsuba(const Digits& x, const Digits& y, std::size_t threshold, long long& digit_products) {
    assert(threshold >= 3);
    std::size_t n = std::max(x.size(), y.size());
    if (x.empty() || y.empty()) {
        return {};
    }
    if (n <= threshold) {
        return school_multiply(x, y, digit_products);
    }
    std::size_t m = n / 2;
    Digits b = slice(x, 0, m);
    Digits a = slice(x, m, x.size());
    Digits d = slice(y, 0, m);
    Digits c = slice(y, m, y.size());

    Digits ac = karatsuba(a, c, threshold, digit_products);
    Digits bd = karatsuba(b, d, threshold, digit_products);
    Digits sums = karatsuba(add(a, b), add(c, d), threshold, digit_products);
    Digits middle = subtract(subtract(sums, ac), bd);  // ad + bc

    Digits result = add(add(shifted(ac, 2 * m), shifted(middle, m)), bd);
    trim(result);
    return result;
}

#endif

// Sigue una vuelta de la mediana de medianas con 23 elementos: grupos, medianas, pivote y conteos.
#include <algorithm>
#include <cstddef>
#include <iostream>
#include <vector>

#include "select.h"

static void print(const std::vector<int>& v) {
    for (std::size_t i = 0; i < v.size(); i++) {
        std::cout << (i == 0 ? "" : " ") << v[i];
    }
}

// Sigue select_linear: una línea por llamada, sin entrar en la búsqueda de cada pivote.
static int trace_select(std::vector<int> a, std::size_t k) {
    long long ignored = 0;
    std::cout << "|A| = " << a.size() << ", k = " << k;
    if (a.size() <= 5) {
        std::sort(a.begin(), a.end());
        std::cout << ": se ordena [";
        print(a);
        std::cout << "] y la respuesta es " << a[k] << "\n";
        return a[k];
    }
    int pivot = median_of_medians(a, ignored);
    Parts parts = partition3(a, pivot, ignored);
    std::cout << ": pivote " << pivot << ", |L| = " << parts.less.size() << ", |E| = " << parts.equal.size()
              << ", |G| = " << parts.greater.size() << "\n";
    if (k < parts.less.size()) {
        return trace_select(parts.less, k);
    }
    if (k < parts.less.size() + parts.equal.size()) {
        std::cout << "k cae en E: la respuesta es " << pivot << "\n";
        return pivot;
    }
    return trace_select(parts.greater, k - parts.less.size() - parts.equal.size());
}

int main() {
    std::vector<int> a = {12, 3, 25, 8, 17, 40, 6, 31, 22, 9, 15, 37, 1, 28, 19, 44, 11, 33, 5, 26, 14, 30, 2};
    long long comparisons = 0;
    std::cout << "n = " << a.size() << "\n";
    for (std::size_t start = 0; start < a.size(); start += 5) {
        std::size_t stop = std::min(start + 5, a.size());
        std::vector<int> group(a.begin() + start, a.begin() + stop);
        std::cout << "grupo: ";
        print(group);
        sort_group(group, comparisons);
        std::cout << "  ordenado: ";
        print(group);
        std::cout << "  mediana: " << group[(group.size() - 1) / 2] << "\n";
    }
    comparisons = 0;
    std::vector<int> medians = group_medians(a, comparisons);
    std::cout << "medianas: ";
    print(medians);
    std::cout << "\n";
    int pivot = median_of_medians(a, comparisons);
    std::cout << "pivote (mediana inferior de las medianas): " << pivot << "\n";
    auto at_most = std::count_if(a.begin(), a.end(), [&](int x) { return x <= pivot; });
    auto at_least = std::count_if(a.begin(), a.end(), [&](int x) { return x >= pivot; });
    std::cout << "elementos <= pivote: " << at_most << "\n";
    std::cout << "elementos >= pivote: " << at_least << "\n";
    std::cout << "cota 3n/10 - 3: " << 3.0 * static_cast<double>(a.size()) / 10 - 3 << "\n";
    long long ignored = 0;
    Parts parts = partition3(a, pivot, ignored);
    std::cout << "menores: " << parts.less.size() << "  iguales: " << parts.equal.size()
              << "  mayores: " << parts.greater.size() << "\n";
    std::cout << "cota 7n/10 + 3: " << 7.0 * static_cast<double>(a.size()) / 10 + 3 << "\n";

    std::cout << "\nla mediana, k = 11\n";
    int answer = trace_select(a, 11);
    std::vector<int> sorted = a;
    std::sort(sorted.begin(), sorted.end());
    std::cout << "ordenando, la posicion 11 es " << sorted[11] << "; la funcion dio " << answer << "\n";
}

// Cuántos elementos deja a cada lado el pivote de la mediana de medianas, con varios tipos de entrada.
#include <algorithm>
#include <cstdio>
#include <random>
#include <string>
#include <utility>
#include <vector>

#include "select.h"

static void report(const std::string& name, const std::vector<int>& a) {
    long long comparisons = 0;
    int pivot = median_of_medians(a, comparisons);
    auto at_most = std::count_if(a.begin(), a.end(), [&](int x) { return x <= pivot; });
    auto at_least = std::count_if(a.begin(), a.end(), [&](int x) { return x >= pivot; });
    std::printf("%-20s %8ld %8ld\n", name.c_str(), static_cast<long>(at_most), static_cast<long>(at_least));
}

int main() {
    const int n = 1000;
    std::mt19937 rng(2026);
    std::vector<std::pair<std::string, std::vector<int>>> inputs;

    std::vector<int> random_values(n);
    for (int& x : random_values) {
        x = static_cast<int>(rng() % 1000000);
    }
    inputs.push_back({"al azar", random_values});

    std::vector<int> sorted(n);
    for (int i = 0; i < n; i++) {
        sorted[i] = i;
    }
    inputs.push_back({"ordenado", sorted});

    inputs.push_back({"todos iguales", std::vector<int>(n, 7)});

    std::vector<int> few(n);
    for (int& x : few) {
        x = static_cast<int>(rng() % 3);
    }
    inputs.push_back({"tres valores", few});

    inputs.push_back({"construida", tight_input(n)});

    std::printf("n = %d, cota 3n/10 - 3 = %d\n", n, 3 * n / 10 - 3);
    std::printf("%-20s %8s %8s\n", "entrada", "<= piv", ">= piv");
    for (const auto& [name, values] : inputs) {
        report(name, values);
    }
}

// Sigue quickselect con el primer elemento como pivote, vuelta por vuelta.
#include <cstddef>
#include <iostream>
#include <vector>

#include "select.h"

static void print(const std::vector<int>& v) {
    std::cout << "[";
    for (std::size_t i = 0; i < v.size(); i++) {
        std::cout << (i == 0 ? "" : " ") << v[i];
    }
    std::cout << "]";
}

int main() {
    std::vector<int> a = {6, 3, 9, 1, 7, 2, 8, 5, 4};
    std::size_t k = 4;
    long long comparisons = 0;
    while (true) {
        std::cout << "k = " << k << "  A = ";
        print(a);
        std::cout << "\n";
        if (a.size() == 1) {
            std::cout << "  un solo elemento: respuesta " << a[0] << "\n";
            break;
        }
        int pivot = a[0];
        Parts parts = partition3(a, pivot, comparisons);
        std::cout << "  pivote " << pivot << "  menores ";
        print(parts.less);
        std::cout << "  iguales ";
        print(parts.equal);
        std::cout << "  mayores ";
        print(parts.greater);
        std::cout << "\n";
        if (k < parts.less.size()) {
            a = parts.less;
        } else if (k < parts.less.size() + parts.equal.size()) {
            std::cout << "  k cae en los iguales: respuesta " << pivot << "\n";
            break;
        } else {
            k -= parts.less.size() + parts.equal.size();
            a = parts.greater;
        }
    }
    std::cout << comparisons << " comparaciones\n\n";

    std::cout << "arreglo ordenado, k = n - 1, primer elemento como pivote\n";
    std::cout << "        n   comparaciones   n(n+1)/2 - 1\n";
    for (long long n = 1000; n <= 8000; n *= 2) {
        std::vector<int> sorted(n);
        for (long long i = 0; i < n; i++) {
            sorted[i] = static_cast<int>(i);
        }
        long long count = 0;
        quickselect_first(sorted, n - 1, count);
        std::cout.width(9);
        std::cout << n;
        std::cout.width(16);
        std::cout << count;
        std::cout.width(15);
        std::cout << n * (n + 1) / 2 - 1 << "\n";
    }
}

// Calcula las recurrencias con grupos de 5 y de 3, tal como salen de contar elementos, y compara T(n)/n.
//   T5(n) = T5(techo(n/5)) + T5(piso(7n/10) + 3) + n    para n >= 80
//   T3(n) = T3(techo(n/3)) + T3(piso(2n/3) + 2) + n     para n >= 80
// Con n < 80 se toma T(n) = n.
#include <cstdio>
#include <vector>

int main() {
    const long long limit = 10240000;
    std::vector<long long> t5(limit + 1);
    std::vector<long long> t3(limit + 1);
    for (long long n = 0; n <= limit; n++) {
        if (n < 80) {
            t5[n] = n;
            t3[n] = n;
        } else {
            t5[n] = t5[(n + 4) / 5] + t5[7 * n / 10 + 3] + n;
            t3[n] = t3[(n + 2) / 3] + t3[2 * n / 3 + 2] + n;
        }
    }
    std::printf("%10s %10s %10s\n", "n", "T5(n)/n", "T3(n)/n");
    for (long long n = 10000; n <= limit; n *= 4) {
        std::printf("%10lld %10.2f %10.2f\n", n,
                    static_cast<double>(t5[n]) / static_cast<double>(n),
                    static_cast<double>(t3[n]) / static_cast<double>(n));
    }
    double worst = 0;
    for (long long n = 1; n <= limit; n++) {
        double ratio = static_cast<double>(t5[n]) / static_cast<double>(n);
        if (ratio > worst) {
            worst = ratio;
        }
    }
    std::printf("mayor T5(n)/n con n hasta %lld: %.2f\n", limit, worst);
}

// Mide selección con entradas al azar: comparaciones exactas y tiempo, con -O2 y sin sanitizers.
// Los tiempos dependen de la máquina. Son la mediana de 7 corridas.
#include <algorithm>
#include <cstdio>
#include <random>
#include <vector>

#include "bench.h"
#include "select.h"

int main() {
    std::printf("%9s | %11s %11s | %9s %9s %9s %9s\n", "n", "cmp.azar/n", "cmp.mdm/n", "sort ms", "nth ms", "azar ms", "mdm ms");
    for (int n = 125000; n <= 8000000; n *= 2) {
        std::mt19937 rng(2026);
        std::vector<int> a(n);
        for (int& x : a) {
            x = static_cast<int>(rng() >> 1);
        }
        std::size_t k = a.size() / 2;
        long long random_count = 0;
        long long linear_count = 0;
        std::mt19937 pivot_rng(7);
        int answer = quickselect_random(a, k, pivot_rng, random_count);
        int same = select_linear(a, k, linear_count);
        assert(answer == same);

        double sort_ms = 1000 * median_seconds(7, [&] {
            std::vector<int> copy = a;
            std::sort(copy.begin(), copy.end());
            keep(copy[k]);
        });
        double nth_ms = 1000 * median_seconds(7, [&] {
            std::vector<int> copy = a;
            std::nth_element(copy.begin(), copy.begin() + static_cast<long>(k), copy.end());
            keep(copy[k]);
        });
        double random_ms = 1000 * median_seconds(7, [&] {
            long long ignored = 0;
            keep(quickselect_random(a, k, pivot_rng, ignored));
        });
        double linear_ms = 1000 * median_seconds(7, [&] {
            long long ignored = 0;
            keep(select_linear(a, k, ignored));
        });
        std::printf("%9d | %11.2f %11.2f | %9.1f %9.1f %9.1f %9.1f\n", n,
                    static_cast<double>(random_count) / n, static_cast<double>(linear_count) / n,
                    sort_ms, nth_ms, random_ms, linear_ms);
    }
}

// Selección: el k-ésimo menor de un arreglo, con k contado desde 0.
// `comparisons` cuenta cada vez que un elemento se compara con un pivote o, dentro de un
// grupo, con otro elemento. Es un conteo exacto: no depende de la máquina.
#ifndef SELECT_H
#define SELECT_H

#include <algorithm>
#include <cassert>
#include <cstddef>
#include <random>
#include <utility>
#include <vector>

struct Parts {
    std::vector<int> less;
    std::vector<int> equal;
    std::vector<int> greater;
};

// Parte en tres: menores, iguales y mayores que el pivote. Cuesta una comparación por elemento.
inline Parts partition3(const std::vector<int>& a, int pivot, long long& comparisons) {
    Parts parts;
    for (int x : a) {
        comparisons++;
        if (x < pivot) {
            parts.less.push_back(x);
        } else if (x > pivot) {
            parts.greater.push_back(x);
        } else {
            parts.equal.push_back(x);
        }
    }
    return parts;
}

// Quickselect con una regla para elegir el pivote.
// Invariante: la respuesta buscada es el k-ésimo menor de `a`.
template <typename PivotRule>
int quickselect(std::vector<int> a, std::size_t k, PivotRule choose_pivot, long long& comparisons) {
    assert(k < a.size());
    while (a.size() > 1) {
        Parts parts = partition3(a, choose_pivot(a), comparisons);
        if (k < parts.less.size()) {
            a = std::move(parts.less);
        } else if (k < parts.less.size() + parts.equal.size()) {
            return parts.equal[0];
        } else {
            k -= parts.less.size() + parts.equal.size();
            a = std::move(parts.greater);
        }
    }
    return a[0];
}

inline int quickselect_first(std::vector<int> a, std::size_t k, long long& comparisons) {
    return quickselect(std::move(a), k, [](const std::vector<int>& v) { return v[0]; }, comparisons);
}

inline int quickselect_random(std::vector<int> a, std::size_t k, std::mt19937& rng, long long& comparisons) {
    auto rule = [&rng](const std::vector<int>& v) { return v[rng() % v.size()]; };
    return quickselect(std::move(a), k, rule, comparisons);
}

// Ordena un grupo pequeño por inserción. En un grupo de 5 son a lo sumo 10 comparaciones.
inline void sort_group(std::vector<int>& group, long long& comparisons) {
    for (std::size_t i = 1; i < group.size(); i++) {
        int value = group[i];
        std::size_t j = i;
        while (j > 0) {
            comparisons++;
            if (group[j - 1] <= value) {
                break;
            }
            group[j] = group[j - 1];
            j--;
        }
        group[j] = value;
    }
}

// Reparte en grupos de 5 por posición y devuelve la mediana inferior de cada uno. El último grupo puede ser más corto.
inline std::vector<int> group_medians(const std::vector<int>& a, long long& comparisons) {
    std::vector<int> medians;
    for (std::size_t start = 0; start < a.size(); start += 5) {
        std::size_t stop = std::min(start + 5, a.size());
        std::vector<int> group(a.begin() + start, a.begin() + stop);
        sort_group(group, comparisons);
        medians.push_back(group[(group.size() - 1) / 2]);
    }
    return medians;
}

int select_linear(std::vector<int> a, std::size_t k, long long& comparisons);

// El pivote de la mediana de medianas: la mediana inferior de las medianas de los grupos.
inline int median_of_medians(const std::vector<int>& a, long long& comparisons) {
    std::vector<int> medians = group_medians(a, comparisons);
    std::size_t middle = (medians.size() - 1) / 2;
    return select_linear(std::move(medians), middle, comparisons);
}

// Selección en tiempo lineal en el peor caso.
inline int select_linear(std::vector<int> a, std::size_t k, long long& comparisons) {
    assert(k < a.size());
    if (a.size() <= 5) {
        sort_group(a, comparisons);
        return a[k];
    }
    int pivot = median_of_medians(a, comparisons);
    Parts parts = partition3(a, pivot, comparisons);
    if (k < parts.less.size()) {
        return select_linear(std::move(parts.less), k, comparisons);
    }
    if (k < parts.less.size() + parts.equal.size()) {
        return pivot;
    }
    return select_linear(std::move(parts.greater), k - parts.less.size() - parts.equal.size(), comparisons);
}

// Una entrada que lleva la mediana de medianas casi hasta su cota: n múltiplo de 5, todos distintos.
// La mitad de los grupos (redondeada hacia arriba) tiene sus tres menores entre los valores
// más pequeños de todos, y el resto de los valores queda por encima del pivote.
inline std::vector<int> tight_input(int n) {
    assert(n % 5 == 0);
    int groups = n / 5;
    int low_groups = (groups + 1) / 2;
    int next_low = 0;
    int next_high = 3 * low_groups;
    std::vector<int> a;
    for (int g = 0; g < groups; g++) {
        int lows = g < low_groups ? 3 : 0;
        for (int i = 0; i < lows; i++) {
            a.push_back(next_low++);
        }
        for (int i = lows; i < 5; i++) {
            a.push_back(next_high++);
        }
    }
    return a;
}

#endif

// Cuatro experimentos con la franja del par más cercano, con fuerza bruta y un solo nivel de la división.
#include <cstdio>
#include <vector>

#include "strip.h"

static void print_step(const RootStep& step) {
    std::printf("linea x = %lld\n", step.line);
    std::printf("mejor izquierda = %lld, mejor derecha = %lld, delta^2 = %lld\n", step.left_best, step.right_best, step.delta2);
    std::printf("franja por y:");
    for (const StripPoint& p : step.strip) {
        std::printf(" (%lld,%lld)%c", p.point.x, p.point.y, p.from_left ? 'I' : 'D');
    }
    std::printf("\n");
}

int main() {
    std::printf("1. Ejemplo de los apuntes\n");
    std::vector<Point> example = {{1, 9}, {2, 3}, {3, 7}, {4, 1}, {5, 5}, {6, 4}, {7, 10}, {8, 8}, {9, 2}, {10, 6}};
    RootStep root = analyze_root(example);
    print_step(root);
    for (std::size_t i = 0; i < root.strip.size(); i++) {
        for (std::size_t j = i + 1; j < root.strip.size(); j++) {
            long long d = dist2(root.strip[i].point, root.strip[j].point);
            std::printf("  posiciones %zu y %zu: distancia^2 %lld%s\n", i, j, d, d < root.delta2 ? "  mejora" : "");
        }
    }
    std::printf("respuesta = %lld, fuerza bruta = %lld\n\n", root.best_in_strip, brute_force(example));

    std::printf("2. No basta con los puntos mas cercanos a la linea en x\n");
    std::vector<Point> near_line = {{-2, 0}, {-1, 100}, {1, -100}, {2, 0}};
    RootStep second = analyze_root(near_line);
    print_step(second);
    std::printf("distancia^2 entre los mas cercanos a la linea: %lld\n", dist2({-1, 100}, {1, -100}));
    std::printf("fuerza bruta = %lld\n\n", brute_force(near_line));

    std::printf("3. No basta con comparar con el siguiente\n");
    std::vector<Point> next_only = {{-9, 50}, {-1, 0}, {0, 2}, {8, 1}};
    RootStep third = analyze_root(next_only);
    print_step(third);
    std::printf("solo con el siguiente: %lld\n", best_with_window(third, 1));
    std::printf("con los 7 siguientes: %lld\n", best_with_window(third, 7));
    std::printf("fuerza bruta = %lld\n\n", brute_force(next_only));

    std::printf("4. La respuesta a 3 posiciones\n");
    std::vector<Point> three_apart = {{-99, 50}, {0, 0}, {1, 99}, {99, 51}};
    RootStep fourth = analyze_root(three_apart);
    print_step(fourth);
    for (std::size_t window = 1; window <= 3; window++) {
        std::printf("con los %zu siguientes: %lld\n", window, best_with_window(fourth, window));
    }
    std::printf("fuerza bruta = %lld\n", brute_force(three_apart));
}

// Herramientas para revisar el lema de la franja del par más cercano. Todo es fuerza bruta:
// estos programas comprueban afirmaciones de los apuntes con pocos puntos, y no sirven para
// muchos puntos porque cuestan Θ(n²).
#ifndef STRIP_H
#define STRIP_H

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

struct Point {
    long long x;
    long long y;
    bool operator==(const Point&) const = default;
};

inline long long dist2(const Point& p, const Point& q) {
    long long dx = p.x - q.x;
    long long dy = p.y - q.y;
    return dx * dx + dy * dy;
}

inline long long brute_force(const std::vector<Point>& points) {
    long long best = std::numeric_limits<long long>::max();
    for (std::size_t i = 0; i < points.size(); i++) {
        for (std::size_t j = i + 1; j < points.size(); j++) {
            best = std::min(best, dist2(points[i], points[j]));
        }
    }
    return best;
}

struct StripPoint {
    Point point;
    bool from_left;
};

struct RootStep {
    long long line;                // ℓ: la x del punto en la posición n div 2
    long long left_best;           // δ² de la mitad izquierda
    long long right_best;          // δ² de la mitad derecha
    long long delta2;              // δ²: el menor de los dos
    std::vector<StripPoint> strip;  // los puntos con (x - ℓ)² < δ², ordenados por y
    long long max_gap;              // el mayor j - i entre puntos de la franja con distancia² < δ²
    long long best_in_strip;        // la menor distancia² entre esos pares; δ² si no hay ninguno
    bool all_crossing;              // ¿todos esos pares tienen un punto de cada mitad?
};

// El primer nivel de la división, calculado por fuerza bruta. `points` ya debe venir ordenado por x
// y el orden entre puntos de igual y en la franja es el que traiga la entrada: `points` puede
// barajarse antes para probar otros desempates.
inline RootStep analyze_root(const std::vector<Point>& points) {
    assert(points.size() >= 4);
    std::size_t mid = points.size() / 2;
    std::vector<Point> left(points.begin(), points.begin() + static_cast<long>(mid));
    std::vector<Point> right(points.begin() + static_cast<long>(mid), points.end());
    RootStep step{};
    step.line = points[mid].x;
    step.left_best = brute_force(left);
    step.right_best = brute_force(right);
    step.delta2 = std::min(step.left_best, step.right_best);
    for (std::size_t i = 0; i < points.size(); i++) {
        long long dx = points[i].x - step.line;
        if (dx * dx < step.delta2) {
            step.strip.push_back({points[i], i < mid});
        }
    }
    std::stable_sort(step.strip.begin(), step.strip.end(),
                     [](const StripPoint& p, const StripPoint& q) { return p.point.y < q.point.y; });
    step.best_in_strip = step.delta2;
    step.all_crossing = true;
    for (std::size_t i = 0; i < step.strip.size(); i++) {
        for (std::size_t j = i + 1; j < step.strip.size(); j++) {
            long long d = dist2(step.strip[i].point, step.strip[j].point);
            if (d < step.delta2) {
                step.max_gap = std::max(step.max_gap, static_cast<long long>(j - i));
                step.best_in_strip = std::min(step.best_in_strip, d);
                if (step.strip[i].from_left == step.strip[j].from_left) {
                    step.all_crossing = false;
                }
            }
        }
    }
    return step;
}

// Compara cada punto de la franja solo con sus `window` siguientes.
inline long long best_with_window(const RootStep& step, std::size_t window) {
    long long best = step.delta2;
    for (std::size_t i = 0; i < step.strip.size(); i++) {
        for (std::size_t j = i + 1; j < step.strip.size() && j <= i + window; j++) {
            best = std::min(best, dist2(step.strip[i].point, step.strip[j].point));
        }
    }
    return best;
}

inline void sort_by_x(std::vector<Point>& points) {
    std::stable_sort(points.begin(), points.end(), [](const Point& p, const Point& q) { return p.x < q.x; });
}

#endif