计算几何模板

基本运算

浮点数的比较

using _T = long long;

constexpr double eps = 1e-8;
constexpr long double pi = 3.1415926535897932384l;
int sign(double x) {
    if (fabs(x) < eps)
        return 0;
    if (x < 0)
        return -1;
    return 1;
}
int cmp(double x, double y) {
    if (fabs(x - y) < eps)
        return 0;
    if (x < y)
        return -1;
    return 1;
}

$Point$ 类

template <class T>
struct point {
	T x, y;
	point(T x = 0, T y = 0) : x(x), y(y) {}
	bool operator<(const point& B) const {
		return x == B.x ? y < B.y : x < B.x;
	}
	bool operator==(const point& B) const {
		return !sign(x - B.x) && !sign(y - B.y);
	}
	point operator+(const point& B) const {
		return point(x + B.x, y + B.y);
	}
	point operator-(const point& B) const {
		return point(x - B.x, y - B.y);
	}
	point operator*(const T a) const {
		return point(x * a, y * a);
	}
	point operator/(const T a) const {
		return point(x / a, y / a);
	}
	T operator*(const point& B) const {
		return x * B.x + y * B.y;
	}
	T operator^(const point& B) const {
		return x * B.y - y * B.x;
	}
	point operator-() const {
		return point(-x, -y);
	}
	double angle() const {
		return atan2(this->y, this->x);
	}
	T length2() const {
		return x * x + y * y;
	}
	double length() const {
		return sqrt(length2());
	}
	point Unit() {	// 单位向量
		return *this / this->length();
	}
	point Normal() {  // 单位法向量
		return point(-y, x) / this->length();
	}
	point trunc(double r) {	 // 化为长度为r的向量
		double l = length();
		if (!sign(l))
			return *this;
		r /= l;
		return point(x * r, y * r);
	}
	point toleft() {
		return point(-y, x);
	}
	point toright() {
		return point(y, -x);
	}
	point rotate(double rad) {	// 逆时针旋转
		return point(x * cos(rad) - y * sin(rad), x * sin(rad) + y * cos(rad));
	}
	friend int relation(const point& a, const point& b, const point& c) {
		return sign((b - a) ^ (c - a));
	}
	friend double get_angle(const point& a, const point& b) {
		return acos((a * b) / a.length() / b.length());
	}
	friend T area(const point& a, const point& b, const point& c) {
		return abs((b - a) ^ (c - a));
	}
	friend T get_dist2(const point& a, const point& b) {
		return (a - b).length2();
	}
	friend double get_dist(const point& a, const point& b) {
		return sqrt(get_dist2(a, b));
	}
	friend double project(point a, point b, point c) {	// 求向量ac在向量ab上的投影长度
		return ((b - a) * (c - a)) / (b - a).length();
	}
	bool up() const {
		return y > 0 || (y == 0 && x >= 0);
	}
};
using Point = point<_T>;

极角排序

/* 按一,二,三,四象限的顺序(逆时针) */
bool argcmp(const Point& A, const Point& B) {
	if (A.up() ^ B.up())
		return A.up() > B.up();
	return (A ^ B) > 0;
}

点绕点旋转

//【将点A绕点B顺时针旋转theta(弧度)】
Point turn_PP(Point a, Point b, double theta) {
    double x = (a.x - b.x) * cos(theta) + (a.y - b.y) * sin(theta) + b.x;
    double y = -(a.x - b.x) * sin(theta) + (a.y - b.y) * cos(theta) + b.y;
    return Point(x, y);
}

复数表示

#include <complex>
using namespace std;
typedef complex<double> Point;
typedef Point Vector;  //复数定义向量后，自动拥有构造函数、加减法和数量积
const double eps = 1e-9;
int sign(double x) {
    if (fabs(x) < eps)
        return 0;
    if (x < 0)
        return -1;
    return 1;
}
double Length(Vector A) {
    return abs(A);
}
double Dot(Vector A, Vector B) {  //conj(a+bi)返回共轭复数a-bi
    return real(conj(A) * B);
}
double Cross(Vector A, Vector B) {
    return imag(conj(A) * B);
}
Vector Rotate(Vector A, double rad) {
    return A * exp(Point(0, rad));  //exp(p)返回以e为底复数的指数
}

直线

$Line$类

struct Line {
	Point p, v;
	double rad;
	Line() {}
	Line(Point P, Point V) : p(P), v(V) {
		rad = atan2(v.y, v.x);
	}
	Point get_point_in_line(double t) const {
		return p + v * t;
	}
	bool operator<(const Line& L) const {
		if (!cmp(rad, L.rad))
			return L.onLeft(p) > 0;
		return rad < L.rad;
	}
	int onLeft(const Point& a) const {	// 点a是否在直线的左边(>0:左  <0:右)
		return relation(p, p + v, a);
	}
	double distance(const Point& a) {  // 点a到直线的距离
		return abs((v ^ (a - p)) / v.length());
	}
	Point foot_point(const Point& a) {	// 点a在直线上的投影(垂足)
		return p + v * ((v * (a - p)) / (v * v));
	}
	Point symmetry_point(const Point& a) {	// 点a关于直线的对称点
		Point q = foot_point(a);
		return Point(2 * q.x - p.x, 2 * q.y - p.y);
	}
	friend Point getIntersection(const Line& a, const Line& b) {  // 两直线交点(不能平行)
		Point u = a.p - b.p;
		double t = (b.v ^ u) / (a.v ^ b.v);
		return a.get_point_in_line(t);
	}
	friend bool on_right(const Line& a, const Line& b, const Line& c) {	 // b,c的交点是否在直线a的右边
		Point o = getIntersection(b, c);
		return a.onLeft(o) <= 0;
	}
};

$Segment$类

struct Segment {
	Point a, b;
	Segment(Point a, Point b) : a(a), b(b) {}
	bool OnSegment(const Point& p) const {	// 点是否在线段上(包含端点)
		return sign((a - p) ^ (b - p)) == 0 && sign((a - p) * (b - p)) <= 0;
	}
	bool SegInterSection(const Segment& seg) const {  // 线段是否相交
		const Point &a1 = a, &a2 = b, &b1 = seg.a, &b2 = seg.b;
		auto c1 = (a2 - a1) ^ (b1 - a1), c2 = (a2 - a1) ^ (b2 - a1);
		auto c3 = (b2 - b1) ^ (a1 - b1), c4 = (b2 - b1) ^ (a2 - b1);
		return sign(c1) * sign(c2) <= 0 && sign(c3) * sign(c4) <= 0;
	}
	bool LineIntersection(const Line& L) const {  // 线段和直线是否相交(包含端点)
		return L.onLeft(a) * L.onLeft(b) <= 0;
	}
	double distancePoint(const Point& p) const {  // 点到线段的距离
		if (sign((p - a) * (b - a)) < 0 || sign((p - b) * (a - b)) < 0)
			return min(get_dist(a, p), get_dist(b, p));
		return Line(a, b - a).distance(p);
	}
	double distanceSeg(const Segment& seg) const {	// 线段到线段的距离
		if (SegInterSection(seg))
			return 0;
		return min({distancePoint(seg.a), distancePoint(seg.b), seg.distancePoint(a), seg.distancePoint(b)});
	}
};

