#include <lemon/unionfind.h>

 #include <lemon/kruskal.h>

 void triangle()

 namespace _delaunay_bits {

 {

   struct Part {

     int prev, curr, next;

     Part(int p, int c, int n) : prev(p), curr(c), next(n) {} 

   };

   inline std::ostream& operator<<(std::ostream& os, const Part& part) {

     os << '(' << part.prev << ',' << part.curr << ',' << part.next << ')';

     return os;

   }

   inline double circle_point(const Point& p, const Point& q, const Point& r) {

     double a = p.x * (q.y - r.y) + q.x * (r.y - p.y) + r.x * (p.y - q.y);

     if (a == 0) return std::numeric_limits<double>::quiet_NaN();

     double d = (p.x * p.x + p.y * p.y) * (q.y - r.y) +

       (q.x * q.x + q.y * q.y) * (r.y - p.y) +

       (r.x * r.x + r.y * r.y) * (p.y - q.y);

     double e = (p.x * p.x + p.y * p.y) * (q.x - r.x) +

       (q.x * q.x + q.y * q.y) * (r.x - p.x) +

       (r.x * r.x + r.y * r.y) * (p.x - q.x);

     double f = (p.x * p.x + p.y * p.y) * (q.x * r.y - r.x * q.y) +

       (q.x * q.x + q.y * q.y) * (r.x * p.y - p.x * r.y) +

       (r.x * r.x + r.y * r.y) * (p.x * q.y - q.x * p.y);

     return d / (2 * a) + sqrt((d * d + e * e) / (4 * a * a) + f / a);

   }

   inline bool circle_form(const Point& p, const Point& q, const Point& r) {

     return rot90(q - p) * (r - q) < 0.0;

   }

   inline double intersection(const Point& p, const Point& q, double sx) {

     const double epsilon = 1e-8;

     if (p.x == q.x) return (p.y + q.y) / 2.0;

     if (sx < p.x + epsilon) return p.y;

     if (sx < q.x + epsilon) return q.y;

     double a = q.x - p.x;

     double b = (q.x - sx) * p.y - (p.x - sx) * q.y;    

     double d = (q.x - sx) * (p.x - sx) * (p - q).normSquare();

     return (b - sqrt(d)) / a;

   }

   struct YLess {

     YLess(const std::vector<Point>& points, double& sweep) 

       : _points(points), _sweep(sweep) {}

     bool operator()(const Part& l, const Part& r) const {

       const double epsilon = 1e-8;

       //      std::cerr << l << " vs " << r << std::endl;

       double lbx = l.prev != -1 ?

 	intersection(_points[l.prev], _points[l.curr], _sweep) :

 	- std::numeric_limits<double>::infinity();

       double rbx = r.prev != -1 ?

 	intersection(_points[r.prev], _points[r.curr], _sweep) :

 	- std::numeric_limits<double>::infinity();

       double lex = l.next != -1 ?

 	intersection(_points[l.curr], _points[l.next], _sweep) :

 	std::numeric_limits<double>::infinity();

       double rex = r.next != -1 ?

 	intersection(_points[r.curr], _points[r.next], _sweep) :

 	std::numeric_limits<double>::infinity();

       if (lbx > lex) std::swap(lbx, lex);

       if (rbx > rex) std::swap(rbx, rex);

       if (lex < epsilon + rex && lbx + epsilon < rex) return true;

       if (rex < epsilon + lex && rbx + epsilon < lex) return false;

       return lex < rex;

     }

     const std::vector<Point>& _points;

     double& _sweep;

   };

   struct BeachIt;

   typedef std::multimap<double, BeachIt> SpikeHeap;

   typedef std::multimap<Part, SpikeHeap::iterator, YLess> Beach;

   struct BeachIt {

     Beach::iterator it;

     BeachIt(Beach::iterator iter) : it(iter) {}

   };

 }

 inline void delaunay() {

   std::vector<Pedge> pedges;

   for(NodeIt n(g);n!=INVALID;++n) 

   using namespace _delaunay_bits;

     for(NodeIt m=++(NodeIt(n));m!=INVALID;++m)

   typedef _delaunay_bits::Part Part;

   typedef std::vector<std::pair<double, int> > SiteHeap;

   std::vector<Point> points;

   std::vector<Node> nodes;

   for (NodeIt it(g); it != INVALID; ++it) {

     nodes.push_back(it);

     points.push_back(coords[it]);

   }

   SiteHeap siteheap(points.size());

   double sweep;

   for (int i = 0; i < int(siteheap.size()); ++i) {

     siteheap[i] = std::make_pair(points[i].x, i);

   }

   std::sort(siteheap.begin(), siteheap.end());

   sweep = siteheap.front().first;

   YLess yless(points, sweep);

   Beach beach(yless);

   SpikeHeap spikeheap;

   std::set<std::pair<int, int> > edges;

   beach.insert(std::make_pair(Part(-1, siteheap[0].second, -1), 

 			      spikeheap.end()));

   int siteindex = 1;

   while (siteindex < int(points.size()) || !spikeheap.empty()) {

     SpikeHeap::iterator spit = spikeheap.begin();

     if (siteindex < int(points.size()) && 

 	(spit == spikeheap.end() || siteheap[siteindex].first < spit->first)) {

       int site = siteheap[siteindex].second;

       sweep = siteheap[siteindex].first;

       Beach::iterator bit = beach.upper_bound(Part(site, site, site));

       if (bit->second != spikeheap.end()) {

 	spikeheap.erase(bit->second);	

       }

       int prev = bit->first.prev;

       int curr = bit->first.curr;

       int next = bit->first.next;

       beach.erase(bit);

       SpikeHeap::iterator pit = spikeheap.end();

       if (prev != -1 && 

 	  circle_form(points[prev], points[curr], points[site])) {

 	double x = circle_point(points[prev], points[curr], points[site]);

 	pit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end())));

 	pit->second.it =

 	  beach.insert(std::make_pair(Part(prev, curr, site), pit));

       } else {

 	beach.insert(std::make_pair(Part(prev, curr, site), pit));

       }

       beach.insert(std::make_pair(Part(curr, site, curr), spikeheap.end()));

       SpikeHeap::iterator nit = spikeheap.end();

       if (next != -1 && 

 	  circle_form(points[site], points[curr],points[next])) {

 	double x = circle_point(points[site], points[curr], points[next]);

 	nit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end())));

 	nit->second.it =

 	  beach.insert(std::make_pair(Part(site, curr, next), nit));

       } else {

 	beach.insert(std::make_pair(Part(site, curr, next), nit));

       }

       ++siteindex;

     } else {

       sweep = spit->first;      

       Beach::iterator bit = spit->second.it;

       int prev = bit->first.prev;

       int curr = bit->first.curr;

       int next = bit->first.next;

 	Pedge p;

 	std::pair<int, int> edge;

 	p.a=n;

 	p.b=m;

 	edge = prev < curr ? 

 	p.len=(coords[m]-coords[n]).normSquare();

 	  std::make_pair(prev, curr) : std::make_pair(curr, prev);

 	pedges.push_back(p);

       }

 	if (edges.find(edge) == edges.end()) {

   std::sort(pedges.begin(),pedges.end(),pedgeLess);

 	  edges.insert(edge);

   for(std::vector<Pedge>::iterator pi=pedges.begin();pi!=pedges.end();++pi)

 	  g.addEdge(nodes[prev], nodes[curr]);

     {

 	  ++cnt;

       Line li(pi->a,pi->b);

 	}

       UEdgeIt e(g);

       for(;e!=INVALID && !cross(e,li);++e) ;

 	edge = curr < next ? 

       UEdge ne;

 	  std::make_pair(curr, next) : std::make_pair(next, curr);

       if(e==INVALID) {

 	ne=g.addEdge(pi->a,pi->b);

 	if (edges.find(edge) == edges.end()) {

 	edges.push_back(ne);

 	  edges.insert(edge);

 	cnt++;

 	  g.addEdge(nodes[curr], nodes[next]);

       }

 	  ++cnt;

     }

 	}

       }

       Beach::iterator pbit = bit; --pbit;

       int ppv = pbit->first.prev;

       Beach::iterator nbit = bit; ++nbit;

       int nnt = nbit->first.next;

       if (bit->second != spikeheap.end()) spikeheap.erase(bit->second);

       if (pbit->second != spikeheap.end()) spikeheap.erase(pbit->second);

       if (nbit->second != spikeheap.end()) spikeheap.erase(nbit->second);

       beach.erase(nbit);

       beach.erase(bit);

       beach.erase(pbit);

       SpikeHeap::iterator pit = spikeheap.end();

       if (ppv != -1 && ppv != next && 

 	  circle_form(points[ppv], points[prev], points[next])) {

 	double x = circle_point(points[ppv], points[prev], points[next]);

 	if (x < sweep) x = sweep;

 	pit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end())));

 	pit->second.it =

 	  beach.insert(std::make_pair(Part(ppv, prev, next), pit));

       } else {

 	beach.insert(std::make_pair(Part(ppv, prev, next), pit));

       }

       SpikeHeap::iterator nit = spikeheap.end();

       if (nnt != -1 && prev != nnt &&

 	  circle_form(points[prev], points[next], points[nnt])) {

 	double x = circle_point(points[prev], points[next], points[nnt]);

 	if (x < sweep) x = sweep;

 	nit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end())));

 	nit->second.it =

 	  beach.insert(std::make_pair(Part(prev, next, nnt), nit));

       } else {

 	beach.insert(std::make_pair(Part(prev, next, nnt), nit));

       }

     }

   }

   for (Beach::iterator it = beach.begin(); it != beach.end(); ++it) {

     int curr = it->first.curr;

     int next = it->first.next;

     if (next == -1) continue;

     std::pair<int, int> edge;

     edge = curr < next ? 

       std::make_pair(curr, next) : std::make_pair(next, curr);

     if (edges.find(edge) == edges.end()) {

       edges.insert(edge);

       g.addEdge(nodes[curr], nodes[next]);

       ++cnt;

     }

   }

   std::cout << T.realTime() << "s: Setting up the edges...\n";

   std::cout << T.realTime() << "s: Creating delaunay triangulation...\n";

   for(NodeIt n(g);n!=INVALID;++n) 

   delaunay();

     {

   std::cout << T.realTime() << "s: Calculating spanning tree...\n";

       for(NodeIt m=++(NodeIt(n));m!=INVALID;++m)

   LengthSquareMap<ListUGraph, ListUGraph::NodeMap<Point> > ls(g, coords);

 	{

   ListUGraph::UEdgeMap<bool> tree(g);

 	  Pedge p;

   kruskal(g, ls, tree);

 	  p.a=n;

   std::cout << T.realTime() << "s: Removing non tree edges...\n";

 	  p.b=m;

   std::vector<UEdge> remove;

 	  p.len=(coords[m]-coords[n]).normSquare();

   for (UEdgeIt e(g); e != INVALID; ++e) {

 	  pedges.push_back(p);

     if (!tree[e]) remove.push_back(e);

 	}

   }

       if(progress && en>=pr*double(N)/100) 

   for(int i = 0; i < int(remove.size()); ++i) {

 	{

     g.erase(remove[i]);

 	  std::cout << pr << "%  \r" << std::flush;

   }

 	  pr++;

 	}

     }

   std::cout << T.realTime() << "s: Sorting the edges...\n";

   std::sort(pedges.begin(),pedges.end(),pedgeLess);

   std::cout << T.realTime() << "s: Creating the tree...\n";

   ListUGraph::NodeMap<int> comp(g);

   UnionFind<ListUGraph::NodeMap<int> > uf(comp);

   for (NodeIt it(g); it != INVALID; ++it)

     uf.insert(it);  

   for(std::vector<Pedge>::iterator pi=pedges.begin();pi!=pedges.end();++pi)

     {

       if ( uf.join(pi->a,pi->b) ) {

 	g.addEdge(pi->a,pi->b);

 	en++;

 	if(en>=N-1)

 	  break;

       }

     }