计算几何

定义

struct Vec {
    double x, y;
    Vec(double x = 0, double y = 0) : x(x), y(y) {}
    Vec operator + (const Vec & a) const { return Vec(x + a.x, y + a.y); }
    Vec operator - (const Vec & a) const { return Vec(x - a.x, y - a.y); }
    Vec operator * (const double a) const { return Vec(x * a, y * a); }
    double operator ^ (const Vec & a) const { return x * a.x + y * a.y; }
    double operator * (const Vec & a) const { return x * a.y - y * a.x; }
};

typedef Vec Pt;

struct Line {
    Pt p; Vec v;
};

点与直线的位置关系

bool online(Pt a, Line b) {
    return fabs((a - b.p) * b.v) < eps ? true : false;
}

两直线的位置关系

int check(Line a, Line b) {
    return fabs(a.v * b.v) < eps ? fabs(a.v * (b.p - a.p)) < eps ? 0 : 1 : 2;
}

求两直线的交点

（下图来自这里）

可推出 \(\frac{S_{四边形ABCD}}{S_{△ACD}} = \frac{AB}{AP}\)。

Pt cross(Line a, Line b) {
    return a.p + a.v * ((b.v * (a.p - b.p)) / (a.v * b.v));
}

旋转矩阵

将向量 \((x,y)\) 绕原点逆时针旋转 \(\theta\) 角：

\[\begin{bmatrix} x'\\ y' \end{bmatrix}= \begin{bmatrix} \cos\theta & -\sin\theta\\ \sin\theta & \cos\theta \end{bmatrix} \begin{bmatrix} x\\ y \end{bmatrix} \]

即

\[x' = x\cos\theta - y\sin\theta\\ y' = x\sin\theta + y\cos\theta \]

一些性质

• 互相垂直的两向量的点积为 \(0\)。
• 在一条直线上的两向量的叉积为 \(0\)。
• 如果 \(\vec{a}\) 与 \(\vec{b}\) 的叉积为正，\(\vec{a}\) 与 \(\vec{c}\) 的叉积为负，则 \(\vec{b}\) 与 \(\vec{c}\) 分别在 \(\vec{a}\) 的两侧。
• 将向量的箭头向上，则左叉右为负，右叉左为正。

最小圆覆盖

#include <cmath>
#include <cstdio>
#include <algorithm>

const double eps = 1e-12;

struct Vec {
	double x, y;
	Vec(double x = 0, double y = 0) : x(x), y(y) {}
	Vec operator + (const Vec & a) const { return Vec(x + a.x, y + a.y); }
	Vec operator - (const Vec & a) const { return Vec(x - a.x, y - a.y); }
	Vec operator * (const double a) const { return Vec(x * a, y * a); }
	Vec operator / (const double a) const { return Vec(x / a, y / a); }
	double operator * (const Vec & a) const { return x * a.y - y * a.x; }
};
typedef Vec Pt;
struct Line { Pt p; Vec v; };
struct Cir { Pt p; double r; } c;
Pt p[100002]; int n;

double dis(Pt a, Pt b) {
	return (a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y);
}
Pt cross(Line a, Line b) {
	double s1 = (b.p - a.p) * (b.p + b.v - a.p), s2 = (b.p + b.v - a.p - a.v) * (b.p - a.p - a.v);
	return a.p + a.v * (s1 / (s1 + s2));
}
Vec rotate(Vec a) {
	return Vec(-a.y, a.x);
}
int main() {
	scanf("%d", &n);
	for (int i = 1; i <= n; ++i) scanf("%lf%lf", &p[i].x, &p[i].y);
	std::random_shuffle(p + 1, p + n + 1);
	c.p = p[1], c.r = 0;
	for (int i = 2; i <= n; ++i) {
		if (dis(c.p, p[i]) <= c.r + eps) continue;
		c.p = p[i], c.r = 0;
		for (int j = 1; j < i; ++j) {
			if (dis(c.p, p[j]) <= c.r + eps) continue;
			c.p = (p[i] + p[j]) / 2, c.r = dis(c.p, p[j]);
			for (int k = 1; k < j; ++k) {
				if (dis(c.p, p[k]) <= c.r + eps) continue;
				Line a, b;
				a.p = (p[i] + p[j]) / 2, b.p = (p[i] + p[k]) / 2;
				a.v = rotate(p[j] - p[i]), b.v = rotate(p[k] - p[i]);
				c.p = cross(a, b), c.r = dis(p[i], c.p);
			}
		}
	}
	printf("%.10lf\n%.10lf %.10lf", sqrt(c.r), c.p.x, c.p.y);
	return 0;
}

极角排序

此题为BZOJ 1132，题意是平面上有 \(n \le 3000\) 个点. 求出所有以这 \(n\) 个点为顶点的三角形的面积和。

利用叉积的分配律，用前缀和计算。

注意一开始的时候要将所有点按照 \(y < a.y || (y == a.y && x < a.x);\) 进行排序，否则后面的计算中会出现负数。

#include <cstdio>
#include <algorithm>

typedef long long LL;