光线反射

/* 光线v照射到平面l后反射 */
Point rotate(const Point& v, const Point& l) {
	Point res;
	Point E = l / l.length();  // 单位向量
	double d = E * v;
	return (E * 2 * d - v);
}

例题：两异面直线间最小距离的平方(以分式形式输出)

struct Point {
    int x, y, z;
    Point(int x = 0, int y = 0, int z = 0) : x(x), y(y), z(z) {}
};
Point read_Point() {
    Point p;
    scanf("%d%d%d", &p.x, &p.y, &p.z);
    return p;
}
int Dot(const Vector& A, const Vector& B) {
    return A.x * B.x + A.y * B.y + A.z * B.z;
}
int Length(const Vector& A) {
    return Dot(A, A);
}
Vector Cross(const Vector& A, const Vector& B) {
    return Vector(A.y * B.z - A.z * B.y, A.z * B.x - A.x * B.z, A.x * B.y - A.y * B.x);
}

LL gcd(LL a, LL b) {
    return b ? gcd(b, a % b) : a;
}
LL lcm(LL a, LL b) {
    return a / gcd(a, b) * b;
}
/* 分数类 */
struct Rat {
    LL a, b;
    Rat(LL a = 0) : a(a), b(1) {}
    Rat(LL x, LL y) : a(x), b(y) {
        if (b < 0)
            a = -a, b = -b;
        LL d = gcd(a, b);
        if (d < 0)
            d = -d;
        a /= d;
        b /= d;
    }
};

Rat operator+(const Rat& A, const Rat& B) {
    LL x = lcm(A.b, B.b);
    return Rat(A.a * (x / A.b) + B.a * (x / B.b), x);
}
Rat operator-(const Rat& A, const Rat& B) {
    return A + Rat(-B.a, B.b);
}
Rat operator*(const Rat& A, const Rat& B) {
    return Rat(A.a * B.a, A.b * B.b);
}
void updatemin(Rat& A, const Rat& B) {
    if (A.a * B.b > B.a * A.b)
        A.a = B.a, A.b = B.b;
}

// 点P到线段AB的距离的平方
Rat Rat_Distance2ToSegment(const Point& P, const Point& A, const Point& B) {
    if (A == B)
        return Length(P - A);
    Vector v1 = B - A, v2 = P - A, v3 = P - B;
    if (Dot(v1, v2) < 0)
        return Length(v2);
    else if (Dot(v1, v3) > 0)
        return Length(v3);
    else
        return Rat(Length(Cross(v1, v2)), Length(v1));
}
// 求异面直线 p1+s*u 和 p2+t*v 的公垂线对应的s。如果平行/重合，返回false
// p1p2 = p2-p1+t*v-s*u,该向量与u,v分别垂直,联立方程求解s,t
bool Rat_LineDistance3D(const Point& p1, const Vector& u, const Point& p2, const Vector& v, Rat& s) {
    LL b = (LL)Dot(u, u) * Dot(v, v) - (LL)Dot(u, v) * Dot(u, v); // 推出公式的分母恰好能判断是否平行/相交
    if (b == 0)
        return false;
    LL a = (LL)Dot(u, v) * Dot(v, p1 - p2) - (LL)Dot(v, v) * Dot(u, p1 - p2);
    s = Rat(a, b);
    return true;
}
// 求公垂线交点的坐标
void Rat_GetPointOnLine(const Point& A, const Point& B, const Rat& t, Rat& x, Rat& y, Rat& z) {
    x = Rat(A.x) + Rat(B.x - A.x) * t;
    y = Rat(A.y) + Rat(B.y - A.y) * t;
    z = Rat(A.z) + Rat(B.z - A.z) * t;
}
// 两垂足间的距离
Rat Rat_Distance2(const Rat& x1, const Rat& y1, const Rat& z1, const Rat& x2, const Rat& y2, const Rat& z2) {
    return (x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2) + (z1 - z2) * (z1 - z2);
}

int main() {
    LL maxx = 0;
    Point A = read_Point();
    Point B = read_Point();
    Point C = read_Point();
    Point D = read_Point();
    Rat s, t;
    bool ok = false;
    Rat ans = Rat(1000000000);
    if (Rat_LineDistance3D(A, B - A, C, D - C, s))
        if (s.a > 0 && s.a < s.b && Rat_LineDistance3D(C, D - C, A, B - A, t))
            if (t.a > 0 && t.a < t.b) {
                ok = true;  // 异面直线/相交直线
                Rat x1, y1, z1, x2, y2, z2;
                Rat_GetPointOnLine(A, B, s, x1, y1, z1);
                Rat_GetPointOnLine(C, D, t, x2, y2, z2);
                ans = Rat_Distance2(x1, y1, z1, x2, y2, z2);
            }
    if (!ok) {  // 平行直线/重合直线
        updatemin(ans, Rat_Distance2ToSegment(A, C, D));
        updatemin(ans, Rat_Distance2ToSegment(B, C, D));
        updatemin(ans, Rat_Distance2ToSegment(C, A, B));
        updatemin(ans, Rat_Distance2ToSegment(D, A, B));
    }
    printf("%lld %lld\n", ans.a, ans.b);
}

例题：平面最近点对

const int N = 200010;
Point q[N], tmp[N];
bool cmpx(const Point& a, const Point& b) {
	return a.x < b.x || (a.x == b.x && a.y < b.y);
}
bool cmpy(const Point& a, const Point& b) {
	return a.y < b.y || (a.y == b.y && a.x < b.x);
}
double Closest_Pair(int l, int r) {
	double res = 1e18;
	if (l == r)
		return res;
	if (l + 1 == r)
		return get_dist(q[l], q[r]);
	int mid = l + r >> 1;
	double d1 = Closest_Pair(l, mid);
	double d2 = Closest_Pair(mid + 1, r);
	res = min(d1, d2);
	int cnt = 0;
	for (int i = l; i <= r; i++) {
		if (fabs(q[mid].x - q[i].x) <= res)
			tmp[++cnt] = q[i];
	}
	sort(tmp + 1, tmp + 1 + cnt, cmpy);
	for (int i = 1; i <= cnt; i++)
		for (int j = i + 1; j <= cnt && tmp[j].y - tmp[i].y < res; j++)
			res = min(res, get_dist(tmp[i], tmp[j]));
	return res;
}
signed main() {
	int n;
	cin >> n;
	for (int i = 1; i <= n; i++)
		cin >> q[i].x >> q[i].y;
	sort(q + 1, q + 1 + n, cmpx);
	cout << fixed << setprecision(4) << Closest_Pair(1, n) << endl;
}

多边形

三角形

