Closest points divide and conquer implementation nightmare (C++)

I'm trying to implement a divide and conquer closest points algorithm. As standard as it gets, yet my head is about to explode, because my code seems to (randomly) give incorrect answers. I wrote a random number generator using stl for testing purposes and the error I keep coming across is that every few runs the algorithm returns a pair that is clearly farther apart than the closest pair (separated by 1 unit of distance, which I've input manually).

Please forgive the global variables, but this is the 3rd time I've rewritten this and I just felt it was slightly easier to work with. Pastebin link for those who like to see more on their screens: http://ift.tt/1GdodKl

#include <vector>
#include <algorithm>
#include <cmath>
#include <iostream>
#include <ctime>
#include <random>

using namespace std;
using point = pair<int, int>;
double MAX = 1000000000.0;
double GLOBAL_minDist = MAX;
pair<point, point> GLOBAL_nearestPoints;

bool Xcmp ( const point &c1, const point &c2 ) {
  return ( c1.first < c2.first );
}

bool Ycmp ( const point &c1, const point &c2 ) {
  return ( c1.second < c2.second );
}

inline ostream& operator<< ( ostream& os, const point& p ) {
  return os << p.first << " " << p.second << "\n";
}

inline ostream& operator<< ( ostream& os, const vector<point> &points ) {
  for( auto itr = points.begin(); itr < points.end(); itr++ ) {
    os << *itr;
  }
  return os;
}

inline ostream& operator<< ( ostream& os, const pair<point, point> nearestPair ) {
  return os << static_cast<int> (nearestPair.second.first) << " " << static_cast<int> (nearestPair.second.second) << "\n"
            << static_cast<int> (nearestPair.first.first) << " " << static_cast<int> (nearestPair.first.second);
}

inline double distance( const point a, const point b ) {
  return sqrt( pow(( a.first - b.first ), 2 ) + pow( a.second - b.second, 2 ));
}

void bruteCP( const vector<point> &Xs ) {
  for( auto it = Xs.begin(); it < Xs.end() - 1; it++ ) {
    for( auto it2 = it + 1; it2 < Xs.end(); it2++ ) {
      double minDist = distance( *it, *it2 );
      if( minDist < GLOBAL_minDist ) {
        cout << minDist << "\n";
        GLOBAL_minDist = minDist;
        GLOBAL_nearestPoints = pair<point, point> ( *it, *it2 );
      }
    }
  }
}

void divConCP( const vector<point>& Xs, const vector<point>& Ys ) {
  int Xsize = Xs.size();
  if( Xsize <= 3 ) { bruteCP( Xs ); return; }

  int mid =  Xsize / 2;
  int median = Xs[mid].first;

  vector<point> leftYs;
  copy_if( Ys.begin(), Ys.end(), back_inserter(leftYs), [median](const point& point) 
          {return point.first <= median;} );
  vector<point>leftXs;
  leftXs.insert( leftXs.end(), Xs.begin(), Xs.begin() + mid );
  divConCP( leftXs, leftYs );

  vector<point> rightYs, rightXs;
  copy_if( Ys.begin(), Ys.end(), back_inserter(leftYs), [median](const point& point) 
          {return point.first > median;} );
  rightXs.insert( rightXs.end(), Xs.begin() + mid, Xs.end() );
  divConCP( rightXs, rightYs );

  vector<point> strip;
  copy_if( Ys.begin(), Ys.end(), back_inserter(strip), [median, GLOBAL_minDist](const point& point) 
          {return abs(point.first - median) < GLOBAL_minDist;} );

  //vector<point> uniqStrip;
  //unique_copy( strip.begin(), strip.end(), uniqStrip.begin() );

  for( auto itr = strip.begin(); itr < strip.end(); itr++ ) {
    for( auto itr2 = itr + 1; itr2 < strip.end() && (*itr2).second < GLOBAL_minDist; itr2++ ) {
      double minDist = distance( *itr, *itr2 );
      if( minDist < GLOBAL_minDist) {
        //cout << minDist << "\n";
        //cout << *itr << " " << *itr2 << "\n";
        GLOBAL_minDist = minDist;
        GLOBAL_nearestPoints = pair<point, point> ( *itr, *itr2 );
      }
    }
  }
}

int main() {
  int n, x, y;
  vector<point> Xs, Ys;
/*
  cin >> n;
  for( int i = 0; i < n; i++ ) {
    cin >> x;
    cin >> y;
   // x = -i;
   // y = -i;

    point xy( x, y );
    Xs.push_back( xy );
    Ys.push_back( xy );
  }
*/
    // DEBUG //

  n = 100000;
  srand(time(0));
  std::default_random_engine gen(time(0));
  std::uniform_int_distribution<int> dis(-20000000, 20000000);
  for( int i = 0; i < n - 2; i++ ) {
    x = dis( gen );
    y = dis( gen );
    //x = i;
    //y = i;
    point xy( x, y );
    Xs.push_back( xy );
    Ys.push_back( xy );
  }

    Xs.push_back( point( 20001, 20001 ));
    Ys.push_back( point( 20001, 20001 ));
    Xs.push_back( point( 20000, 20001 ));
    Ys.push_back( point( 20000, 20001 ));

  // DEBUG //

  sort( Xs.begin(), Xs.end(), Xcmp );
  sort( Ys.begin(), Ys.end(), Ycmp );

  divConCP( Xs, Ys );
  //bruteCP( Xs );
  cout << GLOBAL_minDist << "\n";
  cout << GLOBAL_nearestPoints << "\n";

}