struct Vec {
	LL x, y;
	Vec(LL x = 0, LL y = 0) : x(x), y(y) {}
	Vec operator + (const Vec & a) const { return Vec(x + a.x, y + a.y); }
	Vec operator - (const Vec & a) const { return Vec(x - a.x, y - a.y); }
	double operator * (const Vec & a) const { return x * a.y - y * a.x; }
	bool operator < (const Vec & a) const { return y < a.y || (y == a.y && x < a.x); }
};
typedef Vec Pt; Pt p[3002]; Vec v[3002]; LL ans;
LL read() {
	LL x = 0, f = 1; char c = getchar();
	while (c < '0' || c > '9') { if (c == '-') f = -1; c = getchar(); }
	while (c >= '0' && c <= '9') x = (x << 3) + (x << 1) + (c ^ 48), c = getchar();
	return x * f;
}
bool cmp(Vec a, Vec b) {
	double c = a * b;
	return c == 0 ? a.x > b.x : c > 0;
}
int main() {
	int n = read();
	for (int i = 1; i <= n; ++i) p[i].x = read(), p[i].y = read();
	std::sort(p + 1, p + n + 1);
	for (int i = 1; i <= n - 2; ++i) {
		for (int j = i + 1; j <= n; ++j) v[j] = p[j] - p[i];
		std::sort(v + i + 1, v + n + 1, cmp);
		for (int j = i + 2; j <= n; ++j) ans += v[i + 1] * v[j], v[i + 1] = v[i + 1] + v[j];
	}
	printf("%lld", ans >> 1), puts(ans & 1 ? ".5" : ".0");
	return 0;
}

计算多边形的面积

#include <cstdio>

struct Vec {
	double x, y;
	Vec(double x = 0, double y = 0) : x(x), y(y) {}
	Vec operator - (const Vec & a) const { return Vec(x - a.x, y - a.y); }
	double operator * (const Vec & a) const { return x * a.y - y * a.x; }
};
typedef Vec Pt; Pt p[102]; int n; double ans;

int main() {
	scanf("%d", &n);
	for (int i = 1; i <= n; ++i) scanf("%lf%lf", &p[i].x, &p[i].y);
	for (int i = 2; i < n; ++i) ans += (p[i] - p[1]) * (p[i + 1] - p[1]);
	if (ans < 0) ans = -ans;
	printf("%.0lf\n", ans / 2);
	return 0;
}

旋转卡壳

BZOJ 1069 最大土地面积

枚举对角线，然后求左右两侧的最大三角形面积。

#include <cstdio>
#include <algorithm>

const double eps = 1e-6;

struct Vec {
	double x, y;
	Vec(double x = 0, double y = 0) : x(x), y(y) {}
	Vec operator - (const Vec & a) const { return Vec(x - a.x, y - a.y); }
	double operator * (const Vec & a) const { return x * a.y - y * a.x; }
	bool operator < (const Vec & a) const { return y < a.y || (y == a.y && x < a.x); }
};
typedef Vec Pt; Pt p[2002], sta[2002]; int n, top; double ans;
bool cmp(Pt a, Pt b) {
	double c = (a - p[1]) * (b - p[1]);
	return -eps <= c && c <= eps ? a.x > b.x : c > eps;
}
int main() {
	scanf("%d", &n);
	for (int i = 1; i <= n; ++i) {
		scanf("%lf%lf", &p[i].x, &p[i].y);
		if (p[i] < p[1]) p[0] = p[1], p[1] = p[i], p[i] = p[0];
	}
	std::sort(p + 2, p + n + 1, cmp);
	sta[1] = p[1], sta[2] = p[2], top = 2;
	for (int i = 3; i <= n; ++i) {
		while (top > 1 && (p[i] - sta[top]) * (sta[top - 1] - sta[top]) <= eps) --top;
		sta[++top] = p[i];
	}
	while (top > 1 && (p[1] - sta[top]) * (sta[top - 1] - sta[top]) <= eps) --top;
	for (int i = 1; i <= top - 2; ++i) {
		int x = i + 1, y = i + 2 == top ? 1 : i + 3;
		for (int j = i + 2; j <= top; ++j) {
			while ((sta[x + 1] - sta[i]) * (sta[j] - sta[i]) > (sta[x] - sta[i]) * (sta[j] - sta[i])) if (++x > top) x = 1;
			while ((sta[j] - sta[i]) * (sta[y + 1] - sta[i]) > (sta[j] - sta[i]) * (sta[y] - sta[i])) if (++y > top) y = 1;
			ans = std::max(ans, (sta[x] - sta[i]) * (sta[j] - sta[i]) + (sta[j] - sta[i]) * (sta[y] - sta[i]));
		}
	}
	printf("%.3lf", ans / 2);
	return 0;
}

旋转卡壳还可以计算：

• 两凸边形的最远距离（线1卡住凸边形1的最低点，线2卡住凸边形2的最高点，然后同时逆时针旋转，最远距离一定是点与点）
• 两凸边形的最近距离（线1卡住凸边形1的最低点，线2卡住凸边形2的最高点，然后同时逆时针旋转，最近距离可能是点与点、点与边、边与边）
• 凸边形的最小面积外接矩形（利用性质：矩形的一条边一定与凸边形的一条边重合）
• 凸边形的最小周长外接矩形（利用性质：矩形的一条边一定与凸边形的一条边重合）