**三角形面积：**两边叉乘除以$2$的绝对值，海伦公式$S=sqrt(p(p-a)(p-b)(p-c))$，$S=(absinC)/2$

**重心：**到三角形三顶点距离的平方和最小，到三边距离之积最大

求任意多边形面积

double polygon_area(Point p[], int n) {
    double s = 0;
    for(int i = 1; i + 1 < n; i++)
        s += Cross(p[i] - p[0], p[i + 1] - p[0]);
    return s / 2;
}

判断点是否在多边形内

任意多边形

// (回转数法) 点在多边形内返回1, 点在多边形外返回0, 点在多边形上返回-1
int point_in_polygon(const Point& p, vector<Point>& poly) {
    int wn = 0;
    int n = poly.size();
    for (int i = 0; i < n; i++) {
        if (OnSegment(p, poly[i], poly[(i + 1) % n]))
            return -1;
        int k = sign(Cross(poly[(i + 1) % n] - poly[i], p - poly[i]));
        int d1 = sign(poly[i].y - p.y);
        int d2 = sign(poly[(i + 1) % n].y - p.y);
        if (k > 0 && d1 <= 0 && d2 > 0)
            wn++;
        if (k < 0 && d2 <= 0 && d1 > 0)
            wn--;
    }
    if (wn != 0)
        return 1;
    return 0;
}

凸多边形

bool inConvex(const Point& p, const vector<Point>& a) {   // a为凸包(按顺序排列)
	if (a.empty())
		return false;
	int l = 1, r = a.size() - 1;
	while (l <= r) {
		int mid = l + r >> 1;
		double ls = (a[mid] - a[0]) ^ (p - a[0]);
		double rs = (a[mid + 1] - a[0]) ^ (p - a[0]);
		if (ls >= 0 && rs <= 0) {
			if (((a[mid + 1] - a[mid]) ^ (p - a[mid])) >= 0)
				return true;
			return false;
		} else if (ls < 0) {
			r = mid - 1;
		} else {
			l = mid + 1;
		}
	}
	return false;
}

皮克定理

计算点阵中顶点在格点上的多边形面积（$a$表示多边形内部点数，$b$表示多边形边上点数）

$S=a+b/2-1$

圆

$Circle$类

struct Circle {
    Point p;
    double r;
    Circle(Point _p = Point(0, 0), double _r = 0) : p(_p), r(_r) {}
    // 三角形外接圆
    Circle(Point a, Point b, Point c) {
        Line u = Line({(a + b) / 2}, {(b - a).rotate(pi / 2)});
        Line v = Line({(a + c) / 2}, {(c - a).rotate(pi / 2)});
        p = getIntersection(u, v);
        r = get_dist(p, a);
    }
    // 三角形内切圆(bool t 只是为了与外接圆区别)
    Circle(Point a, Point b, Point c, bool t) {
        Line u, v;
        double m = atan2(b.y - a.y, b.x - a.x), n = atan2(c.y - a.y, c.x - a.x);
        u.p = a;
        u.v = u.p + Point(cos((n + m) / 2), sin((n + m) / 2));
        v.p = b;
        m = atan2(a.y - b.y, a.x - b.x), n = atan2(c.y - b.y, c.x - b.y);
        v.v = v.p + Point(cos((n + m) / 2), sin((n + m) / 2));
        p = getIntersection(u, v);
        r = Line(a, b).distance(p);
    }
    bool operator==(Circle v) {
        return (p == v.p) && sign(r - v.r) == 0;
    }
    bool operator<(Circle v) const {
        return ((p < v.p) || ((p == v.p) && sign(r - v.r) < 0));
    }
    double area() {
        return pi * r * r;
    }
    double length() {
        return 2 * pi * r;
    }
    // 点和圆的关系        -1圆内   0圆上   1圆外
    int relation(Point a) {
        double dist = get_dist(p, a);
        if (sign(dist - r) < 0)
            return -1;
        else if (sign(dist - r) == 0)
            return 0;
        return 1;
    }
    // 直线和圆的关系     -1相交   0相切    1相离
    int line_relation(Line v) {
        double dist = v.distance(p);
        if (sign(dist - r) < 0)
            return -1;
        else if (sign(dist - r) == 0)
            return 0;
        else
            return 1;
    }
    // 两圆的关系  5 相离   4 外切   3 相交   2内切    1 内含
    int circle_relation(Circle v) {
        double dist = get_dist(p, v.p);
        if (sign(dist - r - v.r) > 0)
            return 5;
        if (sign(dist - r - v.r) == 0)
            return 4;
        double l = fabs(r - v.r);
        if (sign(dist - r - v.r) < 0 && sign(dist - l) > 0)
            return 3;
        if (sign(dist - l) == 0)
            return 2;
        if (sign(dist - l) < 0)
            return 1;
        return -1;
    }
    // 求两个圆的交点，并返回交点个数
    int cross_circle(Circle v, Point& p1, Point& p2) {
        int rel = circle_relation(v);
        if (rel == 1 || rel == 5)
            return 0;
        double d = get_dist(p, v.p);
        double l = (d * d + r * r - v.r * v.r) / (d * 2);
        double h = sqrt(r * r - l * l);
        Point tmp = p + (v.p - p).trunc(l);
        p1 = tmp + ((v.p - p).toleft().trunc(h));
        p2 = tmp + ((v.p - p).toright().trunc(h));
        if (rel == 2 || rel == 4)
            return 1;
        return 2;
    }
    // 求直线和圆的交点，返回交点个数
    int cross_line(Line v, Point& p1, Point& p2) {
        if ((*this).line_relation(v) == 1)
            return 0;
        Point a = v.foot_point(p);
        double d = v.distance(p);
        d = sqrt(r * r - d * d);
        if (sign(d) == 0) {
            p1 = a, p2 = a;
            return 1;
        }
        p1 = a + v.v.trunc(d);
        p2 = a - v.v.trunc(d);
        return 2;
    }
    // 过一点作圆的切线(先判断点和圆的关系)
    int tangent_line(Point q, Line& u, Line& v) {
        int x = relation(q);
        if (x == -1)
            return 0;
        if (x == 0) {
            u = Line(q, (p - q).toleft());
            v = u;
            return 1;
        }
        double d = get_dist(p, q);
        double rad = asin(r / d);
        u = Line(q, (p - q).rotate(rad));
        v = Line(q, (p - q).rotate(-rad));
        return 2;
    }
    // 求两圆相交面积
    double circle_cross_area(Circle v) {
        int rel = circle_relation(v);
        if (rel >= 4)
            return 0;
        if (rel <= 2)
            return min(area(), v.area());
        double d = get_dist(p, v.p);
        double hf = (r + v.r + d) / 2;
        double ss = 2 * sqrt(hf * (hf - r) * (hf - v.r) * (hf - d));
        double a1 = acos((r * r + d * d - v.r * v.r) / (2.0 * r * d));
        a1 = a1 * r * r;
        double a2 = acos((v.r * v.r + d * d - r * r) / (2.0 * v.r * d));
        a2 = a2 * v.r * v.r;
        return a1 + a2 - ss;
    }
    // 得到过a,b两点, 半径为r1的两个圆
    friend int get_circle(Point a, Point b, double r1, Circle& c1, Circle& c2) {
        Circle x(a, r1), y(b, r1);
        int t = x.cross_circle(y, c1.p, c2.p);
        if (!t)
            return 0;
        c1.r = c2.r = r1;
        return t;
    }
};

