#include <iostream>
#include <fstream>
#include <string>
#include <cmath>
#include <vector>
#include <algorithm>
using namespace std;

struct point {
    double x, y;
    point() { x = y = 0.0; }
};

double dist_sq( const point& p1, const point& p2 ) {
    return (p1.x - p2.x) * (p1.x - p2.x) + (p1.y - p2.y) * (p1.y - p2.y);
}

// used to sort by x-coordinate
bool comp_x ( const point& p1,const point& p2 ) {
    return p1.x < p2.x;
}

// used to sort and merge by y-coordinate
bool comp_y ( const point& p1,const point& p2 ) {
    return p1.y < p2.y;
}

int min( int a, int b ) {
    if ( a < b )
        return a;
    else
        return b;
}

/* v is a vector of points. close1 and close2 are a closest pair
of points. The closest distance is returned.*/
double closest_pair( vector< point >& p, point& close1, point& close2 );

/* v is a vector of points. close1 and close2 are a closest pair
of points within v[i],...,v[j]. The closest distance is returned.*/
double rec_cl_pair( vector< point >&p, 
                    int i,
                    int j,
                    point& close1,
                    point& close2 );

int main() {
    string infile;
    cout << "Input file? ";
    cin >> infile;
    ifstream in( infile.c_str() );
    if ( !in.is_open() ) {
        cerr << "Unable to open " << infile << '\n';
        return 0;
    }

    vector< point > p;
    point ptemp;
    while ( in >> ptemp.x >> ptemp.y )
        p.push_back( ptemp );
    in.close();

    point close1, close2;
    cout << "delta = " << closest_pair( p, close1, close2 ) << '\n';
    cout << "closest pair: ";
    cout << '(' << close1.x << ',' << close1.y << ")  ";
    cout << '(' << close2.x << ',' << close2.y << ")\n";

    return 0;
}

double closest_pair( vector< point >& p, point& close1, point& close2 ) {
    point* start = &p[ 0 ];
    point* end = start + p.size();
    // partial_sort guarantees worst-case time n log n
    // sort by x-coordinate
    partial_sort( start, end, end, comp_x );
    return rec_cl_pair( p, 0, p.size() - 1, close1, close2 );
}

double rec_cl_pair( vector< point >&p, 
                    int i,
                    int j,
                    point& close1,
                    point& close2 ) {
    if ( j - i < 3 ) {
        // sort by y-coordinate
        sort( &p[ i ], &p[ i ] + ( j - i + 1 ), comp_y );
        if ( j - i == 1 ) { // two points
            close1 = p[ i ];
            close2 = p[ j ];
        }
        else { // three points
            int p1 = i, p2 = i + 1;
            if ( dist_sq( p[ i ], p[ i + 2 ] ) 
                < dist_sq( p[ p1 ], p[ p2 ] ) )
                p2 = i + 2;
            if ( dist_sq( p[ i + 1 ], p[ i + 2 ] ) 
                < dist_sq( p[ p1 ], p[ p2 ] ) ) {
                p1 = i + 1;
                p2 = i + 2;
            }
            close1 = p[ p1 ];
            close2 = p[ p2 ];
        }
        return sqrt( dist_sq( close1, close2 ) );
    }
    int k = ( i + j ) / 2;
    int l = p[ k ].x;
    point cl1L, cl2L, cl1R, cl2R;
    double deltaL, deltaR, deltasq, delta;
    deltaL = rec_cl_pair( p, i, k, cl1L, cl2L );
    deltaR = rec_cl_pair( p, k + 1, j, cl1R, cl2R );
    if ( deltaL < deltaR ) {
        deltasq = deltaL * deltaL;
        close1 = cl1L;
        close2 = cl2L;
    }
    else {
        deltasq = deltaR * deltaR;
        close1 = cl1R;
        close2 = cl2R;
    }
    delta = sqrt( deltasq );

    /* v is a temporary vector first used for merging,
    then for storing the points in the strip */
    vector< point > v( j - i + 1 );

    /* merge p[i],...,p[k] and p[k+1],...,p[j] by y-coordinate with
    the result in v */
    merge( &p[ i ], &p[ k ] + 1, &p[ k ] + 1, &p[ j ] + 1, 
        v.begin(), comp_y );

    // copy result of merge back into p
    int t;
    for ( t = i; t <= j; t++ )
        p[ t ] = v[ t - i ];

    // store points in strip in v
    t = -1;
    for ( k = i; k <= j; k++ )
        if ( p[ k ].x > l - delta && p[ k ].x < l + delta )
            v[ ++t ] = p[ k ];

    // points in strip are v[ 0 ],...,v[ t ]

    // look in strip for closer pairs
    int s;
    float dtemp;
    for ( k = 0; k < t; k++ )
        for ( s = k + 1; s <= min( t, k + 7 ); s++ ) {
            dtemp = dist_sq( v[ k ], v[ s ] );
            if ( dtemp < deltasq ) {
                deltasq = dtemp;
                close1 = v[ k ];
                close2 = v[ s ];
            }
        }
    return sqrt( deltasq );
}