求助 Guibas-Stolfi 分治算法
查看原帖
求助 Guibas-Stolfi 分治算法
448887
cancan123456楼主2022/10/2 11:14
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <algorithm>
using namespace std;
const int N = 100005;
const int M = 3 * N;
const double eps = 1e-9;
int n;
struct Point2D {
	double x, y;
	Point2D(double x_ = 0, double y_ = 0) {
		x = x_;
		y = y_;
	}
	void shake() {
		x += eps * (rand() - RAND_MAX / 2) / (RAND_MAX / 2);
		y += eps * (rand() - RAND_MAX / 2) / (RAND_MAX / 2);
	}
} point[N];
Point2D operator + (const Point2D & a, const Point2D & b) {
	return Point2D(a.x + b.x, a.y + b.y);
}
Point2D operator - (const Point2D & a, const Point2D & b) {
	return Point2D(a.x - b.x, a.y - b.y);
}
Point2D operator * (const Point2D & a, const double & b) {
	return Point2D(a.x * b, a.y * b);
}
double dot(const Point2D & a, const Point2D & b) {
	return a.x * b.x + a.y * b.y;
}
double cross(const Point2D & a, const Point2D & b) {
	return a.x * b.y - a.y * b.x;
}
double length(const Point2D & a) {
	return sqrt(dot(a, a));
}
bool operator < (const Point2D & a, const Point2D & b) {
	return a.x == b.x ? a.y < b.y : a.x < b.x;
}
struct Point3D {
	double x, y, z;
	Point3D(double x_ = 0, double y_ = 0, double z_ = 0) {
		x = x_;
		y = y_;
		z = z_;
	}
};
Point3D operator + (const Point3D & a, const Point3D & b) {
	return Point3D(a.x + b.x, a.y + b.y, a.z + b.z);
}
Point3D operator - (const Point3D & a, const Point3D & b) {
	return Point3D(a.x - b.x, a.y - b.y, a.z - b.z);
}
Point3D operator * (const Point3D & a, const double & b) {
	return Point3D(a.x * b, a.y * b, a.z * b);
}
double dot(const Point3D & a, const Point3D & b) {
	return a.x * b.x + a.y * b.y + a.z * b.z;
}
Point3D cross(const Point3D & a, const Point3D & b) {
	return Point3D(a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x);
}
Point3D magic_function(const Point2D & a) {
	return Point3D(a.x, a.y, dot(a, a));
}
bool in_circle(Point2D a, Point2D b, Point2D c, Point2D d) {
	if (cross(b - a, c - a) < 0) {
		swap(b, c);
	}
	Point3D a_ = magic_function(a);
	Point3D b_ = magic_function(b);
	Point3D c_ = magic_function(c);
	Point3D d_ = magic_function(d);
	Point3D normal = cross(b_ - a_, c_ - a_);
	return dot(normal, d_ - a_) <= -eps;
}
bool intersect(const Point2D & a, const Point2D & b, const Point2D & c, const Point2D & d) {
	return cross(c - a, b - a) * cross(b - a, d - a) >= eps && cross(a - d, c - d) * cross(c - d, b - d) >= eps;
}
struct Edge {
	int u, v, prev, next;
} edge[2 * M];
int head[N];
int cnt = 1;
void add_edge(int u, int v) {
	cnt++;
	edge[cnt].u = u;
	edge[cnt].v = v;
	edge[head[u]].prev = cnt;
	edge[cnt].next = head[u];
	head[u] = cnt;
}
void del_edge(int i) {
	if (edge[i].prev != 0) {
		edge[edge[i].prev].next = edge[i].next;
	}
	if (edge[i].next != 0) {
		edge[edge[i].next].prev = edge[i].prev;
	}
	if (head[edge[i].u] == i) {
		head[edge[i].u] = edge[i].next;
	}
}
int stack[N], top;
pair < int, int > get_base_LR_edge(int l, int r) {
	int mid = (l + r) / 2;
	stack[1] = l;
	stack[2] = l + 1;
	top = 2;
	for (int i = l + 2; i <= r; i++) {
		while (top > 1 && cross(point[stack[top]] - point[stack[top - 1]], point[i] - point[stack[top]]) < eps) {
			top--;
		}
		top++;
		stack[top] = i;
	}
	for (int i = l; i < r; i++) {
		if (stack[i] <= mid && stack[i + 1] > mid) {
			return make_pair(stack[i], stack[i + 1]);
		}
	}
	return make_pair(0, 0);
}
void solve(int l, int r) {
	if (r - l + 1 <= 3) {
		for (int i = l; i <= r; i++) {
			for (int j = l; j <= r; j++) {
				if (i != j) {
					add_edge(i, j);
				}
			}
		}
	} else {
		int mid = (l + r) / 2;
		solve(l, mid);
		solve(mid + 1, r);
		pair < int, int > now_LR_edge = get_base_LR_edge(l, r);
		while (true) {
			add_edge(now_LR_edge.first, now_LR_edge.second);
			int l_next = 0, r_next = 0;
			for (int v, i = head[now_LR_edge.first]; i != 0; i = edge[i].next) {
				v = edge[i].v;
				if (cross(point[now_LR_edge.second] - point[now_LR_edge.first], point[v] - point[now_LR_edge.first]) >= eps) {
					if (l_next == 0 || in_circle(point[now_LR_edge.first], point[now_LR_edge.second], point[l_next], point[v])) {
						l_next = v;
					}
				}
			}
			for (int v, i = head[now_LR_edge.second]; i != 0; i = edge[i].next) {
				v = edge[i].v;
				if (cross(point[v] - point[now_LR_edge.second], point[now_LR_edge.first] - point[now_LR_edge.second]) <= -eps) {
					if (r_next == 0 || in_circle(point[now_LR_edge.first], point[now_LR_edge.second], point[r_next], point[v])) {
						r_next = v;
					}
				}
			}
			if (l_next == 0 && r_next == 0) {
				return;
			}
			if (l_next != 0 && r_next != 0) {
				if (in_circle(point[now_LR_edge.first], point[now_LR_edge.second], point[l_next], point[r_next])) {
					r_next = 0;
				} else {
					l_next = 0;
				}
			}
			if (l_next != 0) {
				for (int v, i = head[now_LR_edge.first]; i != 0; i = edge[i].next) {
					v = edge[i].v;
					if (intersect(point[l_next], point[now_LR_edge.second], point[now_LR_edge.first], point[v])) {
						del_edge(i);
						del_edge(i ^ 1);
					}
				}
				now_LR_edge = make_pair(l_next, now_LR_edge.second);
			} else {
				for (int v, i = head[now_LR_edge.second]; i != 0; i = edge[i].next) {
					v = edge[i].v;
					if (intersect(point[r_next], point[now_LR_edge.first], point[now_LR_edge.second], point[v])) {
						del_edge(i);
						del_edge(i ^ 1);
					}
				}
				now_LR_edge = make_pair(now_LR_edge.first, r_next);
			}
		}
	}
}
struct KruskalEdge {
	int u, v;
	double w;
} res[2 * M];
bool operator < (const KruskalEdge & a, const KruskalEdge & b) {
	return a.w < b.w;
}
int fa[N];
int find(int x) {
	return x == fa[x] ? x : fa[x] = find(fa[x]);
}
void merge(int x, int y) {
	fa[find(x)] = find(y);
}
bool query(int x, int y) {
	return find(x) != find(y);
}
int main() {
	srand(19260817);
	int n;
	scanf("%d", &n);
	for (int i = 1; i <= n; i++) {
		scanf("%lf %lf", &point[i].x, &point[i].y);
		point[i].shake();
	}
	sort(point + 1, point + n + 1);
	solve(1, n);
	int cnt = 0;
	for (int u = 1; u <= n; u++) {
		for (int v, i = head[u]; i != 0; i = edge[i].next) {
			v = edge[i].v;
			cnt++;
			res[cnt].u = u;
			res[cnt].v = v;
			res[cnt].w = length(point[u] - point[v]);
		}
	}
	sort(res + 1, res + cnt + 1);
	for (int i = 1; i <= n; i++) {
		fa[i] = i;
	}
	double ans = 0;
	for (int i = 1; i <= cnt; i++) {
		if (query(res[i].u, res[i].v)) {
			merge(res[i].u, res[i].v);
			ans += res[i].w;
		}
	}
	printf("%lf", ans);
	return 0;
}
2022/10/2 11:14
加载中...