二维平面计算几何模板

  1 //author Eterna
  2 #include<iostream>
  3 #include<algorithm>
  4 #include<cstdio>
  5 #include<vector>
  6 #include<cstring>
  7 #include<string>
  8 #include<cmath>
  9 #include<cstdlib>
 10 #include<utility>
 11 #include<deque>
 12 #include<queue>
 13 using namespace std;
 14 //------------------------------点与直线部分-----------------------------
 15 const double pi = acos(-1.0);
 16 const double eps = 1e-10;
 17 //double比较和0的大小关系
 18 inline int Dcmp(double x) {
 19     if (fabs(x) < eps)return 0;
 20     else return x < 0 ? -1 : 1;
 21 }
 22 struct point {
 23     friend istream& operator >>(istream& in, point& rhs) {
 24         in >> rhs.x >> rhs.y;
 25         return in;
 26     }
 27     friend ostream& operator <<(ostream& os, const point& rhs) {
 28         os << rhs.x << ‘ ‘ << rhs.y;
 29         return os;
 30     }
 31     bool operator == (const point& rhs)const {
 32         return Dcmp(x - rhs.x) == 0 && Dcmp(y - rhs.y) == 0;
 33     }
 34     bool operator < (const point& rhs)const {
 35         return x < rhs.x || (x == rhs.x && y < rhs.y);
 36     }
 37     bool operator >(const point& rhs)const {
 38         return !((*this) < rhs) && !((*this) == rhs);
 39     }
 40     point operator + (const point& p)const {
 41         return point(x + p.x, y + p.y);
 42     }
 43     point operator - (const point& p)const {
 44         return point(x - p.x, y - p.y);
 45     }
 46     point operator * (double p)const {
 47         return point(x * p, y * p);
 48     }
 49     point operator / (double p)const {
 50         return point(x / p, y / p);
 51     }
 52     //点乘
 53     double dot(const point& p) const {
 54         return x * p.x + y * p.y;
 55     }
 56     //叉乘
 57     double det(const point& p)const {
 58         return x * p.y - y * p.x;
 59     }
 60     point() {}
 61     point(double _x, double _y) :x(_x), y(_y) {}
 62     double x, y;
 63 };
 64 struct line {
 65     friend istream& operator >>(istream& in, line& rhs) {
 66         in >> rhs.p >> rhs.v;
 67         return in;
 68     }
 69     friend ostream& operator <<(ostream& os, const line& rhs) {
 70         os << rhs.p << ‘ ‘ << rhs.v;
 71         return os;
 72     }
 73     bool operator < (const line& rhs)const {
 74         if (Dcmp(ang - rhs.ang))return ang < rhs.ang;
 75         else return Dcmp(v.det(rhs.v)) == -1;
 76     }
 77     point Get_Point(double t)const {
 78         return p + v * t;
 79     }
 80     bool IsParallel(const line& rhs)const { return Dcmp(v.det(rhs.v)) == 0; }
 81     point p, v;
 82     double ang;
 83     line() {}
 84     line(point _p, point _v) :p(_p), v(_v) { ang = atan2(_v.y, _v.x); }
 85 };
 86 //不想写scanf
 87 inline point read_point() {
 88     double x, y;
 89     scanf("%lf %lf", &x, &y);
 90     return point(x, y);
 91 }
 92 //向量长度
 93 inline double length(const point& a) { return sqrt(a.dot(a)); }
 94 //两点之间长度的平方
 95 inline double Dist(const point& a, const point& b) { return (a - b).dot(a - b); }
 96 //两向量的夹角
 97 inline double Angle(const point& a, const point& b) { return acos(a.dot(b) / (length(a) * length(b))); }
 98 //三点无向面积公式
 99 inline double Area2(const point& a, const point& b, const point& c) { return fabs((b - a).det(c - a)) / 2; }
