零. 三态函数
const double eps = 1e-8;
const double Pi = acos(-1.0);
int sgn(double x) {
if(fabs(x) < eps) return 0;
return x < 0 ? -1 : 1;
}
一. 点类及其常见操作
typedef struct Point {
double x, y;
Point(double x=0, double y=0): x(x), y(y) {}
void input() { scanf("%lf%lf", &x, &y); }
double Len() { return sqrt(x*x + y*y); } // 模长
double Len2() { return x*x + y*y; } // 模长平方
Point operator + (Point r) { return Point(x+r.x, y+r.y); } // 向量相加
Point operator - (Point r) { return Point(x-r.x, y-r.y); } // 向量相减
Point operator * (double r) { return Point(x*r, y*r); } // 向量模改
Point operator / (double r) { return Point(x/r, y/r); } // 向量模改
double operator ^ (Point r) { return x*r.y - r.x*y; } // 叉积
double operator * (Point r) { return x*r.x + y*r.y; } // 点积
bool operator == (Point r) { return sgn(x-r.x)==0 && sgn(y-r.y)==0; }
bool operator < (Point r) { return sgn(y-r.y)<0 || (sgn(y-r.y)==0&&sgn(x-r.x)<0); }
void transXY(double rad) { // 点绕原点旋转 rad 弧度
double tx = x, ty = y;
x = tx*cos(rad) - ty*sin(rad);
y = tx*sin(rad) + ty*cos(rad);
}
}Vector;
-
点积
double Dot(Point a, Point b) { return a.x*b.x + a.y*b.y; } -
叉积
double Cross(Point a, Point b) { return a.x*b.y - b.x*a.y; } double Cross(Point pt, Point a, Point b) { return Cross(a-pt, b-pt); } -
向量的旋转（逆时针为正方向）
Vector Ratate(Vector v, double rad) { return v.transXY(rad), v; } -
向量的垂直向量
Vector Normal(Vector) { return Vector(-v.y, v.x); }
二. 线类（两点式）
struct Line {
Point s, e;
Line(Point s=o, Point e=o): s(s), e(e) {}
void input() { s.input(), e.input(); }
pair<int, Point> operator & (Line r) {
Point res = s;
if(sgn(Cross(r.e - r.s, e - s)) == 0) {
if(sgn(Cross(r.e - r.s, e - r.s)) == 0) return make_pair(1, res);
return make_pair(2, res);
}
double t = Cross(s-r.s, r.e-r.s)/Cross(r.e-r.s, e-s);
res.x += (e.x-s.x) * t;
res.y += (e.y-s.y) * t;
return make_pair(0, res);
}
};
-
线段相交
bool Seg_inter_Seg(Line A, Line B) { return max(A.s.x, A.e.x) >= min(B.s.x, B.e.x) && max(A.s.y, A.e.y) >= min(B.s.y, B.e.y) && max(B.s.x, B.e.x) >= min(A.s.x, A.e.x) && max(B.s.y, B.e.y) >= min(A.s.y, A.e.y) && sgn(Cross(A.e-A.s, B.e-A.s)) * sgn(Cross(A.e-A.s, B.s-A.s)) <= 0 && sgn(Cross(B.e-B.s, A.e-B.s)) * sgn(Cross(B.e-B.s, A.s-B.s)) <= 0; } -
点是否在向量上
bool OnSeg(Point pt, Line L) { return sgn(Cross(L.s - pt, L.e - pt)) == 0 && sgn(Dot(L.s - pt, L.e - pt)) <= 0; } -
判断直线和线段是否相交
bool Seg_inter_Line(Line seg, Line L) { return sgn(Cross(seg.s-L.s, L.e-L.s)) * sgn(Cross(seg.e-L.s, L.e-L.s)) <= 0; } -
两点间的距离
double Dist(Point A, Point B) { return (A-B).Len(); } -
点到直线的距离
double Point_to_Line(Point pt, Line L) { return fabs(Cross(pt-L.s, L.s-L.e)/Dist(L.s, L.e)); } -
点到线段的距离
double Point_to_Seg(Point pt, Line L) { if(sgn(Dot(pt-L.s, L.e-L.s)) < 0) return Dist(pt, L.s); if(sgn(Dot(pt-L.e, L.s-L.e)) < 0) return Dist(pt, L.e); return Point_to_Line(pt, L); } -
两线段的距离
double Seg_to_Seg(Line A, Line B) { if (Seg_inter_Seg(A, B)) return 0; return min(Point_to_Seg(A.e, B), Point_to_Seg(A.s, B)); }
三. 线类（点向式）
struct Plane {
Point p; // 直线上的一点
Vector v; // 方向向量
double ang; // 从 x 轴转向 v 方向的夹角
Plane(Point p=o, Vector v=Vector(1,0)): p(p), v(v) { ang = atan2(v.y, v.x); }
void Angle() { ang = atan2(v.y, v.x); }
bool operator < (Plane r) { return sgn(ang - r.ang) > 0; }
};
-
判断点 p p p 在直线 L L L 的左边，即点 p p p 在 L L L 的外边。
bool OnLeft(Plane L, Point p) { return sgn(Cross(L.v, p - L.p)) > 0; } -
判断点是否在直线上
bool OnLine(Plane L, Point pt) { return sgn(Cross(pt - L.p, L.v)) == 0; } -
判断点是否不在直线的右边
bool NotOnRight(Plane L, Point pt) { return sgn(Cross(L.v, pt - L.p)) >= 0; } -
两直线的交点
Point Cross_point(Plane A, Plane B) { Vector u = A.p - B.p; double t = Cross(B.v, u)/Cross(A.v, B.v); return A.p + A.v * t; } -
计算到固定点的距离和
double Dist_and(Point pt, Point* p, int n) { double res = 0; for(int i = 0; i < n; i ++) res += Dist(pt, p[i]); return res; }
四. 圆类
struct Circle{
Point p; double r;
Circle(Point p=o, double r=0): p(p), r(r) {}
void input() { p.input(); scanf("%lf", &r); }
};
-
过圆外一点与圆的切线
Plane Tangent(Point p, Circle c, int clock) { Plane res; Vector v = c.p - p; double d = Dist(c.p, p); double rad = asin(c.r/d); if(clock > 0) res.v = Ratate(v, rad); else res.v = Ratate(v, -rad); res.p = p; return res; }
五. 凸包相关
-
计算凸包的周长
double Perimeter_polygon(Point* p, int n) { p[n] = p[0]; double ans = 0; for(int i = 0; i < n; i ++) ans += Dist(p[i], p[i+1]); return ans; } -
计算凸包的面积
double Area_polygon(Point* p, int n) { p[n] = p[0]; double ans = 0; for(int i = 0; i < n; i ++) ans += Cross(p[i], p[i+1]); return ans; } -
判断圆是否在凸多边形内部
bool Circle_in_polygon(Circle O, Point* p, int n) { p[n] = p[0]; for(int i = 0; i < n; i ++) { if(sgn(Point_to_Seg(O.p, Line(p[i], p[i+1])) - O.r) < 0) return false; } return true; } -
判断是否为凸包，方向不一定是逆时针
bool Judge_Convex_hull(Point* p, int n){ int dir = 0; p[n] = p[0], p[n+1] = p[1]; for(int i = 1; i <= n; i ++) { int u = sgn(Cross(p[i] - p[i-1], p[i+1] - p[i])); if(! dir) dir = u; if(u * dir < 0) return false; } return true; } -
将凸包变为逆时针顺序
void Anti_clockwise(Point* p, int n) { p[n] = p[0]; for(int i = 0; i < n - 1; i ++) { if(sgn(Cross(p[i], p[i+1], p[i+2])) > 0) return ; else if(sgn(Cross(p[i], p[i+1], p[i+2])) < 0) { for(int j = 0; j < n/2; j ++) { Point pt = p[j]; p[j] = p[n-j-1], p[n-j-1] = pt; } return ; } } } -
判断点是否在多边形（不限于凸包）内部
int in_polygon(Point pt, Point* p, int n) { for(int i = 0; i < n; i ++) if(pt == p[i]) return 3; for(int i = 0; i < n; i ++) { int j = i + 1 == n ? 0 : i + 1; if(OnSeg(pt, Line(p[i], p[j]))) return 2; } int num = 0; for(int i = 0; i < n; i ++) { int j = i + 1 == n ? 0 : i + 1; int c = sgn(Cross(pt-p[j], p[i]-p[j])); int u = sgn(p[j].y - pt.y); int v = sgn(p[i].y - pt.y); if(c > 0 && u >= 0 && v < 0) num ++; if(c < 0 && u < 0 && v >= 0) num --; } return num != 0; } -
求凸包（Andrew算法）
int Convex_hull(Point* p, Point* ch, int n) { sort(p, p + n); int v = 0; for(int i = 0; i < n; i ++) { while(v > 1 && sgn(Cross(ch[v-1]-ch[v-2], p[i]-ch[v-2])) <= 0) v --; ch[v ++] = p[i]; } int j = v; for(int i = n - 2; ~ i; i --) { while(v > j && sgn(Cross(ch[v-1]-ch[v-2], p[i]-ch[v-2])) <= 0) v --; ch[v ++] = p[i]; } if(n > 1) v --; return v; } -
两个多边形的面积交，两个凸包都是逆时针
double CPI_Area(Point* a, int na, Point* b, int nb) { static Point p[N], tmp[N]; int tn, sflag, eflag; a[na] = a[0], b[nb] = b[0]; memcpy(p, b, sizeof(Point)*(nb+1)); /// 整个过程类似半平面交 /// 枚举一个凸包的所有边对另一个凸包进行切割 for(int i = 0; i < na && nb > 2; i ++) { sflag = sgn(Cross(a[i], a[i+1], p[0])); for(int j = tn = 0; j < nb; j ++, sflag = eflag) { if(sflag >= 0) tmp[tn ++] = p[j]; /// 当前顶点不在右侧, 直接保留 eflag = sgn(Cross(a[i], a[i+1], p[j+1])); /// 异号说明两直线相交, 计算交点切割边 if((sflag ^ eflag) == -2) tmp[tn ++] = (Line(a[i], a[i+1]) & Line(p[j], p[j+1])).second; } memcpy(p, tmp, sizeof(Point)*tn); nb = tn, p[nb] = p[0]; } if(nb < 3) return 0; /// 返回切割后的面积就是面积交 return Area_polygon(p, nb); } -
两个多边形的面积并
double SPI_Area(Point* a, int na, Point* b, int nb) { Point t1[4], t2[4]; double res = 0;int num1, num2; a[na] = t1[0] = a[0], b[nb] = t2[0] = b[0]; /// 以各自的 0 处的点作为基本点, 并保证整个三角形都是都是逆时针 for(int i = 2; i < na; i ++) { t1[1] = a[i-1], t1[2] = a[i]; num1 = sgn(Cross(t1[0], t1[1], t1[2])); if(num1 < 0) swap(t1[1], t1[2]); for(int j = 2; j < nb; j ++) { t2[1] = b[j-1], t2[2] = b[j]; num2 = sgn(Cross(t2[0], t2[1], t2[2])); if(num2 < 0) swap(t2[1], t2[2]); res += CPI_Area(t1, 3, t2, 3) * num1 * num2; } } return fabs(res); }
六. 半平面交
-
算法实现
int HPI(Plane* L, Point *ch, int n) { sort(L, L + n); /// 极角排序 int h = 0, t = 0, v = 0; static Plane q[N]; /// 双端队列 static Point p[N]; /// 两个相邻半平面的交点 q[h] = L[0]; for(int i = 1; i < n; i ++) { /// 删除队尾的半平面 while(h < t && ! OnLeft(L[i], p[t-1])) t --; /// 删除队首的半平面 while(h < t && ! OnLeft(L[i], p[h])) h ++; q[++ t] = L[i]; /// 将当前的半平面加入双端队列的尾部 /// 极角相同的两个半平面保留左边 if(sgn(Cross(q[t].v, q[t-1].v)) == 0) { t --; if(OnLeft(q[t], L[i].p)) q[t] = L[i]; } /// 计算队列尾部的半平面交点 if(h < t) p[t - 1] = Cross_point(q[t-1], q[t]); } /// 删除队列尾部无用的半平面 while(h < t && ! OnLeft(q[h], p[t-1])) t --; if(t - h <= 1) return 0; /// 空集 p[t] = Cross_point(q[t], q[h]); /// 计算队列首尾部的交点 for(int i = h; i <= t; i ++) ch[v ++] = p[i]; /// 复制 return v; /// 返回凸多边形大小 } -
判断一个点是否有核（单独一个点也算）
bool HPI_Core(Plane* L, int n){ sort(L, L + n); int h = 0, t = 0; static Plane q[N]; static Point p[N]; q[0] = L[0]; for(int i = 1; i < n; i ++) { while(h < t && ! OnLeft(L[i], p[t-1])) t --; while(h < t && ! OnLeft(L[i], p[h])) h ++; q[++ t] = L[i]; if(sgn(Cross(q[t].v, q[t-1].v)) == 0) { t --; if(OnLeft(q[t], L[i].p)) q[t] = L[i]; } if(h < t) p[t-1] = Cross_point(q[t], q[t-1]); } while(h < t && ! OnLeft(q[h], p[t-1])) t --; return t - h + 1 > 2; }
七. 旋转卡壳
-
旋转卡壳求凸包直径（宽度）
double Rotating_cailper(Point* p, int n) { double res = 0; if(n == 2) return Dist(p[0], p[1]); p[n] = p[0]; for(int i = 0, j = 2; i < n; i ++) { while(sgn(Cross(p[i], p[i+1], p[j]) - Cross(p[i], p[i+1], p[j+1])) < 0) if(++ j >= n) j -= n; /// 计算凸包宽度时取对踵点到直线的距离 /// res = max(res, Point_to_Line(p[j], Line(p[i], p[i+1]))); res = max(res, max(Dist(p[j], p[i]), Dist(p[j], p[i+1]))); } return res; } -
旋转卡壳求两个凸包的最近距离
double Nearest_dist(Point* p, int n, Point* q, int m) { p[n] = p[0], q[m] = q[0]; double tmp, res = 1e18; int u = 0, v = 0; for(int i = 0; i < n; i ++) if(p[i].y < p[u].y) u = i; for(int i = 0; i < m; i ++) if(q[i].y > q[v].y) v = i; for(int i = 0; i < n; i ++) { while(sgn(tmp = Cross(p[u], p[u+1], q[v+1]) - Cross(p[u], p[u+1], q[v])) > 0) if(++ v >= m) v -= m; if(sgn(tmp) < 0) res = min(res, Point_to_Seg(q[v], Line(p[u], p[u+1]))); else res = min(res, Seg_to_Seg(Line(p[u], p[u+1]), Line(q[v], q[v+1]))); if(++ u >= n) u -= n; } return res; }