#include<vector>
#include<cmath>
#include<algorithm>
#include<limits>


struct Vector2D {
    double x;
    double y;
    
    Vector2D operator-(Vector2D r) {
        return {x - r.x, y - r.y};
    }
    double operator*(Vector2D r) {
        return x * r.x + y * r.y;
    }
    
    double length() {
        return hypot(x, y);
    }
};


double closest_pair_divide_and_conquer(std::vector<Vector2D> points) {
    if(points.size() <= 1) {
        return std::numeric_limits<double>::infinity();
    } else if (points.size() == 2) {
        return (points[1] - points[0]).length();
    }
    
    std::sort(points.begin(), points.end(), [](Vector2D& left, Vector2D& right) {
        return left.x < right.x;
    });  // Sort points by x coordinates
    
    int mid_i = points.size() / 2;
    double min_dist = std::min(
        closest_pair_divide_and_conquer(std::vector<Vector2D>(points.begin(), points.begin()+mid_i)),
        closest_pair_divide_and_conquer(std::vector<Vector2D>(points.begin()+mid_i, points.end()))
    );  // Do the divide and conquer recursive call
    
    // Now we get the points that are in strips 
    double mid_x = (points[mid_i-1].x + points[mid_i].x) / 2.0f;
    std::vector<Vector2D> left_strip(mid_i);
    std::vector<Vector2D> right_strip(points.size()-mid_i);
    
    left_strip.resize(std::distance(
        left_strip.begin(),
        std::copy_if(points.begin(), points.begin()+mid_i, left_strip.begin(), [&](Vector2D point) {
            return mid_x - point.x < min_dist;
        })
    ));  // Populate left strip with points within min_dist
    right_strip.resize(std::distance(
        right_strip.begin(),
        std::copy_if(points.begin()+mid_i, points.end(), right_strip.begin(), [&](Vector2D point) {
            return point.x - mid_x < min_dist;
        })
    ));  // Populate right strip with points within min_dist
    
    // Now iterate over the strips and see if there is a closer pair across the mid_x line
    std::sort(left_strip.begin(), left_strip.end(), [](Vector2D& left, Vector2D& right) {
        return left.y < right.y;
    });  // Sort left strip by y coordinate
    std::sort(right_strip.begin(), right_strip.end(), [](Vector2D& left, Vector2D& right) {
        return left.y < right.y;
    });  // Sort right strip by y coordinate
    {
        // First point of right strip with y-coordinate within min_dist of the current point on left strip
        int begin = 0;
        // First point of right strip with y-coordinate more than min_dist above the current point on the left strip
        int end = 0;
        double strip_min_dist = min_dist;
        
        for(auto point : left_strip) {
            while(begin < (int)right_strip.size() && right_strip[begin].y <= point.y - min_dist) {
                begin++;
            }
            while(end < (int)right_strip.size() && right_strip[end].y < point.y + min_dist) {
                end++;
            }
            // Compare the point on the left strip to points on the right strip that can be within min_dist distance
            for(int i=begin; i<end; i++) {
                strip_min_dist = std::min(strip_min_dist, (right_strip[i] - point).length());
            }
        }
        min_dist = std::min(min_dist, strip_min_dist);
    }
    return min_dist;
}


int main() {
    std::vector<Vector2D> points = {
        {1, 2},
        {5, 6},
        {4, 8},
        {10, 4}
    };
    double result = closest_pair_divide_and_conquer(points);
    return 0;
}