凸包

$Convex$类

template <class T>
struct convex {
	vector<Point> q;
	convex() {}
	convex(vector<Point>& B) : q(B) {}
	convex(const convex& B) : q(B.q) {}
	convex& operator=(const convex& B) {
		q = B.q;
		return *this;
	}
	Point& operator[](int x) noexcept {
		return q[x];
	}
	int size() const {
		return q.size();
	}
	int ne(int x) const {
		return x == size() - 1 ? 0 : x + 1;
	}
	int pre(int x) const {
		return x == 0 ? size() - 1 : x - 1;
	}
	void init(vector<Point>& v) {
		sort(v.begin(), v.end());
		int n = v.size(), top = 0;
		vector<int> st(n + 10);
		for (int i = 0; i < n; i++) {
			while (top > 1 && sign((v[st[top]] - v[st[top - 1]]) ^ (v[i] - v[st[top - 1]])) <= 0)
				top--;
			st[++top] = i;
		}
		int k = top;
		for (int i = n - 2; i >= 0; i--) {
			while (top > k && sign((v[st[top]] - v[st[top - 1]]) ^ (v[i] - v[st[top - 1]])) <= 0)
				top--;
			st[++top] = i;
		}
		for (int i = 1; i < top; i++)
			q.push_back(v[st[i]]);
		return;
	}
	double get_length() {
		double res = 0;
		for (int i = 0; i < size(); i++)
			res += get_dist(q[i], q[ne(i)]);
		return res;
	}
	T get_area2() {
		T res = 0;
		for (int i = 0; i < size(); i++)
			res += (q[i] ^ q[ne(i)]);
		return abs(res);
	}
    Point getBaryCentre() const {  // 重心
		Point res(0, 0);
		double are = 0;
		const int sz = size();
		for (int i = 1; i < sz - 1; i++) {
			double tmp = (q[i] - q[0]) ^ (q[i + 1] - q[0]);
			if (!sign(tmp))
				continue;
			are += tmp;
			res.x += (q[0].x + q[i].x + q[i + 1].x) / 3 * tmp;
			res.y += (q[0].y + q[i].y + q[i + 1].y) / 3 * tmp;
		}
		if (sign(are))
			res = res / are;
		return res;
	}
    int is_in(const Point& p) const {  // -1:外  0:在边上  1:内
		int wn = 0;
		int n = q.size();
		for (int i = 0; i < n; i++) {
			if (Segment(q[i], q[(i + 1) % n]).OnSegment(p))
				return 0;
			int k = sign((q[(i + 1) % n] - q[i]) ^ (p - q[i]));
			int d1 = sign(q[i].y - p.y);
			int d2 = sign(q[(i + 1) % n].y - p.y);
			if (k > 0 && d1 <= 0 && d2 > 0)
				wn++;
			if (k < 0 && d2 <= 0 && d1 > 0)
				wn--;
		}
		if (wn != 0)
			return 1;
		return -1;
	}
    vector<T> sum;
	void get_sum() {
		vector<T> a(q.size());
		for (int i = 0; i < q.size(); i++)
			a[i] = q[pre(i)] ^ q[i];
		sum.resize(q.size());
		partial_sum(a.begin(), a.end(), sum.begin());
	}
	T query_sum(int l, int r) const {
		if (l <= r)
			return sum[r] - sum[l] + (q[r] ^ q[l]);
		return sum[size() - 1] - sum[l] + sum[r] + (q[r] ^ q[l]);
	}

	// 闵可夫斯基和
	convex operator+(const convex& B) const {
		const auto& a = this->q;
		const auto& b = B.q;
		int n = q.size(), m = b.size();
		Point sa = q[0], sb = b[0];
		for (int i = 0; i < n; i++) {
			if (a[i].y < sa.y || (a[i].y == sa.y && a[i].x < sa.x))
				sa = a[i];
		}
		for (int i = 0; i < m; i++) {
			if (b[i].y < sb.y || (b[i].y == sb.y && b[i].x < sb.x)) {
				sb = b[i];
			}
		}
		auto s = sa + sb;
		vector<Point> d(n + m);
		for (int i = 0; i < n; i++)
			d[i] = a[(i + 1) % n] - a[i];
		for (int i = 0; i < m; i++)
			d[n + i] = b[(i + 1) % m] - b[i];
		sort(d.begin(), d.end(), [&](const Point& A, const Point& B) {
			if (A.up() ^ B.up())
				return A.up() > B.up();
			return (A ^ B) > 0;
		});
		vector<Point> c(n + m);
		c[0] = s;
		for (int i = 0; i < n + m - 1; i++)
			c[i + 1] = c[i] + d[i];
		return c;
	}

	// 旋转卡壳
	template <class F>
	void rotating_calipres(const F& func) const {
		for (int i = 0, j = 1; i < q.size(); i++) {
			auto d = q[i], e = q[ne(i)];
			func(d, e, q[j]);
			while (area(d, e, q[j]) <= area(d, e, q[ne(j)])) {
				j = ne(j);
				func(d, e, q[j]);
			}
		}
	}

	// 凸包直径(平方)
	T diameter2() const {
		if (q.size() == 1)
			return 0;
		if (q.size() == 2)
			return get_dist2(q[0], q[1]);
		T ans = 0;
		auto func = [&](const Point& a, const Point& b, const Point& c) {
			ans = max({ans, get_dist2(a, c), get_dist2(b, c)});
		};
		rotating_calipres(func);
		return ans;
	}

	// 凸多边形关于某一方向的极点, 复杂度 O(logn)
	template <class F>
	int extreme(const F& dir) const {
		const auto check = [&](const int i) {
			return sign(dir(q[i]) ^ (q[ne(i)] - q[i])) >= 0;
		};
		const auto dir0 = dir(q[0]);
		const auto check0 = check(0);
		if (!check0 && check(this->size() - 1))
			return 0;
		const auto cmp = [&](const Point& v) {
			const int vi = &v - q.data();
			if (vi == 0)
				return 1;
			const auto checkv = check(vi);
			const auto t = sign(dir0 ^ (v - q[0]));
			if (vi == 1 && checkv == check0 && sign(dir0 ^ (v - q[0])) == 0)
				return 1;
			return checkv ^ (checkv == check0 && t <= 0);
		};
		return partition_point(q.begin(), q.end(), cmp) - q.begin();
	}

