#ifndef NAN2INF_EXPERIMENT1_ALGORITHMS_HPP
#define NAN2INF_EXPERIMENT1_ALGORITHMS_HPP

#include <algorithm>
#include <cstddef>
#include <cstdint>
#include <limits>
#include <optional>
#include <random>
#include <stdexcept>
#include <utility>
#include <vector>

namespace experiment1 {

// BEGIN closest
struct Point {
    std::int32_t x;
    std::int32_t y;
};

inline std::int64_t squared_distance(Point a, Point b) {
    const auto dx = std::int64_t(a.x) - std::int64_t(b.x);
    const auto dy = std::int64_t(a.y) - std::int64_t(b.y);
    return dx * dx + dy * dy;
}

inline bool by_y(Point a, Point b) {
    return a.y < b.y || (a.y == b.y && a.x < b.x);
}

// Entry: x order on [l,r). Return: y order on the same interval.
inline std::int64_t closest_range(std::vector<Point>& a,
                                 std::vector<Point>& scratch,
                                 std::size_t l, std::size_t r) {
    auto best = std::numeric_limits<std::int64_t>::max();
    if (r - l <= 3) {
        for (auto i = l; i < r; ++i)
            for (auto j = i + 1; j < r; ++j)
                best = std::min(best, squared_distance(a[i], a[j]));
        std::sort(a.begin() + l, a.begin() + r, by_y);
        return best;
    }
    const auto mid = l + (r - l) / 2;
    const auto px = a[mid].x;
    const auto left = closest_range(a, scratch, l, mid);
    const auto right = closest_range(a, scratch, mid, r);
    best = std::min(left, right);
    std::merge(a.begin() + l, a.begin() + mid,
               a.begin() + mid, a.begin() + r,
               scratch.begin() + l, by_y);
    std::copy(scratch.begin() + l, scratch.begin() + r, a.begin() + l);
    if (best == 0) return 0;

    std::size_t count = 0;
    for (auto i = l; i < r; ++i) {
        const auto dx = std::int64_t(a[i].x) - std::int64_t(px);
        if (dx * dx >= best) continue;
        for (auto j = count; j > 0; --j) {
            const auto dy = std::int64_t(a[i].y) - std::int64_t(scratch[j - 1].y);
            if (dy * dy >= best) break;
            best = std::min(best, squared_distance(a[i], scratch[j - 1]));
        }
        scratch[count++] = a[i];
    }
    return best;
}

inline std::optional<std::int64_t> closest_squared(std::vector<Point> points) {
    if (points.size() > 400000) throw std::invalid_argument("too many points");
    for (const auto p : points) {
        if (p.x < -10000000 || p.x > 10000000 ||
            p.y < -10000000 || p.y > 10000000)
            throw std::invalid_argument("coordinate outside [-10^7,10^7]");
    }
    if (points.size() < 2) return std::nullopt;
    std::sort(points.begin(), points.end(), [](Point a, Point b) {
        return a.x < b.x || (a.x == b.x && a.y < b.y);
    });
    std::vector<Point> scratch(points.size());
    return closest_range(points, scratch, 0, points.size());
}
// END closest

// BEGIN insertion
inline void Insertion(std::vector<int>& num) {
    for (std::size_t i = 1; i < num.size(); ++i) {
        const int value = num[i];
        auto j = i;
        while (j > 0 && value < num[j - 1]) {
            num[j] = num[j - 1];
            --j;
        }
        num[j] = value;
    }
}
// END insertion

// BEGIN selection
inline void Selection(std::vector<int>& num) {
    for (std::size_t i = 0; i + 1 < num.size(); ++i) {
        auto pos = i;
        for (auto j = i + 1; j < num.size(); ++j)
            if (num[j] < num[pos]) pos = j;
        std::swap(num[i], num[pos]);
    }
}
// END selection

// BEGIN shell
inline std::vector<std::size_t> shell_gaps(std::size_t n, int mode) {
    if (mode < 1 || mode > 3) throw std::invalid_argument("invalid Shell mode");
    const std::size_t multiplier = mode == 3 ? 3 : 2;
    const std::size_t add = mode == 1 ? 0 : 1;
    std::vector<std::size_t> gaps;
    for (std::size_t h = 1; h < n;) {
        gaps.push_back(h);
        if (h > (std::numeric_limits<std::size_t>::max() - add) / multiplier)
            break;
        h = multiplier * h + add;
    }
    return gaps;
}

inline void Shell(std::vector<int>& num, int mode) {
    const auto gaps = shell_gaps(num.size(), mode);
    for (auto it = gaps.rbegin(); it != gaps.rend(); ++it) {
        const auto h = *it;
        for (auto i = h; i < num.size(); ++i) {
            const int value = num[i];
            auto j = i;
            while (j >= h && value < num[j - h]) {
                num[j] = num[j - h];
                j -= h;
            }
            num[j] = value;
        }
    }
}

inline void Shell_1(std::vector<int>& num) { Shell(num, 1); }
inline void Shell_2(std::vector<int>& num) { Shell(num, 2); }
inline void Shell_3(std::vector<int>& num) { Shell(num, 3); }
// END shell

// BEGIN quick
inline std::size_t partition(std::vector<int>& num, std::size_t l,
                             std::size_t r, std::mt19937& rng) {
    std::uniform_int_distribution<std::size_t> choose(l, r - 1);
    std::swap(num[choose(rng)], num[r - 1]);
    const int pivot = num[r - 1];
    auto boundary = l;
    for (auto j = l; j + 1 < r; ++j)
        if (num[j] < pivot) std::swap(num[boundary++], num[j]);
    std::swap(num[boundary], num[r - 1]);
    return boundary;
}

inline void quick_range(std::vector<int>& num, std::size_t l, std::size_t r,
                        std::size_t depth, std::size_t& max_depth,
                        std::mt19937& rng) {
    while (l < r) {
        max_depth = std::max(max_depth, depth);
        if (r - l == 1) return;
        const auto p = partition(num, l, r, rng);
        // Recurse only into the smaller side; the larger side uses this frame.
        if (p - l < r - (p + 1)) {
            quick_range(num, l, p, depth + 1, max_depth, rng);
            l = p + 1;
        } else {
            quick_range(num, p + 1, r, depth + 1, max_depth, rng);
            r = p;
        }
        ++depth;
    }
}

inline std::size_t Quicksort(std::vector<int>& num,
                             std::uint32_t seed = 20260911U) {
    std::mt19937 rng(seed);
    std::size_t max_depth = 0;
    quick_range(num, 0, num.size(), 1, max_depth, rng);
    return max_depth;
}
// END quick

// BEGIN merge
inline void merge_range(std::vector<int>& num, std::vector<int>& scratch,
                        std::size_t l, std::size_t r) {
    if (r - l < 2) return;
    const auto mid = l + (r - l) / 2;
    merge_range(num, scratch, l, mid);
    merge_range(num, scratch, mid, r);
    auto i = l;
    auto j = mid;
    auto k = l;
    while (i < mid && j < r)
        scratch[k++] = num[i] <= num[j] ? num[i++] : num[j++];
    while (i < mid) scratch[k++] = num[i++];
    while (j < r) scratch[k++] = num[j++];
    std::copy(scratch.begin() + l, scratch.begin() + r, num.begin() + l);
}

inline void Mergesort(std::vector<int>& num) {
    if (num.size() < 2) return;
    std::vector<int> scratch(num.size());
    merge_range(num, scratch, 0, num.size());
}
// END merge

// BEGIN match
struct EqualBand { std::size_t begin, end; };

template<class Item, class RelativeOrder>
EqualBand cross_partition(std::vector<Item>& a, std::size_t l, std::size_t r,
                          RelativeOrder order) {
    auto low = l;
    auto i = l;
    auto high = r;
    while (i < high) {
        const int sign = order(a[i]);
        if (sign < 0) std::swap(a[low++], a[i++]);
        else if (sign > 0) std::swap(a[i], a[--high]);
        else ++i;
    }
    return {low, high};
}

template<class Nut, class Bolt, class Compare>
void match_range(std::vector<Nut>& nuts, std::vector<Bolt>& bolts,
                 std::size_t l, std::size_t r, Compare& compare,
                 std::mt19937& rng) {
    while (l < r) {
        if (r - l == 1) {
            if (compare(nuts[l], bolts[l]) != 0)
                throw std::invalid_argument("unmatched singleton");
            return;
        }
        std::uniform_int_distribution<std::size_t> choose(l, r - 1);
        const Nut pivot_nut = nuts[choose(rng)];
        const auto b = cross_partition(bolts, l, r, [&](const Bolt& bolt) {
            const int sign = compare(pivot_nut, bolt);
            return sign > 0 ? -1 : sign < 0 ? 1 : 0;
        });
        if (b.end - b.begin != 1)
            throw std::invalid_argument("missing or duplicate matching bolt");
        const Bolt pivot_bolt = bolts[b.begin];
        const auto n = cross_partition(nuts, l, r, [&](const Nut& nut) {
            return compare(nut, pivot_bolt);
        });
        if (n.end - n.begin != 1 || n.begin != b.begin)
            throw std::invalid_argument("duplicate nut or unequal size sets");
        const auto p = n.begin;
        if (p - l < r - (p + 1)) {
            match_range(nuts, bolts, l, p, compare, rng);
            l = p + 1;
        } else {
            match_range(nuts, bolts, p + 1, r, compare, rng);
            r = p;
        }
    }
}

template<class Nut, class Bolt, class Compare>
void MatchNutsBolts(std::vector<Nut>& nuts, std::vector<Bolt>& bolts,
                   Compare compare, std::uint32_t seed = 20260911U) {
    if (nuts.size() != bolts.size())
        throw std::invalid_argument("different counts");
    std::mt19937 rng(seed);
    match_range(nuts, bolts, 0, nuts.size(), compare, rng);
}
// END match

} // namespace experiment1
#endif