100 //判断p q两点是否在线段ab的同侧
101 inline bool SameSide(const point& a, const point& b, const point& p, const point& q) {
102     point ab = b - a, ac = p - a, ap = q - a;
103     double res = ab.det(ac) * ab.det(ap);
104     return Dcmp(res) >= 0;
105 }
106 //判断点是否在三角形内
107 inline bool InTriangle(const point& a, const point& b, const point& c, const point& p) { return SameSide(a, b, c, p) && SameSide(b, c, a, p) && SameSide(c, a, b, p); }
108 //逆时针旋转向量
109 inline point Rotate(const point& rhs, double rad) { return point(rhs.x * cos(rad) - rhs.y *sin(rad), rhs.x * sin(rad) + rhs.y * cos(rad)); }
110 //求向量单位法向量
111 inline point Normal(const point& rhs) {
112     if (Dcmp(rhs.x) == 0 && Dcmp(rhs.y) == 0)return point(0.0, 0.0);
113     double len = length(rhs);
114     return point(-rhs.y / len, rhs.x / len);
115 }
116 //求单位向量
117 inline point Unit_vector(const point& rhs) {
118     double len = length(rhs);
119     return point(rhs.x / len, rhs.y / len);
120 }
121 //判断q点是否在线段p1, p1上
122 inline bool On_seg(const point& p1, const point& p2, const point& q) { return Dcmp((p1 - q).det(p2 - q)) == 0 && Dcmp((p1 - q).dot(p2 - q)) < 0; }
123 //判断p是否在直线左边
124 inline bool On_Left(const line& L, const point& p) { return L.v.det(p - L.p) > 0; }
125 //通过4点坐标求两直线交点
126 inline point Intersection_point(const point& p1, const point& p2, const point& q1, const point& q2) {
127     return p1 + (p2 - p1) * ((q2 - q1).det(q1 - p1) / (q2 - q1).det(p2 - p1));
128 }
129 //通过直线上一点及其方向向量求交点
130 inline point Intersection_line(const line& L1, const line& L2) {
131     point u = L1.p - L2.p;
132     double t = L2.v.det(u) / L1.v.det(L2.v);
133     return L1.p + L1.v * t;
134 }
135 //判断两线段是否规范相交
136 inline bool SegmentProperIntersection(const point& a1, const point& a2, const point& b1, const point& b2) {
137     double c1 = (a2 - a1).det(b1 - a1), c2 = (a2 - a1).det(b2 - a1), c3 = (b2 - b1).det(a1 - b1), c4 = (b2 - b1).det(a2 - b1);
138     return Dcmp(c1) * Dcmp(c2) < 0 && Dcmp(c3) * Dcmp(c4) < 0;
139 }
140 //判断两线段是否相交，判断所有情况
141 inline bool SegmentIntersection(const point& a1, const point& a2, const point& b1, const point& b2) {
142     if (Dcmp((a2 - a1).det(b2 - b1)) == 0)return On_seg(a1, a2, b1) || On_seg(a1, a2, b2) || On_seg(b1, b2, a1) || On_seg(b2, b2, a2);
143     double c1 = (a2 - a1).det(b1 - a1), c2 = (a2 - a1).det(b2 - a1), c3 = (b2 - b1).det(a1 - b1), c4 = (b2 - b1).det(a2 - b1);
144     return Dcmp(c1) * Dcmp(c2) < 0 && Dcmp(c3) * Dcmp(c4) < 0;
145 }
146 //求点p到直线ab的垂直距离
147 inline double DistanceToLine(const point& p, const point& a, const point& b) { return 2 * Area2(p, a, b) / length(a - b); }
148 //求点p到线段ab的距离
149 inline double DistanceToSegment(const point& p, const point& a, const point& b) {
150     if (a == b)return length(p - a);
151     point v1 = b - a, v2 = p - a, v3 = p - b;
152     if (Dcmp(v1.dot(v2) < 0))return length(v2);
153     else if (Dcmp(v1.dot(v3) > 0))return length(v3);
154     else return DistanceToLine(p, a, b);
155 }
156 //求点p在直线ab上的投影点
157 inline point GetLineProjection(const point& p, const point& a, const point& b) {
158     point v = b - a;
159     return a + v * (v.dot(p - a) / v.dot(v));
160 }
161 //计算凸多边形面积
162 double ConvexPolygonArea(point* p, int n) {
163     double res = 0;
164     for (int i = 1; i < n - 1; ++i)res += Area2(p[0], p[i], p[i + 1]);
165     return res;
166 }
167 //计算任意多边形面积公式
168 double PolygonArea(point* p, int n) {
169     double res = 0;
170     for (int i = 1; i < n - 1; ++i)res += (p[i] - p[0]).det(p[i + 1] - p[0]);
171     return res / 2;
172 }
173 //判断点是否在多边形内
174 int isPointinPolygon(const point& p, vector<point>& v) {
175     int cnt = 0, n = v.size();
176     for (int i = 0; i < n; ++i) {
177         if (On_seg(v[i], v[(i + 1) % n], p))return -1;
178         int k = Dcmp((v[(i + 1) % n] - v[i]).det(p - v[i]));
179         int d1 = Dcmp(v[i].y - p.y), d2 = Dcmp(v[(i + 1) % n].y - p.y);
180         if (k > 0 && d1 <= 0 && d2 > 0)++cnt;
181         if (k < 0 && d2 <= 0 && d1 > 0)--cnt;
182     }
183     if (cnt)return 1;//内部
184     else return 0;//外部
185 }
186 vector<point> v;
187 //构造凸包
188 inline void Get_convex_hall(point* arr, int n) {
189     v.clear();
190     v.resize(n << 1);
191     sort(arr + 1, arr + 1 + n);
192     int cnt = 0;
193     for (int i = 1; i <= n; ++i) {
194         while (cnt > 1 && (v[cnt - 1] - v[cnt - 2]).det(arr[i] - v[cnt - 1]) <= 0)--cnt;
195         v[cnt++] = arr[i];
196     }
197     for (int i = n - 1, t = cnt; i > 0; --i) {
198         while (cnt > t && (v[cnt - 1] - v[cnt - 2]).det(arr[i] - v[cnt - 1]) <= 0)--cnt;
199         v[cnt++] = arr[i];
200     }
201     v.resize(cnt - 1);
202 }
203 //旋转卡壳
204 double Rotate_Calipers() {
205     double res = 0.0;
206     int n = v.size();
207     v.push_back(v[0]);
208     int now = 1;
209     for (int i = 0; i < n; ++i) {
210         while ((v[now] - v[i + 1]).det(v[i] - v[i + 1]) < (v[now + 1] - v[i + 1]).det(v[i] - v[i + 1]))now = (now + 1) % n;
211         res = max(res, max(Dist(v[now], v[i]), Dist(v[now + 1], v[i + 1])));
212     }
213     v.pop_back();
214     return res;
215 }
216 //求半平面交
217 bool HalfPlaneIntersection(line* arr, int n, vector<point>& v) {
218     deque<line> qLine;
219     deque<point> qPoint;
220     sort(arr + 1, arr + 1 + n);
221     qLine.push_back(arr[1]);
222     for (int i = 2; i <= n; ++i) {
223         if (Dcmp(arr[i].ang - arr[i - 1].ang)) {
224             if (qLine.size() > 1 && (qLine[0].IsParallel(qLine[1]) || qLine[qLine.size() - 1].IsParallel(qLine[qLine.size() - 2])))return false;
225             while (qLine.size() > 1 && !On_Left(arr[i], qPoint[qPoint.size() - 1]))qPoint.pop_back(), qLine.pop_back();
226             while (qLine.size() > 1 && !On_Left(arr[i], qPoint[0]))qPoint.pop_front(), qLine.pop_front();
227             qLine.push_back(arr[i]);
228             if (qLine.size() > 1)qPoint.push_back(Intersection_line(qLine[qLine.size() - 1], qLine[qLine.size() - 2]));
229         }
230     }
231     while (qLine.size() > 1 && !On_Left(qLine[0], qPoint[qPoint.size() - 1]))qLine.pop_back(), qPoint.pop_back();
232     while (qLine.size() > 1 && !On_Left(qLine[qLine.size() - 1], qPoint[0]))qLine.pop_front(), qPoint.pop_front();
233     if (qPoint.size() < 2)return false;
234     qPoint.push_back(Intersection_line(qLine[0], qLine[qLine.size() - 1]));
235     for (int i = 0; i < qPoint.size(); ++i)v.push_back(qPoint[i]);
236     return true;
237 }
238 //------------------------------点与直线部分-----------------------------
239 
240 //------------------------------圆与球部分-------------------------------
241 struct circle {
242     friend istream& operator >>(istream& in, circle& rhs) {
243         in >> rhs.c >> rhs.r;
244         return in;
245     }
246     friend ostream& operator <<(ostream& os, const circle& rhs) {
247         os << rhs.c << " " << rhs.r;
248         return os;
249     }
250     point c;
251     double r;
252     circle() {}
253     circle(point _c, double _r) :c(_c), r(_r) {}
254     //通过圆心角求圆上点
255     point Get_point(double a)const { return point(c.x + r * cos(a), c.y + r * sin(a)); }
256 };
257 //求极角
258 inline double polar_angle(const point& rhs) { return atan2(rhs.y, rhs.x); }
259 //将角度转换为弧度
260 inline double torad(double deg) { return deg / 180. * pi; }
261 //解方程求直线与圆的交点
262 int GetLineCircleIntersection(const line& L, const circle& C, double& t1, double& t2, vector<point>& v) {
263     double a = L.v.x, b = L.p.x - C.c.x, c = L.v.y, d = L.p.y - C.c.y;
264     double e = a * a + c * c, f = 2 * (a * b + c * d), g = b * b + d * d - C.r * C.r;
265     double delta = f * f - 4 * e * g;
266     if (Dcmp(delta) < 0)return 0;
267     if (Dcmp(delta) == 0) {
268         t1 = t2 = -f / (2 * e);
269         v.push_back(L.Get_Point(t1));
270         return 1;
271     }
272     t1 = (-f - sqrt(delta)) / (2 * e);
273     t2 = (-f + sqrt(delta)) / (2 * e);
274     v.push_back(L.Get_Point(t1));
275     v.push_back(L.Get_Point(t2));
276     return 2;
277 }
278 //勾股定理求直线与圆的交点
279 int GetLineCircleIntersection(const point& A, const point& B, const circle&C, double& L, vector<point>& v) {
280     double d = DistanceToLine(C.c, A, B);
281     if (Dcmp(C.r - d) < 0)return 0;
282     point p = GetLineProjection(C.c, A, B);
283     double len = sqrt(C.r * C.r - d * d);
284     if (Dcmp(len) == 0) {
285         L = 0;
286         v.push_back(p);
287         return 1;
288     }
289     point unit = Unit_vector(B - A);
290     L = len;
291     v.push_back(p - unit * len);
292     v.push_back(p + unit * len);
293     return 2;
294 }
295 //求两圆交点
296 int GetCircleCircleIntersection(const circle& C1, const circle& C2, vector<point>& v) {
297     double d = length(C1.c - C2.c);
298     if (Dcmp(d) == 0) {
299         if (Dcmp(C1.r - C2.r) == 0)return -1;//两圆重合
300         return 0;
301     }
302     if (Dcmp(C1.r + C2.r - d) < 0 || Dcmp(fabs(C1.r - C1.r) - d) > 0)return 0;
303     double a = polar_angle(C2.c - C1.c);
304     double da = acos((C1.r * C1.r + d * d - C2.r * C2.r) / (2 * C1.r * d));
305     point p1 = C1.Get_point(a - da), p2 = C1.Get_point(a + da);
306     v.push_back(p1);
307     if (p1 == p2)return 1;
308     v.push_back(p2);
309     return 2;
310 }
311 //求点到圆的切线
312 int GetTangents(const point& p, const circle& c, vector<point>& v) {
313     point u = c.c - p;
314     double dist = length(u);
315     if (dist < c.r)return 0;
316     else if (Dcmp(dist - c.r) == 0) {
317         v.push_back(Rotate(u, pi / 2));
318         return 1;
319     }
320     else {
321         double ang = asin(c.r / dist);
322         v.push_back(Rotate(u, -ang));
323         v.push_back(Rotate(u, +ang));
324         return 2;
325     }
326 }
327 //求两圆的切线，返回切线个数
328 int GetTangents(circle A, circle B, vector<point>& a, vector<point>& b) {
329     if (Dcmp(A.r - B.r) < 0) {
330         swap(A, B);
331         swap(a, b);
332     }
333     double d2 = Dist(A.c, B.c), rdiff = A.r - B.r, rsum = A.r + B.r;
334     if (Dcmp(d2 - rdiff * rdiff) == -1)return 0;//内含
335     double base = atan2(B.c.y - A.c.y, B.c.x - A.c.x);
336     if (Dcmp(d2) == 0 && Dcmp(A.r - B.r) == 0)return -1;
337     if (Dcmp(d2 - rdiff * rdiff) == 0) {
338         a.push_back(A.Get_point(base));
339         b.push_back(B.Get_point(base));
340         return 1;
341     }//内切
342     double ang = acos((A.r - B.r) / sqrt(d2));
343     a.push_back(A.Get_point(base + ang)), a.push_back(A.Get_point(base - ang));
344     b.push_back(B.Get_point(base + ang)), b.push_back(B.Get_point(base - ang));
345     if (Dcmp(d2 - rsum * rsum) == 0) {
346         a.push_back(A.Get_point(base));
347         b.push_back(B.Get_point(pi + base));
348     }
349     else if (Dcmp(d2 - rsum * rsum) > 0) {
350         ang = acos((A.r + B.r) / sqrt(d2));
351         a.push_back(A.Get_point(base + ang)), a.push_back(A.Get_point(base - ang));
352         b.push_back(B.Get_point(pi + base + ang)), b.push_back(B.Get_point(pi + base - ang));
353     }
354     return (int)a.size();
355 }
356 //求三角形外接圆
357 circle CircumscribedCircle(const point& a, const point& b, const point& c) {
358     double Bx = b.x - a.x, By = b.y - a.y, Cx = c.x - a.x, Cy = c.y - a.y;
359     double D = 2 * (Bx * Cy - Cx * By);
360     double cx = (Cy * (Bx * Bx + By * By) - By *(Cx * Cx + Cy * Cy)) / D + a.x;
361     double cy = (Bx * (Cx * Cx + Cy * Cy) - Cx *(Bx * Bx + By * By)) / D + a.y;
362     point p = point(cx, cy);
363     return circle(p, length(a - p));
364 }
365 //求三角形内接圆
366 circle InscribedCircle(const point& a, const point& b, const point& c) {
367     double lena = length(b - c), lenb = length(c - a), lenc = length(a - b);
368     point p = (a * lena + b * lenb + c * lenc) / (lena + lenb + lenc);
369     return circle(p, DistanceToLine(p, a, b));
370 }
371 //------------------------------圆与球部分-------------------------------