	// 过凸多边形外一点求凸多边形的切线, 返回切点下标, 复杂度 O(logn)
	// 必须保证点在多边形外
	pair<int, int> tangent(const Point& a) const {
		const int i = extreme([&](const Point& u) {
			return u - a;
		});
		const int j = extreme([&](const Point& u) {
			return a - u;
		});
		return {i, j};
	}

	// 求平行于给定直线的凸多边形的切线, 返回切点下标, 复杂度 O(logn)
	pair<int, int> tangent(const Line& a) const {
		const int i = extreme([&](...) {
			return a.v;
		});
		const int j = extreme([&](...) {
			return -a.v;
		});
		return {i, j};
	}
};
using Convex = convex<_T>;

最小最大三角形面积

pair<_T, _T> minmax_triangle(const vector<Point>& vec) {
	if (vec.size() <= 2)
		return {0, 0};
	const _T tmpans = abs((vec[0] - vec[1]) ^ (vec[0] - vec[2]));
	_T maxans = tmpans, minans = tmpans;

	vector<pair<int, int>> evt;
	evt.reserve(vec.size() * vec.size());

	for (signed i = 0; i < vec.size(); i++) {
		for (signed j = 0; j < vec.size(); j++) {
			if (i == j || vec[i] == vec[j])
				continue;
			evt.push_back({i, j});
		}
	}
	sort(evt.begin(), evt.end(), [&](const pair<int, int>& u, const pair<int, int>& v) {
		const Point du = vec[u.second] - vec[u.first], dv = vec[v.second] - vec[v.first];
		return argcmp({du.y, -du.x}, {dv.y, -dv.x});
	});
	vector<signed> vx(vec.size()), pos(vec.size());
	for (signed i = 0; i < vec.size(); i++)
		vx[i] = i;
	sort(vx.begin(), vx.end(), [&](int x, int y) {
		return vec[x] < vec[y];
	});
	for (signed i = 0; i < vx.size(); i++)
		pos[vx[i]] = i;
	for (auto [u, v] : evt) {
		signed i = pos[u], j = pos[v];
		if (i > j)
			swap(u, v), swap(i, j);
		const Point vecu = vec[u], vecv = vec[v];
		if (i > 0)
			minans = min(minans, abs((vec[vx[i - 1]] - vecu) ^ (vec[vx[i - 1]] - vecv)));
		if (j < vx.size() - 1)
			minans = min(minans, abs((vec[vx[j + 1]] - vecu) ^ (vec[vx[j + 1]] - vecv)));
		maxans = max({maxans, abs((vec[vx[0]] - vecu) ^ (vec[vx[0]] - vecv)), abs((vec[vx.back()] - vecu) ^ (vec[vx.back()] - vecv))});
		swap(vx[i], vx[j]);
		pos[u] = j, pos[v] = i;
	}
	return {minans, maxans};
}

动态凸包

/* 1.插入一个点   2.查询某个点是否被凸包包含(包括边界) */
struct Hull {
	set<Point> su, sd;
	bool query(set<Point>& s, Point u, int flag) {
		auto l = s.upper_bound(u), r = s.lower_bound(u);
		if (r == s.end() || l == s.begin())
			return false;
		l--;
		return relation(*l, u, *r) * flag >= 0;
	}
	void insert(set<Point>& s, Point u, int flag) {
		if (query(s, u, flag))
			return;
		auto it = s.insert(u).first;
		while (true) {
			auto mid = it;
			if (mid == s.begin())
				break;
			auto l = --mid;
			if (l == s.begin())
				break;
			l--;
			if (relation(*l, *mid, u) * flag <= 0)
				break;
			s.erase(mid);
		}
		while (true) {
			auto mid = it;
			mid++;
			if (mid == s.end())
				break;
			auto r = mid;
			r++;
			if (r == s.end())
				break;
			if (relation(u, *mid, *r) * flag <= 0)
				break;
			s.erase(mid);
		}
	}
	void insert(Point u) {
		insert(su, u, 1), insert(sd, u, -1);
	}
	int check(Point u) {
		return query(su, u, 1) && query(sd, u, -1);
	}
};

半平面交

Line line[N];           /* 求多个凸包的面积的交集 */
int n, m, idx;             // O(nlogn)
Point all[N], a[55];
int q[N];
double half_plane() {
    sort(line, line + idx);
    int hh = 0, tt = -1;
    for (int i = 0; i < idx; i++) {
        if (i && !cmp(line[i].rad, line[i - 1].rad))   // 斜率相同的直线只选靠左的
            continue;
        while (hh + 1 <= tt && on_right(line[i], line[q[tt - 1]], line[q[tt]]))
            tt--;
        while (hh + 1 <= tt && on_right(line[i], line[q[hh]], line[q[hh + 1]]))
            hh++;
        q[++tt] = i;
    }
    while (hh + 1 <= tt && on_right(line[q[hh]], line[q[tt - 1]], line[q[tt]]))
        tt--;
    while (hh + 1 <= tt && on_right(line[q[tt]], line[q[hh]], line[q[hh + 1]]))
        hh++;
    q[++tt] = q[hh];
    int k = 0;
    for (int i = hh; i < tt; i++)
        all[k++] = getIntersection(line[q[i]], line[q[i + 1]]);  // 求出最后围成的凸包的点
    double res = 0;
    for (int i = 1; i + 1 < k; i++)
        res += area(all[0], all[i], all[i + 1]);
    return res / 2;
}

/* 返回围成的凸包的顶点 */
vector<Point> half_plane(vector<Line> line) {
	sort(line.begin(), line.end());
	int hh = 0, tt = -1, num = line.size();
	vector<int> q(num + 10);
	for (int i = 0; i < num; i++) {
		if (i && !cmp(line[i].rad, line[i - 1].rad))
			continue;
		while (hh + 1 <= tt && on_right(line[i], line[q[tt - 1]], line[q[tt]]))
			tt--;
		while (hh + 1 <= tt && on_right(line[i], line[q[hh]], line[q[hh + 1]]))
			hh++;
		q[++tt] = i;
	}
	while (hh + 1 <= tt && on_right(line[q[hh]], line[q[tt - 1]], line[q[tt]]))
		tt--;
	while (hh + 1 <= tt && on_right(line[q[tt]], line[q[hh]], line[q[hh + 1]]))
		hh++;
	if (tt - hh + 1 <= 2)
		return {};
	q[++tt] = q[hh];
	vector<Point> res;
	for (int i = hh; i < tt; i++)
		res.push_back(getIntersection(line[q[i]], line[q[i + 1]]));
	return res;
}

$Jiangly$版

