// Comprueba la desigualdad triangular en los 26 lugares de capitales.csv, con tres formas de guardar
// la distancia, y muestra por qué redondear al entero más cercano puede romperla.
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <iomanip>
#include <iostream>
#include <numbers>
#include <string>
#include <vector>

#include "geo.h"

namespace {

using geo::at;
using geo::Place;

struct Check {
    long long triples = 0;
    long long broken = 0;
    double worst = 0;  // cuánto le sobra a un lado, en metros, en el peor caso
};

template <typename Distance>
Check check_triples(const std::vector<Place>& places, Distance distance) {
    Check check;
    const int n = static_cast<int>(places.size());
    for (int a = 0; a < n; a++) {
        for (int b = 0; b < n; b++) {
            for (int c = 0; c < n; c++) {
                if (a == b || b == c || a == c) continue;
                check.triples++;
                const double excess = static_cast<double>(distance(places[at(a)], places[at(c)])) -
                                      static_cast<double>(distance(places[at(a)], places[at(b)])) -
                                      static_cast<double>(distance(places[at(b)], places[at(c)]));
                if (excess > 0) {
                    check.broken++;
                    check.worst = std::max(check.worst, excess);
                }
            }
        }
    }
    return check;
}

void show(const std::string& title, const Check& check) {
    std::cout << "  " << title << ": " << check.broken << " de " << check.triples << " ternas la rompen";
    if (check.broken > 0) std::cout << ", peor exceso " << check.worst << " m";
    std::cout << '\n';
}

}  // namespace

int main() {
    const auto places = geo::read_places("capitales.csv");
    const int n = static_cast<int>(places.size());
    std::cout << std::fixed << std::setprecision(3);

    double nearest = 1e12;
    double farthest = 0;
    int near_a = 0;
    int near_b = 0;
    int far_a = 0;
    int far_b = 0;
    double to_integer = 1;
    long double formula_gap = 0;
    for (int a = 0; a < n; a++) {
        for (int b = a + 1; b < n; b++) {
            const double value = geo::great_circle_m(places[at(a)], places[at(b)]);
            if (value < nearest) {
                nearest = value;
                near_a = a;
                near_b = b;
            }
            if (value > farthest) {
                farthest = value;
                far_a = a;
                far_b = b;
            }
            to_integer = std::min({to_integer, value - std::floor(value), std::ceil(value) - value});
            const double flat = static_cast<double>(geo::flat_m(places[at(a)], places[at(b)]));
            formula_gap = std::max<long double>(formula_gap, std::abs(flat - value) / value);
        }
    }
    std::cout << n << " lugares, " << n * (n - 1) / 2 << " pares\n";
    std::cout << "más cercanos: " << places[at(near_a)].name << " y " << places[at(near_b)].name << ", " << nearest << " m\n";
    std::cout << "más lejanos: " << places[at(far_a)].name << " y " << places[at(far_b)].name << ", " << farthest << " m\n";
    std::cout << std::setprecision(6) << "lo más cerca de un entero que queda una distancia: " << to_integer << " m\n"
              << std::setprecision(3);
    std::cout << "peor error relativo de tratar los grados como un plano: " << 100.0L * formula_gap << " %\n\n";

    std::cout << "Desigualdad triangular, d(a,c) <= d(a,b) + d(b,c), en todas las ternas:\n";
    show("double, en metros", check_triples(places, geo::great_circle_m));
    show("entero, redondeado hacia arriba", check_triples(places, geo::distance_m));
    show("entero, al más cercano", check_triples(places, [](const Place& a, const Place& b) {
        return std::llround(geo::great_circle_m(a, b));
    }));

    double tightest = 1e12;
    int ta = 0;
    int tb = 0;
    int tc = 0;
    for (int a = 0; a < n; a++) {
        for (int b = 0; b < n; b++) {
            for (int c = 0; c < n; c++) {
                if (a == b || b == c || a == c) continue;
                const double slack = geo::great_circle_m(places[at(a)], places[at(b)]) +
                                     geo::great_circle_m(places[at(b)], places[at(c)]) -
                                     geo::great_circle_m(places[at(a)], places[at(c)]);
                if (slack < tightest) {
                    tightest = slack;
                    ta = a;
                    tb = b;
                    tc = c;
                }
            }
        }
    }
    std::cout << "la terna más justa: " << places[at(ta)].name << ", " << places[at(tb)].name << ", "
              << places[at(tc)].name << ", le sobran " << tightest << " m\n\n";

    // Tres puntos sobre el ecuador, separados 1000,4 m: el redondeo al más cercano rompe la desigualdad.
    const double degrees = 1000.4 / (geo::EARTH_RADIUS_M * std::numbers::pi / 180.0);
    const Place a{"a", 0, 0};
    const Place b{"b", 0, degrees};
    const Place c{"c", 0, 2 * degrees};
    std::cout << "Tres puntos sobre el ecuador:\n";
    std::cout << "  reales: ab = " << geo::great_circle_m(a, b) << ", bc = " << geo::great_circle_m(b, c)
              << ", ac = " << geo::great_circle_m(a, c) << '\n';
    std::cout << "  al más cercano: " << std::llround(geo::great_circle_m(a, b)) << " + "
              << std::llround(geo::great_circle_m(b, c)) << " < " << std::llround(geo::great_circle_m(a, c)) << '\n';
    std::cout << "  hacia arriba: " << geo::distance_m(a, b) << " + " << geo::distance_m(b, c) << " >= "
              << geo::distance_m(a, c) << '\n';
}