计算几何

向量定义及运算

typedef double DB;

const DB eps = 1e-12;

int sgn(DB x) { return fabs(x) < eps ? 0 : (x > 0 ? 1 : -1); } // 精度判断

struct Point { // 点的定义，也可看作向量
	DB x, y;
	Point(DB x = 0, DB y = 0) : x(x), y(y) {}
	Point operator + (Point a) { return Point(x + a.x, y + a.y); } // 向量加法
	Point operator - (Point a) { return Point(x - a.x, y - a.y); } // 向量减法
	Point operator - () { return Point(-x, -y); } // 相反向量
	friend Point operator * (DB k, Point a) { return Point(k * a.x, k * a.y); } // 数量乘
	Point operator / (DB k) { return Point(x / k, y / k); } // 数量除
	DB operator % (Point a) { return x * a.x + y * a.y; } // 点乘，为了简化使用%，下同
	DB operator / (Point a) { return x * a.y - y * a.x; } // 叉乘的模长
	operator DB() { return sqrt(x*x + y*y); } // 模长，可使用类型转换
	bool operator < (const Point &a) const { int d1 = sgn(x-a.x), d2 = sgn(y-a.y); return d1 < 0 || (d1 == 0 && d2 < 0); } // 用于排序
};

点到直线投影点

(Prob)：给定三个点(P,P_1,P_1)，求(P)在(P_1P_2)上的投影点。

(Sol)：设投影点为(vec P_1+tvec v)，(vec v)由(P_1)指向(P_2)，则由垂直得到

[(vec P_1+tvec v-P)cdotvec v=0 ]

解得

[t=dfrac{(vec P-vec P_1)cdotvec v}{vec v^2} ]

代入即可得到投影点。

// Projection
Point Proj(Point P, Point P1, Point P2) {
	Point v = P2-P1; DB t = (P-P1)%v / (v%v);
	return P1 + t*v;
}

关于直线对称点