vector<Point> half_plane(vector<Line> line) {
	sort(line.begin(), line.end());
	deque<Line> ls;
	deque<Point> ps;
	for (auto l : line) {
		if (ls.empty()) {
			ls.push_back(l);
			continue;
		}
		while (!ps.empty() && l.onLeft(ps.back()) <= 0) {
			ps.pop_back();
			ls.pop_back();
		}
		while (!ps.empty() && l.onLeft(ps[0]) <= 0) {
			ps.pop_front();
			ls.pop_front();
		}
		if (!sign(l.v ^ ls.back().v)) {
			if ((l.v * ls.back().v) > 0) {
				if (l.onLeft(ls.back().p) <= 0) {
					assert(ls.size() == 1);
					ls[0] = l;
				}
				continue;
			}
			return {};
		}
		ps.push_back(getIntersection(ls.back(), l));
		ls.push_back(l);
	}
	while (!ps.empty() && ls[0].onLeft(ps.back()) <= 0) {
		ps.pop_back();
		ls.pop_back();
	}
	if (ls.size() <= 2)
		return {};
	ps.push_back(getIntersection(ls[0], ls.back()));

	return vector<Point>(ps.begin(), ps.end());
}

最小圆覆盖

/* 最小的圆覆盖所有的点 */
Point q[N];
Line get_line(Point a, Point b) {  // 求中垂线
    return Line({(a + b) / 2}, {(b - a).rotate(-pi / 2)});
}
Circle get_circle(Point a, Point b, Point c) {
    Line l = get_line(a, b), r = get_line(a, c);
    Point p = getIntersection(l, r);
    return {p, get_dist(p, a)};
}
signed main() {
    cin >> n;
    for (int i = 0; i < n; i++)
        cin >> q[i].x >> q[i].y;
    random_shuffle(q, q + n);
    Circle c = {q[0], 0};
    for (int i = 1; i < n; i++) {    // O(n)
        if (cmp(c.r, get_dist(c.p, q[i])) < 0) {
            c = {q[i], 0};
            for (int j = 0; j < i; j++) {
                if (cmp(c.r, get_dist(c.p, q[j])) < 0) {
                    c = {(q[i] + q[j]) / 2, get_dist(q[i], q[j]) / 2};
                    for (int k = 0; k < j; k++) {
                        if (cmp(c.r, get_dist(c.p, q[k])) < 0)
                            c = get_circle(q[i], q[j], q[k]);
                    }
                }
            }
        }
    }
}

/* 最小椭圆覆盖 */
cin >> n;
for (int i = 0; i < n; i++)
    cin >> q[i].x >> q[i].y;
double a, p;     // 椭圆逆时针旋转的角度, 长轴/短轴的值
cin >> a >> p;
for (int i = 0; i < n; i++) {
    q[i] = q[i].rotate(-a / 180 * pi);  // 顺时针旋转
    q[i].x /= p;         // 把坐标放缩, 将椭圆转化为求最小圆覆盖
}
// 接下来求最小圆覆盖即可

旋转卡壳

平面最远点对

const int N = 50010;
int n;
Point q[N];
int st[N], top;
void andrew() {   // O(NlogN)   瓶颈
	sort(q, q + n);
	for (int i = 0; i < n; i++) {
		while (top > 1 && relation(q[st[top - 1]], q[st[top]], q[i]) <= 0)
			top--;
		st[++top] = i;
	}
	int k = top;
	for (int i = n - 2; i >= 0; i--) {
		while (top > k && relation(q[st[top - 1]], q[st[top]], q[i]) <= 0)
			top--;
		st[++top] = i;
	}
	top--;  // 首尾是相同的, 尾部前移
}
int rotating_calipres() {   // O(N)
	if (top <= 2)
		return get_dist(q[0], q[n - 1]);
	int res = 0;
	for (int i = 1, j = 3; i <= top; i++) {
		auto d = q[st[i]], e = q[st[i + 1]];  // 此处不用取模, q[st[top+1]]就是q[st[1]]
		while (area(d, e, q[st[j]]) < area(d, e, q[st[j + 1]])) {
			j++;
			j = j > top ? 1 : j;
		}
		res = max({res, get_dist(d, q[st[j]]), get_dist(e, q[st[j]])});
	}
	return res;
}
signed main() {
	cin.tie(nullptr)->sync_with_stdio(false);
	cin >> n;
	for (int i = 0; i < n; i++)
		cin >> q[i].x >> q[i].y;
	andrew();
	cout << rotating_calipres() << endl;

	return (0 ^ 0);
}

最小矩形覆盖

$①$  每条边上都有凸包上的点

$②$  矩形上的边与凸包上的某条边共线

const int N = 50010;
const double INF = 1e20, pi = acos(-1);
int n;
Point q[N], ans[5];
int st[N], top;
double min_area = INF;
void andrew() {
	sort(q, q + n);
	for (int i = 0; i < n; i++) {
		while (top > 1 && relation(q[st[top - 1]], q[st[top]], q[i]) <= 0)
			top--;
		st[++top] = i;
	}
	int k = top;
	for (int i = n - 2; i >= 0; i--) {
		while (top > k && relation(q[st[top - 1]], q[st[top]], q[i]) <= 0)
			top--;
		st[++top] = i;
	}
	top--;
}
void rotating_calipers() {
	for (int i = 1, a = 3, b = 2, c = 3; i <= top; i++) {
		auto d = q[st[i]], e = q[st[i + 1]];
		while (cmp(area(d, e, q[st[a]]), area(d, e, q[st[a + 1]])) < 0)
			a = a % top + 1;
		while (cmp(project(d, e, q[st[b]]), project(d, e, q[st[b + 1]])) < 0)
			b = b % top + 1;
		if (i == 1)
			c = a;
		while (cmp(project(d, e, q[st[c]]), project(d, e, q[st[c + 1]])) > 0)
			c = c % top + 1;
		auto x = q[st[a]], y = q[st[b]], z = q[st[c]];
		auto h = area(d, e, x) / (e - d).length();
		auto w = (y - z) * (e - d) / (e - d).length();
		if (h * w < min_area) {
			min_area = h * w;
			ans[0] = d + (e - d).Unit() * project(d, e, y);
			ans[3] = d + (e - d).Unit() * project(d, e, z);
			auto u = (e - d).rotate(pi / 2).Unit();
			ans[1] = ans[0] + u * h;
			ans[2] = ans[3] + u * h;
		}
	}
}
signed main() {
	cin.tie(nullptr)->sync_with_stdio(false);
	cin >> n;
	for (int i = 0; i < n; i++)
		cin >> q[i].x >> q[i].y;
	andrew();
	rotating_calipers();

	int k = 0;
	for (int i = 1; i < 4; i++)
		if (cmp(ans[i].y, ans[k].y) < 0 || !cmp(ans[i].y, ans[k].y) && cmp(ans[i].x, ans[k].x) < 0)
			k = i;

	cout << fixed << setprecision(5) << min_area << endl;
	for (int i = 0; i < 4; i++, k++) {
		auto x = ans[k % 4].x, y = ans[k % 4].y;
		if (!sign(x))
			x = 0;
		if (!sign(y))
			y = 0;
		cout << x << ' ' << y << endl;
	}

	return (0 ^ 0);
}

