EuclideanMST.h

/**
 * Author: Simon Lindholm
 * Date: 2018-07-19
 * License: CC0
 * Source: own work
 * Description: Euclidean minimum spanning tree.
 * Add an "index" member to Point if you need indices returned.
 * Usage:
 * Q q; q.ps = ...;
 * q.init(0,0,1 << 30); // if 0 <= x,y < 2^30
 * q.rec({&q}); // outputs MST candidates into 'out'
 * Time: O(n \log n \log 10^9); 3s for $n=10^5$
 * Status: fuzz-tested
 */
#pragma once

#include "Point.h"
typedef Point<ll> P;

int sign(ll x) { return x > 0 ? 2 : !x; }
vector<pair<P,P>> out;
struct Q {
	ll x,y,s; vector<P> ps; Q* ch = 0;
	void init(ll mx, ll my, ll ms) {
		x=mx, y=my, s=ms/2;
		if (sz(ps) <= 1) return;
		ch = new Q[4];
		trav(p,ps) ch[(p.x >= x+s)+2*(p.y >= y+s)].ps.push_back(p);
		rep(i,0,4) ch[i].init(x+i%2*s, y+i/2*s, s);
	}
	void rec(vector<Q*> A) {
		vector<Q*> near;
		trav(a, A) {
			if (ch && a->ch && max(abs(a->x-x), abs(a->y-y)) <= 2*s)
				rep(i,0,4) near.push_back(&a->ch[i]);
			else if (this < a) test(ps, a->ps);
		}
		if (ch) rep(i,0,4) ch[i].rec(near);
	}
	/// called to find the shortest path, OR when sz(b) == 1 to find all paths
	void test(vector<P>& a, vector<P>& b) {
		if (sz(a) > sz(b)) { test(b, a); return; }
		if (sz(a) < 4) {
			trav(x, a) trav(y, b) out.push_back({x, y}); return; }
		P p = a[rand() % sz(a)];
		pair<ll, P> best{LLONG_MAX, P{}};
		trav(y, b) best = min(best, {(p-y).dist2(), y});
		P q = best.second;
		vector<P> A[3], B[3] = {{q},{},{q}};
		trav(x, a) A[sign(p.cross(q, x))].push_back(x);
		trav(x, b) B[sign(p.cross(q, x))].push_back(x);
		trav(x, A[1]) out.emplace_back(x, q);
		test(A[0], B[0]); test(A[2], B[2]);
	}
};