二维计算几何模板整理
/***********三分法求函数极值*************/void solve(){ double L, R, m, mm, mv, mmv; while (L + eps < R) { m = (L + R) / 2; mm = (m + R) / 2; mv = calc(m); mmv = calc(mm); if (mv <= mmv) R = mm; //三分法求最大值时改为mv>=mmv else L = m; }}/*************基础***********/int dcmp(double x) { if(fabs(x) < eps) return 0; else return x < 0 ? -1 : 1;}struct Point { double x, y; Point(double x=0, double y=0):x(x),y(y) { }};typedef Point Vector;Vector operator + (Vector A, Vector B) { return Vector(A.x+B.x, A.y+B.y); }Vector operator - (Point A, Point B) { return Vector(A.x-B.x, A.y-B.y); }Vector operator * (Vector A, double p) { return Vector(A.x*p, A.y*p); }Vector operator / (Vector A, double p) { return Vector(A.x/p, A.y/p); }bool operator < (const Point& a, const Point& b) { return a.x < b.x || (a.x == b.x && a.y < b.y);}bool operator == (const Point& a, const Point &b) { return dcmp(a.x-b.x) == 0 && dcmp(a.y-b.y) == 0;}double Dot(Vector A, Vector B) { return A.x*B.x + A.y*B.y; }double Length(Vector A) { return sqrt(Dot(A, A)); }double Angle(Vector A, Vector B) { return acos(Dot(A, B) / Length(A) / Length(B)); }double angle(Vector v) { return atan2(v.y, v.x); }double Cross(Vector A, Vector B) { return A.x*B.y - A.y*B.x; }Vector vecunit(Vector x){ return x / Length(x);} //单位向量 Vector Normal(Vector x) { return Point(-x.y, x.x) / Length(x);} //垂直法向量 Vector Rotate(Vector A, double rad) { return Vector(A.x*cos(rad)-A.y*sin(rad), A.x*sin(rad)+A.y*cos(rad));}double Area2(const Point A, const Point B, const Point C) { return Length(Cross(B-A, C-A)); }/****************直线与线段**************///求直线p+tv和q+tw的交点 Cross(v, w) == 0无交点 Point GetLineIntersection(Point p, Vector v, Point q, Vector w) { Vector u = p-q; double t = Cross(w, u) / Cross(v, w); return p + v*t; } //点p在直线ab的投影Point GetLineProjection(Point P, Point A, Point B) { Vector v = B-A; return A+v*(Dot(v, P-A) / Dot(v, v));}//点到直线距离double DistanceToLine(Point P, Point A, Point B) { Vector v1 = B - A, v2 = P - A; return fabs(Cross(v1, v2)) / Length(v1); // 如果不取绝对值，得到的是有向距离}//点在p线段上bool OnSegment(Point p, Point a1, Point a2) { return dcmp(Cross(a1-p, a2-p)) == 0 && dcmp(Dot(a1-p, a2-p)) < 0;}// 过两点p1, p2的直线一般方程ax+by+c=0// (x2-x1)(y-y1) = (y2-y1)(x-x1)void getLineGeneralEquation(const Point& p1, const Point& p2, double& a, double& b, double &c) { a = p2.y-p1.y; b = p1.x-p2.x; c = -a*p1.x - b*p1.y;}//点到线段距离double DistanceToSegment(Point p, Point a, Point b) { if(a == b) return Length(p-a); Vector v1 = b-a, v2 = p-a, v3 = p-b; if(dcmp(Dot(v1, v2)) < 0) return Length(v2); else if(dcmp(Dot(v1, v3)) > 0) return Length(v3); else return fabs(Cross(v1, v2)) / Length(v1); } //两线段最近距离double dis_pair_seg(Point p1, Point p2, Point p3, Point p4) { return min(min(DistanceToSegment(p1, p3, p4), DistanceToSegment(p2, p3, p4)), min(DistanceToSegment(p3, p1, p2), DistanceToSegment(p4, p1, p2))); }//线段相交判定 bool SegmentItersection(Point a1, Point a2, Point b1, Point b2) { double c1 = Cross(a2-a1, b1-a1), c2 = Cross(a2-a1, b2-a1), c3 = Cross(b2-b1, a1-b1), c4 = Cross(b2-b1, a2-b1); return dcmp(c1)*dcmp(c2) < 0 && dcmp(c3)*dcmp(c4) < 0; } // 有向直线。它的左边就是对应的半平面struct Line { Point P; // 直线上任意一点 Vector v; // 方向向量 double ang; // 极角，即从x正半轴旋转到向量v所需要的角（弧度） Line() {} Line(Point P, Vector v):P(P),v(v){ ang = atan2(v.y, v.x); } bool operator < (const Line& L) const { return ang < L.ang; }};//两直线交点Point GetLineIntersection(Line a, Line b) { return GetLineIntersection(a.p, a.v, b.p, b.v);}// 点p在有向直线L的左边（线上不算）bool OnLeft(const Line& L, const Point& p) { return Cross(L.v, p-L.P) > 0;}// 二直线交点，假定交点惟一存在Point GetLineIntersection(const Line& a, const Line& b) { Vector u = a.P-b.P; double t = Cross(b.v, u) / Cross(a.v, b.v); return a.P+a.v*t;}// 半平面交主过程vector<Point> HalfplaneIntersection(vector<Line> L) { int n = L.size(); sort(L.begin(), L.end()); // 按极角排序 int first, last; // 双端队列的第一个元素和最后一个元素的下标 vector<Point> p(n); // p[i]为q[i]和q[i+1]的交点 vector<Line> q(n); // 双端队列 vector<Point> ans; // 结果 q[first=last=0] = L[0]; // 双端队列初始化为只有一个半平面L[0] for(int i = 1; i < n; i++) { while(first < last && !OnLeft(L[i], p[last-1])) last--; while(first < last && !OnLeft(L[i], p[first])) first++; q[++last] = L[i]; if(fabs(Cross(q[last].v, q[last-1].v)) < eps) { // 两向量平行且同向，取内侧的一个 last--; if(OnLeft(q[last], L[i].P)) q[last] = L[i]; } if(first < last) p[last-1] = GetLineIntersection(q[last-1], q[last]); } while(first < last && !OnLeft(q[first], p[last-1])) last--; // 删除无用平面 if(last - first <= 1) return ans; // 空集 p[last] = GetLineIntersection(q[last], q[first]); // 计算首尾两个半平面的交点 // 从deque复制到输出中 for(int i = first; i <= last; i++) ans.push_back(p[i]); return ans;}/***********多边形**************///多边形有向面积double PolygonArea(vector<Point> p) { int n = p.size(); double area = 0; for(int i = 1; i < n-1; i++) area += Cross(p[i]-p[0], p[i+1]-p[0]); return area/2;}// 点集凸包// 如果不希望在凸包的边上有输入点，把两个 <= 改成 <// 注意：输入点集会被修改vector<Point> ConvexHull(vector<Point>& p) { // 预处理，删除重复点 sort(p.begin(), p.end()); p.erase(unique(p.begin(), p.end()), p.end()); int n = p.size(); int m = 0; vector<Point> ch(n+1); for(int i = 0; i < n; i++) { while(m > 1 && Cross(ch[m-1]-ch[m-2], p[i]-ch[m-2]) <= 0) m--; ch[m++] = p[i]; } int k = m; for(int i = n-2; i >= 0; i--) { while(m > k && Cross(ch[m-1]-ch[m-2], p[i]-ch[m-2]) <= 0) m--; ch[m++] = p[i]; } if(n > 1) m--; ch.resize(m); return ch;}// 凸包直径返回 点集直径的平方int diameter2(vector<Point>& points) { vector<Point> p = ConvexHull(points); int n = p.size(); if(n == 1) return 0; if(n == 2) return Dist2(p[0], p[1]); p.push_back(p[0]); // 免得取模 int ans = 0; for(int u = 0, v = 1; u < n; u++) { // 一条直线贴住边p[u]-p[u+1] for(;;) { int diff = Cross(p[u+1]-p[u], p[v+1]-p[v]); if(diff <= 0) { ans = max(ans, Dist2(p[u], p[v])); // u和v是对踵点 if(diff == 0) ans = max(ans, Dist2(p[u], p[v+1])); // diff == 0时u和v+1也是对踵点 break; } v = (v + 1) % n; } } return ans;}//两凸包最近距离double RC_Distance(Point *ch1, Point *ch2, int n, int m) { int q=0, p=0; REP(i, n) if(ch1[i].y-ch1[p].y < -eps) p=i; REP(i, m) if(ch2[i].y-ch2[q].y > eps) q=i; ch1[n]=ch1[0]; ch2[m]=ch2[0]; double tmp, ans=1e100; REP(i, n) { while((tmp = Cross(ch1[p+1]-ch1[p], ch2[q+1]-ch1[p]) - Cross(ch1[p+1]-ch1[p], ch2[q]- ch1[p])) > eps) q=(q+1)%m; if(tmp < -eps) ans = min(ans,DistanceToSegment(ch2[q],ch1[p],ch1[p+1])); else ans = min(ans,dis_pair_seg(ch1[p],ch1[p+1],ch2[q],ch2[q+1])); p=(p+1)%n; } return ans; } //凸包最大内接三角形double RC_Triangle(Point* res,int n)// 凸包最大内接三角形 { if(n<3) return 0; double ans=0, tmp; res[n] = res[0]; int j, k; REP(i, n) { j = (i+1)%n; k = (j+1)%n; while((j != k) && (k != i)) { while(Cross(res[j] - res[i], res[k+1] - res[i]) > Cross(res[j] - res[i], res[k] - res[i])) k= (k+1)%n; tmp = Cross(res[j] - res[i], res[k] - res[i]);if(tmp > ans) ans = tmp; j = (j+1)%n; } } return ans; } //模拟退火求费马点 保存在ptres中double fermat_point(Point *pt, int n, Point& ptres) { Point u, v; double step = 0.0, curlen, explen, minlen; int i, j, k, idx; bool flag; u.x = u.y = v.x = v.y = 0.0; REP(i, n) { step += fabs(pt[i].x) + fabs(pt[i].y); u.x += pt[i].x; u.y += pt[i].y; } u.x /= n; u.y /= n; flag = 0; while(step > eps) { for(k = 0; k < 10; step /= 2, ++k) for(i = -1; i <= 1; ++i) for(j = -1; j <= 1; ++j) { v.x = u.x + step*i; v.y = u.y + step*j; curlen = explen = 0.0; REP(idx, n) { curlen += dist(u, pt[idx]); explen += dist(v, pt[idx]); } if(curlen > explen) { u = v; minlen = explen; flag = 1; } } } ptres = u; return flag ? minlen : curlen; } //最近点对bool cmpxy(const Point& a, const Point& b){ if(a.x != b.x) return a.x < b.x; return a.y < b.y;}bool cmpy(const int& a, const int& b){ return point[a].y < point[b].y;}double Closest_Pair(int left, int right){ double d = INF; if(left==right) return d; if(left + 1 == right) return dis(left, right); int mid = (left+right)>>1; double d1 = Closest_Pair(left,mid); double d2 = Closest_Pair(mid+1,right); d = min(d1,d2); int i,j,k=0; //分离出宽度为d的区间 for(i = left; i <= right; i++) { if(fabs(point[mid].x-point[i].x) <= d) tmpt[k++] = i; } sort(tmpt,tmpt+k,cmpy); //线性扫描 for(i = 0; i < k; i++) { for(j = i+1; j < k && point[tmpt[j]].y-point[tmpt[i]].y<d; j++) { double d3 = dis(tmpt[i],tmpt[j]); if(d > d3) d = d3; } } return d;}/************圆************/struct Circle { Point c; double r; Circle(){} Circle(Point c, double r):c(c), r(r){} Point point(double a) //根据圆心角求点坐标 { return Point(c.x+cos(a)*r, c.y+sin(a)*r); } };//直线与圆交点 返回个数int getLineCircleIntersection(Line L, Circle C, double& t1, double& t2, vector<Point>& sol){ double a = L.v.x, b = L.p.x - C.c.x, c = L.v.y, d = L.p.y - C.c.y; double e = a*a + c*c, f = 2*(a*b + c*d), g = b*b + d*d - C.r*C.r; double delta = f*f - 4*e*g; // 判别式 if(dcmp(delta) < 0) return 0; // 相离 if(dcmp(delta) == 0) { // 相切 t1 = t2 = -f / (2 * e); sol.push_back(L.point(t1)); return 1; } // 相交 t1 = (-f - sqrt(delta)) / (2 * e); sol.push_back(L.point(t1)); t2 = (-f + sqrt(delta)) / (2 * e); sol.push_back(L.point(t2)); return 2;}//两圆交点 返回个数int getCircleCircleIntersection(Circle C1, Circle C2, vector<Point>& sol) { double d = Length(C1.c - C2.c); if(dcmp(d) == 0) { if(dcmp(C1.r - C2.r) == 0) return -1; // 重合，无穷多交点 return 0; } if(dcmp(C1.r + C2.r - d) < 0) return 0; if(dcmp(fabs(C1.r-C2.r) - d) > 0) return 0; double a = angle(C2.c - C1.c); double da = acos((C1.r*C1.r + d*d - C2.r*C2.r) / (2*C1.r*d)); Point p1 = C1.point(a-da), p2 = C1.point(a+da); sol.push_back(p1); if(p1 == p2) return 1; sol.push_back(p2); return 2;}//三角形外接圆Circle CircumscribedCircle(Point p1, Point p2, Point p3) { double Bx = p2.x-p1.x, By = p2.y-p1.y; double Cx = p3.x-p1.x, Cy = p3.y-p1.y; double D = 2*(Bx*Cy-By*Cx); double cx = (Cy*(Bx*Bx+By*By) - By*(Cx*Cx+Cy*Cy))/D + p1.x; double cy = (Bx*(Cx*Cx+Cy*Cy) - Cx*(Bx*Bx+By*By))/D + p1.y; Point p = Point(cx, cy); return Circle(p, Length(p1-p));}//三角形内切圆Circle InscribedCircle(Point p1, Point p2, Point p3) { double a = Length(p2-p3); double b = Length(p3-p1); double c = Length(p1-p2); Point p = (p1*a+p2*b+p3*c)/(a+b+c); return Circle(p, DistanceToLine(p, p1, p2));}// 过点p到圆C的切线。v[i]是第i条切线的向量。返回切线条数int getTangents(Point p, Circle C, Vector* v) { Vector u = C.c - p; double dist = Length(u); if(dist < C.r) return 0; else if(dcmp(dist - C.r) == 0) { // p在圆上，只有一条切线 v[0] = Rotate(u, PI/2); return 1; } else { double ang = asin(C.r / dist); v[0] = Rotate(u, -ang); v[1] = Rotate(u, +ang); return 2; }}//所有经过点p 半径为r 且与直线L相切的圆心vector<Point> CircleThroughPointTangentToLineGivenRadius(Point p, Line L, double r) { vector<Point> ans; double t1, t2; getLineCircleIntersection(L.move(-r), Circle(p, r), t1, t2, ans); getLineCircleIntersection(L.move(r), Circle(p, r), t1, t2, ans); return ans;}//半径为r 与a b两直线相切的圆心vector<Point> CircleTangentToLinesGivenRadius(Line a, Line b, double r) { vector<Point> ans; Line L1 = a.move(-r), L2 = a.move(r); Line L3 = b.move(-r), L4 = b.move(r); ans.push_back(GetLineIntersection(L1, L3)); ans.push_back(GetLineIntersection(L1, L4)); ans.push_back(GetLineIntersection(L2, L3)); ans.push_back(GetLineIntersection(L2, L4)); return ans;}//与两圆相切 半径为r的所有圆心vector<Point> CircleTangentToTwoDisjointCirclesWithRadius(Circle c1, Circle c2, double r) { vector<Point> ans; Vector v = c2.c - c1.c; double dist = Length(v); int d = dcmp(dist - c1.r -c2.r - r*2); if(d > 0) return ans; getCircleCircleIntersection(Circle(c1.c, c1.r+r), Circle(c2.c, c2.r+r), ans); return ans;}//多边形与圆相交面积Point GetIntersection(Line a, Line b) //线段交点 { Vector u = a.p-b.p; double t = Cross(b.v, u) / Cross(a.v, b.v); return a.p + a.v*t; } bool OnSegment(Point p, Point a1, Point a2) { return dcmp(Cross(a1-p, a2-p)) == 0 && dcmp(Dot(a1-p, a2-p)) < 0; } bool InCircle(Point x, Circle c) { return dcmp(sqr(c.r) - sqr(Length(c.c - x))) >= 0;} bool OnCircle(Point x, Circle c) { return dcmp(sqr(c.r) - sqr(Length(c.c - x))) == 0;}//线段与圆的交点 int getSegCircleIntersection(Line L, Circle C, Point* sol) { Vector nor = normal(L.v); Line pl = Line(C.c, nor); Point ip = GetIntersection(pl, L); double dis = Length(ip - C.c); if (dcmp(dis - C.r) > 0) return 0; Point dxy = vecunit(L.v) * sqrt(sqr(C.r) - sqr(dis)); int ret = 0; sol[ret] = ip + dxy; if (OnSegment(sol[ret], L.p, L.point(1))) ret++; sol[ret] = ip - dxy; if (OnSegment(sol[ret], L.p, L.point(1))) ret++; return ret; }double SegCircleArea(Circle C, Point a, Point b) //线段切割圆 { double a1 = angle(a - C.c); double a2 = angle(b - C.c); double da = fabs(a1 - a2); if (da > PI) da = PI * 2.0 - da; return dcmp(Cross(b - C.c, a - C.c)) * da * sqr(C.r) / 2.0; } double PolyCiclrArea(Circle C, Point *p, int n)//多边形与圆相交面积 { double ret = 0.0; Point sol[2]; p[n] = p[0]; REP(i, n) { double t1, t2; int cnt = getSegCircleIntersection(Line(p[i], p[i+1]-p[i]), C, sol); if (cnt == 0) { if (!InCircle(p[i], C) || !InCircle(p[i+1], C)) ret += SegCircleArea(C, p[i], p[i+1]); else ret += Cross(p[i+1] - C.c, p[i] - C.c) / 2.0; } if (cnt == 1) { if (InCircle(p[i], C) && !InCircle(p[i+1], C)) ret += Cross(sol[0] - C.c, p[i] - C.c) / 2.0, ret += SegCircleArea(C, sol[0], p[i+1]); else ret += SegCircleArea(C, p[i], sol[0]), ret += Cross(p[i+1] - C.c, sol[0] - C.c) / 2.0; } if (cnt == 2) { if ((p[i] < p[i + 1]) ^ (sol[0] < sol[1])) swap(sol[0], sol[1]); ret += SegCircleArea(C, p[i], sol[0]); ret += Cross(sol[1] - C.c, sol[0] - C.c) / 2.0; ret += SegCircleArea(C, sol[1], p[i+1]); } } return fabs(ret); }