/* 使用Convex类 */
template <class F>
void rotating_calipres(const F& func) const {
    const auto area = [](const Point& u, const Point& v, const Point& w) {
        return (v - u) ^ (w - u);
    };
    for (int i = 0, j = 1, r = 1, l = -1; i < q.size(); i++) {
        int nxti = ne(i);
        const Point &pn1 = q[i], &pn2 = ((q[nxti] - q[i]).toleft() + q[i]);
        while (area(q[i], q[nxti], q[ne(j)]) >= area(q[i], q[nxti], q[j])) {
            j = ne(j);
        }
        if (l == -1)
            l = j;
        while (area(pn1, pn2, q[ne(r)]) <= area(pn1, pn2, q[r])) {
            r = ne(r);
        }
        while (area(pn1, pn2, q[ne(l)]) >= area(pn1, pn2, q[l])) {
            l = ne(l);
        }
        func(q[i], q[nxti], q[j], pn1, pn2, q[l], q[r]);
    }
}

T rectarea() const {
    if (q.size() <= 2)
        return 0;
    T ans = 2e9;
    auto func = [&](const Point& a, const Point& b, const Point& c, const Point& d, const Point& e, const Point& f, const Point& g) {
        Line l1(a, b - a), l2(d, e - d);
        ans = min(ans, abs(l1.distance(c)) * (abs(l2.distance(f)) + abs(l2.distance(g))));
    };
    rotating_calipres(func);
    return ans;
}

非严格第二大多边形面积

给$n$个点，求围成的面积非严格第二大的多边形的面积（不一定是凸多边形，边不能交叉）

最大的面积显然是凸包，次大就枚举删点和删边。

删边时要在内部维护一个凸包，然后枚举外凸包的边和内凸包的点进行旋转卡壳

struct Point {
	ll x, y;
	Point(ll x = 0, ll y = 0) : x(x), y(y) {}
	friend int relation(Point a, Point b, Point c) {
		ll flag = ((b - a) ^ (c - a));
		if (flag > 0)
			return 1;
		else if (flag < 0)
			return -1;
		return 0;
	}
	friend ll area(Point a, Point b, Point c) {
		return abs((b - a) ^ (c - a));
	}
};
const int N = 1e5 + 10;
Point q[N], t[N];
int n, m, st1[N], top1, st2[N], top2;
void andrew(Point q[], int st[], int n, int& top) {
	sort(q, q + n);
	for (int i = 0; i < n; i++) {
		while (top > 1 && relation(q[st[top - 1]], q[st[top]], q[i]) < 0)
			top--;
		st[++top] = i;
	}
	int k = top;
	for (int i = n - 2; i >= 0; i--) {
		while (top > k && relation(q[st[top - 1]], q[st[top]], q[i]) < 0)
			top--;
		st[++top] = i;
	}
	if (n > 1)
		top--;
}
signed main() {
	cin.tie(nullptr)->sync_with_stdio(false);
	cin >> n;
	for (int i = 0; i < n; i++)
		cin >> q[i].x >> q[i].y;
	andrew(q, st1, n, top1);  // 外凸包
	set<int> out;
	for (int i = 1; i <= top1; i++)
		out.insert(st1[i]);
	for (int i = 0; i < n; i++) {
		if (!out.count(i))
			t[m++] = q[i];
	}
	andrew(t, st2, m, top2);  // 内凸包

	ll max_area = 0;
	for (int i = 1; i <= top1; i++)
		max_area += (q[st1[i]] ^ q[st1[i + 1]]);  // 外凸包的面积, 即最大面积

	ll res = LONG_LONG_MAX;			   // 切掉的最小面积
	for (int i = 1; i <= top1; i++) {  // 枚举删掉外凸包上的每个点
		int l = i - 1, r = i + 1;
		if (!l)
			l = top1;
		if (r > top1)
			r = 1;
		res = min(res, area(q[st1[l]], q[st1[i]], q[st1[r]]));
	}

	if (top2) {
		int beg = 1;  // 进行旋转卡壳时, 内凸包的起点
		ll minx = area(q[st1[1]], q[st1[2]], t[st2[1]]);
		for (int i = 2; i <= top2; i++) {
			ll tmp = area(q[st1[1]], q[st1[2]], t[st2[i]]);
			if (tmp < minx) {
				minx = tmp;
				beg = i;
			}
		}

		res = min(res, minx);
		for (int i = 2, j = beg; i <= top1; i++) {
			auto d = q[st1[i]], e = q[st1[i + 1]];
			while (area(d, e, t[st2[j]]) > area(d, e, t[st2[j + 1]]))
				j = j % top2 + 1;
			res = min(res, area(d, e, t[st2[j]]));
		}
	}

	cout << max_area - res << endl;

	return (0 ^ 0);
}

三角剖分

圆的圆心在原点，半径为$R$，一个有$n$个顶点的简单多边形，求圆与多边形的面积并。

/* 读入时相邻两行所描述的顶点在多边形中也是相邻的,但顺/逆时针不确定 */
const int N = 55;
double R;
int n;
Point q[N], r;  // r为圆心(0, 0)
// 点p是否在线段ab上
bool on_segment(Point p, Point a, Point b) {
	return !sign((a - p) ^ (b - p)) && sign((a - p) * (b - p)) < 0;
}
// 求圆和线段的交点, 返回圆心到线段最短距离
double get_circle_line_intersection(Point a, Point b, Point& pa, Point& pb) {
	Point e = getIntersection((Line){a, b - a}, (Line){r, (b - a).rotate(-pi / 2)});
	double mind = get_dist(r, e);
	if (!on_segment(e, a, b))
		mind = min(get_dist(r, a), get_dist(r, b));
	if (cmp(R, mind) <= 0)
		return mind;
	double len = sqrt(R * R - get_dist(r, e) * get_dist(r, e));
	pa = e + (a - b).Unit() * len;
	pb = e + (b - a).Unit() * len;
	return mind;
}
// 求oab扇形的面积(a,b在圆外)
double get_sector(Point a, Point b) {
	double angle = acos((a * b) / a.length() / b.length());
	if (sign(a ^ b) < 0)
		angle = -angle;
	return R * R * angle / 2;
}
// 求圆O和三角形OAB的面积并
double get_circle_triangle_area(const Point& a, const Point& b) {
	double da = get_dist(a, r), db = get_dist(b, r);
	if (cmp(R, da) >= 0 && cmp(R, db) >= 0)
		return (a ^ b) / 2;
	if (!sign(a ^ b))
		return 0;
	Point pa, pb;
	double mind = get_circle_line_intersection(a, b, pa, pb);
	if (cmp(R, mind) <= 0)
		return get_sector(a, b);
	if (cmp(R, da) >= 0)
		return (a ^ pb) / 2 + get_sector(pb, b);
	if (cmp(R, db) >= 0)
		return get_sector(a, pa) + (pa ^ b) / 2;
	return get_sector(a, pa) + (pa ^ pb) / 2 + get_sector(pb, b);
}
void solve() {
	double res = 0;
	for (int i = 0; i < n; i++)
		res += get_circle_triangle_area(q[i], q[(i + 1) % n]);
	cout << fixed << setprecision(2) << fabs(res) << endl;
}
signed main() {
	cin.tie(nullptr)->sync_with_stdio(false);
	while (cin >> R >> n) {
		for (int i = 0; i < n; i++)
			cin >> q[i].x >> q[i].y;
		solve();
	}

	return (0 ^ 0);
}