(Prob)：给定三个点(P,P_1,P_2)，求出(P)在直线(P_1P_2)上的对称点(P')。

(Sol)：求出(P)在(P_1P_2)上的投影点，然后中点公式解出来点(P')。

// Reflection
Point Ref(Point P, Point P1, Point P2) {
	return 2*Proj(P, P1, P2) - P;
}

判断向量方向

(Prob)：给定向量(overrightarrow{PP_1})和(overrightarrow{PP_2})，判断(overrightarrow{PP_2})相对于(overrightarrow{PP_1})的方向

// Direct1判断左右，Direct2判断前后
// 1:<90度 0:=90度（垂直） -1:>90度
int Direct1(Point P, Point P1, Point P2) { return sgn((P1-P)*(P2-P)); }
// 1:逆时针 0:共线 -1:顺时针
int Direct2(Point P, Point P1, Point P2) { return sgn((P1-P)/(P2-P)); }

这个用处很大，可以用来做很多判定。

判断线段是否相交

(Prob)：以点对的形式给定两条线段，判断其是否相交（端点相交合法）

(Sol)：通过“快速排斥实验”和“跨立实验”就能判断其是否相交。快速排斥实验指以线段为对角线构成的矩形判断是否有交；如果不交，显然线段不相交，这里取平行于直角坐标系的矩形即可。否则通过跨立实验判断线段所在直线是否相交，也就是判断端点是否在直线的两端即可。

// IsSegmentIntersection
bool IsSegInter(Point P1, Point P2, Point Q1, Point Q2) {
	if (min(P1.x, P2.x) <= max(Q1.x, Q2.x) && max(P1.x, P2.x) >= min(Q1.x, Q2.x) && min(P1.y, P2.y) <= max(Q1.y, Q2.y) && max(P1.y, P2.y) >= min(Q1.y, Q2.y)) // 快速排斥实验
		if (Direct2(P1, P2, Q1) * Direct2(P1, P2, Q2) <= 0 && Direct2(Q1, Q2, P1) * Direct2(Q1, Q2, P2) <= 0) return 1; // 跨立实验
	return 0;
}

直线求交

(Prob)：以点对的形式给定两条直线，求交点。

(Sol)：用(vec P+tvec v)表示直线1，(vec Q+tvec w)表示直线2，设交点在(vec P+t_0vec v)，显然有

[(vec P+t_0vec v-Q) imesvec w=0 ]

解得

[t_0=frac{(vec Q-vec P) imesvec w}{vec v imesvec w} ]

代入即可

// GetLineIntersection
Point GetLineInter(Point P1, Point P2, Point Q1, Point Q2) {
	Point v = P2-P1, w = Q2-Q1; DB t = (Q1-P1)/w / (v/w);
	return P1 + t*v;
}

点到直线距离

(Sol)：显然可以先求出投影点，然后直接求。也可以用点乘的性质解出来，这里选择前者。

// PointLineDistance
DB PLDist(Point P, Point P1, Point P2) {
	return DB(Proj(P, P1, P2)-P);
}

点到线段距离

(Sol)：如果投影点在线段上，就是点到直线距离；否则为点到线段两段距离的最小值。

// PointSegmentDistance
DB PSDist(Point P, Point P1, Point P2) {
	if (Direct1(P1, P2, P) * Direct1(P2, P1, P) >= 0) return PLDist(P, P1, P2);
	// 存在钝角说明投影点不在线段上
	return min(DB(P-P1), DB(P-P2));
}

线段间最小距离

(Prob)：给定两条线段(l_1,l_2)，(Pin l_1)且(Qin l_2)，求出(|PQ|)的最小值

(Sol)：如果相交，显然最小距离为0；否则(PQ)一定有一端在(l_1)或(l_2)的端点，求所有端点到线段距离最小值。

// SegmentDistance
DB SegDist(Point P1, Point P2, Point Q1, Point Q2) {
	if (IsSegInter(P1, P2, Q1, Q2)) return 0;
	return min(min(PSDist(P1, Q1, Q2), PSDist(P2, Q1, Q2)), min(PSDist(Q1, P1, P2), PSDist(Q2, P1, P2)));
}

多边形结构体

// 博主个性写法
struct Poly {
	std::vector<Point> a;
	Poly(int n = 0) { a.resize(n); rep0(i, n) a[i] = gp(); } // 读取n边形
	int size() { return a.size(); } // 获得n
	void add(Point P) { a.push_back(P); } // 手动加点
	Point &operator [] (int i) { int n = size(); return a[i >= n ? i-n : i]; } // 返回第i个点
};

多边形面积

(Sol)：选择任意起点(O)，将所有边(A_iA_{i+1})运用叉乘算出来(OA_i imes OA_{i+1})，取绝对值即能得到多边形面积的二倍。一般将(O)设在多边形的一个顶点上。

DB Area(Poly A) {
	DB res = 0;
	rep(i, 2, A.size()-1) res += (A[i]-A[0])/(A[i-1]-A[0]);
	return fabs(res)/2;
}

判断多边形是否为凸包

(Sol)：当且仅当沿着多边形走一圈，叉乘一直非负（或非正），说明旋转方向不变，是凸包。

bool IsConvex(Poly A) {
	bool f1 = 1, f2 = 1;
	rep0(i, A.size()) {
		int t = Direct2(A[i], A[i+1], A[i+2]);
		t && (t == 1 ? f1 = 0 : f2 = 0);
	}
	return f1 || f2;
}

判断点是否在多边形中

(Sol)：随机引入一条射线，如果与多边形相交次数为奇数，则在多边形内，偶数说明在多边形外。同时要判断在不在多边形上。注意采用随机引入一条射线能很大概率上避免交与顶点出而带来的麻烦。

// 1:在多边形中 0:不在多边形中 -1:在多边形上
int IsInPoly(Poly A, Point P) {
	DB angle = 2.0*PI*rand()/RAND_MAX; Point Q = P+1e5*Point(cos(angle), sin(angle));
	int cnt = 0;
	rep0(i, A.size()) {
		if (!Direct2(P, A[i], A[i+1]) && Direct1(A[i], A[i+1], P) * Direct1(A[i+1], A[i], P) >= 0) return -1;
		cnt ^= IsSegInter(A[i], A[i+1], P, Q);
	}
	return cnt;
}

求点集凸包

(Sol)：栈维护上下凸包，运用叉乘判断是否方向一致。最后拼起来即可。

Poly GetConvex(int n, Point P[]) { // 这里P已经按照x或者y轴排过序；凸包包含平角
	std::vector<int> s1, s2;
	rep(i, 1, n) {
		while (s1.size() > 1 && Direct2(P[s1[s1.size()-2]], P[s1[s1.size()-1]], P[i]) < 0) s1.pop_back();
		while (s2.size() > 1 && Direct2(P[s2[s2.size()-2]], P[s2[s2.size()-1]], P[i]) > 0) s2.pop_back();
		s1.push_back(i), s2.push_back(i);
	}
	Poly A;
	rep(i, 0, s1.size()-2) A.add(P[s1[i]]);
	per(i, s2.size()-1, 1) A.add(P[s2[i]]);
	return A;
}

平面最远点对（旋转卡壳）

(Sol)：求出所有点的凸包，跑旋转卡壳。旋转卡壳本质是相对于点(P_i)的点(P_j)距离上凸，同时转移(i)后最大值点往后走，可以采用双指针维护。

// 旋转卡壳，这里保证A是凸包
DB MaxConvexDist(Poly A) {
	DB ans = 0;
	for (int l = 0, r = 0; l < A.size(); l++) {
		while (r < A.size()-1 && DB(A[r]-A[l]) <= DB(A[r+1]-A[l])) r++;
		chkmax(ans, DB(A[r]-A[l]));
	}
	return ans;
}

半平面交

（咕咕咕）

三角形内、外接圆

由初中知识易知内切圆圆心为三条角平分线交点，外接圆圆心为三条中垂线交点。这里只需要两条直线即可。内切圆半径为圆心到任意一边的距离，外接圆半径为圆心到任意一顶点的距离。

// 内切圆
Circle GetInnerCir(Point A, Point B, Point C) {
	Point v1 = B-A, v2 = C-A, v3 = C-B;
	v1 = v1/DB(v1), v2 = v2/DB(v2), v3 = v3/DB(v3);
	Point P = GetLineInter(A, A + (v1+v2), B, B + (v3-v1));
	return (Circle){P, DB(P-Proj(P, A, B))};
}
// 外接圆
Circle GetOuterCir(Point A, Point B, Point C) {
	Point M = 0.5*(A+B), N = 0.5*(A+C), v = B-A, w = C-A;
	Point P = GetLineInter(M, M + Point(v.y, -v.x), N, N + Point(w.y, -w.x));
	return (Circle){P, DB(P-A)};
}

圆和直线交点

(Prob)：给定一个圆(O)和直线(l)，求直线交点。

(Sol)：先求出弦心距，然后根据垂径定理解出来两个交点。

// GetCircleLineIntersection
// 这里保证有交点
std::pair<Point, Point> GetCirLineInter(Circle C, Point P1, Point P2) {
	Point P = Proj(C.O, P1, P2), v = P2-P1; v = v/DB(v);
	DB d = DB(C.O-P);
	if (!sgn(d-C.r)) return std::make_pair(P, P); // 判断是否相等，防止开根号出现误差导致错误
	DB t = sqrt(C.r*C.r - d*d);
	return std::make_pair(P+t*v, P-t*v);
}

圆和圆的交点

(Prob)：给定两个圆，求交点。

(Sol)：相当于给定了三角形三边长以及两个点的坐标，利用余弦定理定出来第三个点投影点到其中一个点的距离，然后运用勾股定理解出来高，即可定位第三个点。

// GetCircleIntersection
// 保证有交点
std::pair<Point, Point> GetCirInter(Circle C1, Circle C2) {
	Point A = C1.O, B = C2.O, v = B-A; DB a = DB(v), x = (a*a+C1.r*C1.r-C2.r*C2.r)/(2*a); v = v/a;
	Point P = A+x*v;
	if (!sgn(C1.r-x)) return std::make_pair(P, P);
	DB h = sqrt(C1.r*C1.r - x*x); Point w = Point(v.y, -v.x);
	return std::make_pair(P+h*w, P-h*w);
}

反演变换

给定一个点(O)和常数(k)，点(P)的变换对应点就是在以(O)开始的射线(OP)上的一点(P')使得(|OP||OP'|=k)。

反演的重要性质：点<->点；过O的直线<->过O的直线；不过O的直线<->过O的圆；不过O的圆<->不过O的圆。

以上性质均可通过构造相似三角形证明。