闵可夫斯基和

从原点向图形$A$内部的每一个点做向量，将图形$B$沿每个向量移动，所有的最终位置的并便是闵可夫斯基和（具有交换律）

vector<Point> MinkowskiSum(const vector<Point>& a, const vector<Point>& b) {
	int n = a.size(), m = b.size();
	Point sa = a[0], sb = b[0];
	for (int i = 0; i < n; i++) {
		if (a[i].y < sa.y || (a[i].y == sa.y && a[i].x < sa.x))
			sa = a[i];
	}
	for (int i = 0; i < m; i++) {
		if (b[i].y < sb.y || (b[i].y == sb.y && b[i].x < sb.x)) {
			sb = b[i];
		}
	}
	auto s = sa + sb;
	vector<Point> d(n + m);
	for (int i = 0; i < n; i++)
		d[i] = a[(i + 1) % n] - a[i];
	for (int i = 0; i < m; i++)
		d[n + i] = b[(i + 1) % m] - b[i];
	sort(d.begin(), d.end(), [&](const Point& A, const Point& B) {
		if (A.up() ^ B.up())
			return A.up() > B.up();
		return (A ^ B) > 0;
	});
	vector<Point> c(n + m);
	c[0] = s;
	for (int i = 0; i < n + m - 1; i++)
		c[i + 1] = c[i] + d[i];
	return c;
}

三维计算几何

多面体欧拉定理

平面图欧拉定理：$V-E+F=k+1$（$V$:点数，$E$:边数，$F$:面数，$k$:连通分量）

顶点数  -  棱长数  +  表面数  =  $2$

点数为$n$，则边数最多为$3n-6$，面数最多为$2n-4$

点类

struct Point {
    double x, y, z;
    Point(double x = 0, double y = 0, double z = 0) : x(x), y(y), z(z) {}
    bool operator<(const Point& B) const {
        if (cmp(x, B.x))
            return x < B.x;
        if (cmp(y, B.y))
            return y < B.y;
        return z < B.z;
    }
    bool operator==(const Point& B) const {
        return !cmp(x, B.x) && !cmp(y, B.y) && !cmp(z, B.z);
    }
    Point operator+(const Point& B) const {
        return Point(x + B.x, y + B.y, z + B.z);
    }
    Point operator-(const Point& B) const {
        return Point(x - B.x, y - B.y, z - B.z);
    }
    Point operator*(const double a) const {
        return Point(x * a, y * a, z * a);
    }
    Point operator/(const double a) const {
        return Point(x / a, y / a, z / a);
    }
    double operator*(const Point& B) const {
        return x * B.x + y * B.y + z * B.z;
    }
    Point operator^(const Point& B) const {
        return {y * B.z - B.y * z, z * B.x - x * B.z, x * B.y - y * B.x};
    }
    double length() {
        return sqrt(x * x + y * y + z * z);
    }
    double rand_eps() {
        return ((double)rand() / RAND_MAX - 0.5) * eps;
    }
    void shake() {  // 微小扰动, 防止出现四点共面的情况
        x += rand_eps(), y += rand_eps(), z += rand_eps();
    }
};

直线类

struct Line {
	Point p, v;
	Line() {}
	Line(const Point& A, const Point& B) : p(A), v(B) {}
	Point foot_point(Point a) {     // 垂足
		return p + v * ((v * (a - p)) / (v * v));
	}
	double distance(const Point& a) {   // 点到直线的距离
		return fabs((v ^ (a - p)).length() / v.length());
	}
};

平面类

Point q[N];
struct Plane {
    int v[3];
    Point norm() {  // 平面法向量
        return (q[v[1]] - q[v[0]]) ^ (q[v[2]] - q[v[0]]);
    }
    double area() {  // 面积
        return norm().length() / 2;
    }
    bool above(Point a) {  // 点a是否在平面的上方
        return ((a - q[v[0]]) * norm()) >= 0;
    }
};

三维凸包

Point q[N];
Plane plane[N], np[N];  // 通过微小扰动, 防止四点共面
int n, m;
bool g[N][N];           // 某条边是否是分界线
void get_convex_3d() {  // O(n^2)
    plane[m++] = {0, 1, 2}, plane[m++] = {2, 1, 0};
    for (int i = 3; i < n; i++) {
        int cnt = 0;
        for (int j = 0; j < m; j++) {
            bool t = plane[j].above(q[i]);
            if (!t)
                np[cnt++] = plane[j];
            for (int k = 0; k < 3; k++)
                g[plane[j].v[k]][plane[j].v[(k + 1) % 3]] = t;
        }
        for (int j = 0; j < m; j++) {
            for (int k = 0; k < 3; k++) {
                int a = plane[j].v[k], b = plane[j].v[(k + 1) % 3];
                if (g[a][b] && !g[b][a])
                    np[cnt++] = {a, b, i};
            }
        }
        m = cnt;
        for (int j = 0; j < m; j++)
            plane[j] = np[j];
    }
}
signed main() {
    cin >> n;
    for (int i = 0; i < n; i++) {
        cin >> q[i].x >> q[i].y >> q[i].z;
        q[i].shake();
    }
    get_convex_3d();
    double ans = 0;
    for (int i = 0; i < m; i++)
        ans += plane[i].area();
    cout << fixed << setprecision(6) << ans << endl;
}

两平面的交线

Line planeIntersection(Point a, Point v1, Point b, Point v2) {	// 平面AB的交线 (已知平面上一点及其法向量)
	Point e = (v1 ^ v2), u = (v1 ^ e);
	double d = v2 * u;
	if (!sign(d))
		return Line(Point(0, 0), Point(0, 0));
	Point q = a + u * ((v2 * (b - a)) / d);
	return Line(q, e);
}

球

经纬度转化为空间坐标

double torad(double a) {   // 角度转化为弧度
	return a / 180 * pi;
}
Point getPoint(double R, double a, double b) {  // a经度,b维度,均为角度
	a = torad(a), b = torad(b);
	return Point(R * cos(a) * cos(b), R * cos(a) * sin(b), R * sin(a));
}

两点球面距离

double ball_dist(double R, Point a, Point b) {
	double d = (a - b).length();
	return 2 * asin(d / 2 / R) * R;
}

直线与球的交点

int lineInterBall(Line L, Point O, double R, Point& p1, Point& p2) {
	Point foot = L.foot_point(O);
	double distOF = (O - foot).length();
	if (cmp(distOF, R) > 0)
		return 0;
	else if (!cmp(distOF, R)) {
		p1 = p2 = foot;
		return 1;
	}
	Point delta = L.v.trunc(sqrt(R * R - distOF * distOF));
	p1 = foot + delta, p2 = foot - delta;
	return 2;
}