1. 点、线、凸边形
1 /*******************************************************
2 二维几何基础
3 【注意】数组下标从1开始。
4 *******************************************************/
5
6 #include <iostream>
7 #include <cstdio>
8 #include <cstdlib>
9 #include <cmath>
10 #include <iomanip>
11 #include <algorithm>
12 #include <string>
13 #define re register
14 #define il inline
15 #define ll long long
16 #define ld long double
17 using namespace std;
18 const ll MAXN = 1e5+5;
19 const ll INF = 1e9;
20 const ld EPS = 1e-9;
21
22 //点坐标
23 struct POINT
24 {
25 ld x, y;
26 POINT() : x(0), y(0) {}
27 POINT(ld _x, ld _y) : x(_x), y(_y) {}
28 };
29 //向量
30 typedef POINT VECTOR;
31
32 POINT xy[MAXN]; //顺时针多边形顶点存储
33 POINT xyB[MAXN]; //判断点存储
34
35
36 //符号函数
37 ll sgn(ld x) {return x < -EPS ? -1 : x > EPS;}
38 //向量+向量=向量，点+向量=点
39 VECTOR operator + (VECTOR a, VECTOR b) {return {a.x+b.x,a.y+b.y};}
40 //向量-向量=向量，点-点=向量
41 VECTOR operator - (VECTOR a, VECTOR b) {return {a.x-b.x,a.y-b.y};}
42 //向量*数=向量
43 VECTOR operator * (VECTOR a, ld k) {return {a.x*k,a.y*k};}
44 VECTOR operator * (ld k, VECTOR a) {return {a.x*k,a.y*k};}
45 //向量/数=向量
46 VECTOR operator / (VECTOR a, ld k) {return {a.x/k,a.y/k};}
47 //向量==向量，点==点
48 bool operator == (const VECTOR a, const VECTOR b) {return !sgn(a.x-b.x) && !sgn(a.y-b.y);}
49 //向量偏序关系（先x轴再y轴）
50 bool operator < (const VECTOR a, const VECTOR b) {return !sgn(a.x-b.x) ? a.y < b.y : a.x < b.x;}
51 //向量旋转（逆时针）
52 VECTOR rot(VECTOR a, ld sita) {return {a.x*cos(sita)-a.y*sin(sita),a.x*sin(sita)+a.y*cos(sita)};}
53 //取下端点
54 POINT min(POINT p1, POINT p2) {return p1.y<p2.y ? p1:p2;}
55 //取上端点
56 POINT max(POINT p1, POINT p2) {return p1.y<p2.y ? p2:p1;}
57 //点积
58 ld dot(POINT p1, POINT p2, POINT p) {return (p.x-p1.x)*(p2.x-p1.x)+(p.y-p1.y)*(p2.y-p1.y);} //（三点式：共p1端点）
59 ld dot(VECTOR a, VECTOR b) {return a.x*b.x+a.y*b.y;} //向量式
60 //叉积
61 ld cross(POINT p1, POINT p2, POINT p) {return (p.x-p1.x)*(p2.y-p1.y)-(p.y-p1.y)*(p2.x-p1.x);} //（三点式：共p1端点）
62 ld cross(VECTOR a, VECTOR b) {return a.x*b.y-a.y*b.x;} //向量式（aXb）
63 //向量长度
64 ld vlen(VECTOR a) {return sqrt(dot(a,a));}
65 //向量夹角余弦值
66 ld vcos(VECTOR a, VECTOR b) {return dot(a,b)/(vlen(a)*vlen(b));}
67 //向量夹角正弦值
68 ld vsin(VECTOR a, VECTOR b) {return fabs(cross(a,b))/(vlen(a)*vlen(b));}
69 //向量极角
70 ld vagl(VECTOR a, VECTOR b) {return acos(dot(a,b)/(vlen(a)*vlen(b)));} //向量间
71 ld vagl(VECTOR a) {return acos(a.x/vlen(a));}
72 //求直线斜率【注意】确保斜率存在
73 ld slope(VECTOR a) {return a.y/a.x;} //向量式
74 ld slope(POINT p, POINT q) {return (p.y-q.y)/(p.x-q.x);} //两点式
75 //单位向量【注意】确保非零向量
76 VECTOR unit(VECTOR a) {return a/vlen(a);}
77 //两直线交点
78 POINT intersectline(POINT p, VECTOR v, POINT q, VECTOR w) {return p+v*cross(w,p-q)/cross(v,w);} //（参数式：P=P0+t*v，P0为直线上某一点，v为直线的方向向量）
79 //点在直线上的投影
80 POINT proline(POINT a, POINT b, POINT p) {return a+(b-a)*(dot(b-a,p-a)/dot(b-a,b-a));}
81 //点关于直线的对称点
82 POINT refline(POINT a, POINT b, POINT p) {return proline(a,b,p)*2-p;}
83 //判断两直线是否平行
84 bool parallel(POINT p1, POINT p2, POINT q1, POINT q2) {return !sgn(cross(p2-p1,q2-q1)) && sgn(cross(p1-q1,p2-q1));}
85 //判断两直线重合
86 bool superposition(POINT p1, POINT p2, POINT q1, POINT q2) {return !sgn(cross(p2-p1,q2-q1)) && !sgn(cross(p1-q1,p2-q1));}
87 //点到直线距离
88 ld disline(POINT a, POINT b, POINT p) {return fabs(cross(b-a,p-a))/vlen(b-a);} //不取绝对值得到的是有向距离
89 //点到线段距离（两种情况：点的投影在线段上，则为垂直距离；点的投影不在线段上，则为到两端距离的最小值）
90 ld disseg(POINT a, POINT b, POINT p) {return a==b ? vlen(p-a):dot(b-a,p-a)<0 ? vlen(p-a):dot(b-a,p-b)>0 ? vlen(p-b):fabs(cross(b-a,p-a))/vlen(b-a);}
91 //线段相交判断（严格）
92 bool intersectseg(POINT p1, POINT p2, POINT q1, POINT q2) {return cross(p1-q1,q2-q1)*cross(p2-q1,q2-q1)<0 && cross(q1-p1,p2-p1)*cross(q2-p1,p2-p1)<0;}
93 //判断点p(x,y)是否在线段p1p2上（两种情况：1.包括端点；2.不包含端点）
94 bool onseg(POINT p1, POINT p2, POINT p) {return !sgn(cross(p1,p2,p)) && sgn(dot(p,p1,p2))<=0;} //包含端点
95 bool onseg_strict(POINT p1, POINT p2, POINT p) {return !sgn(cross(p1,p2,p)) && sgn(dot(p,p1,p2))<0;} //不包含端点
96
97 //点射法判断点是否在多边形内部（边界也算在内部）
98 //复杂度O(n)（不论凹凸）
99 bool inpolygon_dot(ld x, ld y, ll n)
100 {
101 ll cnt = 0;
102 for(re ll i = 0; i < n; ++i)
103 {
104 POINT p1 = min(xy[i+1], xy[(i+1)%n+1]); //取下端点
105 POINT p2 = max(xy[i+1], xy[(i+1)%n+1]); //取上端点
106 //特判点在线段上
107 if(onseg(p1,p2,{x,y})) return true;
108 //从点(x,y)向x反方向引出射线
109 //计算点射到的多边形的边数边数
110 if(sgn(p1.y-y)<0 && sgn(y-p2.y) <=0 && sgn(cross(p1,p2,{x,y}))>0) ++cnt;
111 }
112 if(cnt%2) return true;
113 else return false;
114 }
115
116 //二分法判断点是否在多边形内部（不包括边界）
117 //复杂度O(logn)（要求凸多边形）
118 bool inpolygon_bis(ld x, ld y, ll n)
119 {
120 POINT p = {x,y};
121 ll l = 2, r = n;
122 //判断是否在幅角最小和最大的两条边上
123 if(onseg(xy[1],xy[l],p) || onseg(xy[1],xy[n],p)) return false;
124 //二分法判断是否在多边形内部
125 while(l < r)
126 {
127 ll mid = (l+r)>>1; //【注意】如此取中点最终：r==l或r==l-1（只有此两种情况）
128 ll d = sgn(cross(xy[mid]-xy[1],p-xy[1]));
129 if(d < 0) l = mid+1;
130 else if(d > 0) r = mid-1;
131 else
132 {
133 if(onseg_strict(xy[1],xy[mid],p))
134 {
135 return true;
136 }
137 else
138 {
139 return false;
140 }
141 }
142 }
143 //判断在最终二分得到的对角线两端的边是否都满足严格在内的条件
144 if(l >= r && (sgn(cross(xy[l]-xy[l-1],p-xy[l-1]))>=0 || sgn(cross(xy[l%n+1]-xy[l],p-xy[l]))>=0))
145 {
146 return false;
147 }
148 return true;
149 }
150
151 //计算多边形面积（凹凸均可）
152 ld polygonarea(ll n)
153 {
154 ld aera = 0;
155 for(re ll i = 0; i < n; ++i)
156 {
157 aera += cross(xy[i+1], xy[(i+1)%n+1]);
158 }
159 //计算出来的aera为有向面积
160 return fabs(aera)/2;
161 }
162
163 int main()
164 {
165 std::ios::sync_with_stdio(false);
166 //判断直线平行/重合/相交
167 ll n;
168 cin >> n;
169 cout << "INTERSECTING LINES OUTPUT" << endl;
170 for(re ll i = 1; i <= n; ++i)
171 {
172 POINT p1, p2, p3, p4;
173 cin >> p1.x >> p1.y >> p2.x >> p2.y;
174 cin >> p3.x >> p3.y >> p4.x >> p4.y;
175 if(parallel(p1,p2,p3,p4)) cout << "NONE" << endl;
176 else if(superposition(p1,p2,p3,p4)) cout << "LINE" << endl;
177 else
178 {
179 POINT q = intersectline(p1,p2-p1,p3,p4-p3);
180 cout << "POINT " << fixed << setprecision(2) << q.x << " " << q.y << endl;
181 }
182 }
183 cout << "END OF OUTPUT" << endl;
184 return 0;
185 /*
186 //计算多边形的面积
187 ll n;
188 while(true)
189 {
190 std::cin >> n;
191 if(!n) break;
192 for(re ll i = 1; i <= n; ++i)
193 {
194 std::cin >> xy[i].x;
195 std::cin >> xy[i].y;
196 }
197 std::cout << std::fixed << std::setprecision(1) << polygonarea(n) << std::endl;
198 }
199 return 0;
200 */
201 /*
202 //判断点是否在多边形内（不包括边界）
203 ll n;
204 std::cin >> n;
205 for(re ll i = 1; i <= n; ++i)
206 {
207 std::cin >> xy[i].x;
208 std::cin >> xy[i].y;
209 }
210 ll m;
211 cin >> m;
212 for(re ll i = 1; i <= m; ++i)
213 {
214 cin >> xyB[i].x;
215 cin >> xyB[i].y;
216 }
217 for(re ll i = 1; i <= m; ++i)
218 {
219 if(!inpolygon_bis(xyB[i].x, xyB[i].y, n))
220 {
221 printf("NO\n");
222 return 0;
223 }
224 }
225 printf("YES\n");
226 return 0;
227 */
228 /*
229 //判断点是否在多边形内（包括边界）
230 ll N, M, k = 0;
231 while(!(std::cin >> N).eof())
232 {
233 if(!N) break;
234 std::cin >> M;
235 if(k) printf("\n");
236 ++k;
237 for(re ll i = 1; i <= N; ++i)
238 {
239 std::cin >> xy[i].x;
240 std::cin >> xy[i].y;
241 //printf("(%f,%f)\n", xy[i].x, xy[i].y);
242 }
243 printf("Problem %lld:\n", k);
244 for(re ll i = 1; i <= M; ++i)
245 {
246 ld x, y;
247 std::cin >> x;
248 std::cin >> y;
249 //printf("(%f,%f)\n", x, y);
250 if(inpolygon(x, y, N))
251 {
252 printf("Within\n");
253 }
254 else
255 {
256 printf("Outside\n");
257 }
258 }
259 }
260 return 0;
261 */
262 }
2. 圆
1 /*****************************************************************
2 圆
3 【注意】三角与反三角最多用一次（损失精度大）。
4 *****************************************************************/
5
6 #include <bits/stdc++.h>
7 #include <iomanip>
8 #include <stack>
9 #include <cstdio>
10 #include <cmath>
11 #include <algorithm>
12 #define re register
13 #define il inline
14 #define ll long long
15 #define ld long double
16 using namespace std;
17 const ld PI = 3.14159265358979323846;
18 const ll MAXN = 5e5+5;
19 const ld INF = 1e9;
20 const ld EPS = 1e-10;
21
22 //点坐标
23 struct POINT
24 {
25 ld x, y;
26 POINT() : x(0), y(0) {}
27 POINT(ld _x, ld _y) : x(_x), y(_y) {}
28 };
29 //向量
30 typedef POINT VECTOR;
31
32 //符号函数
33 ll sgn(ld x) {return x < -EPS ? -1 : x > EPS;}
34 //向量+向量=向量，点+向量=点
35 VECTOR operator + (VECTOR a, VECTOR b) {return {a.x+b.x,a.y+b.y};}
36 //向量-向量=向量，点-点=向量
37 VECTOR operator - (VECTOR a, VECTOR b) {return {a.x-b.x,a.y-b.y};}
38 //向量*数=向量
39 VECTOR operator * (VECTOR a, ld k) {return {a.x*k,a.y*k};}
40 VECTOR operator * (ld k, VECTOR a) {return {a.x*k,a.y*k};}
41 //向量/数=向量
42 VECTOR operator / (VECTOR a, ld k) {return {a.x/k,a.y/k};}
43 //向量==向量，点==点
44 bool operator == (const VECTOR a, const VECTOR b) {return !sgn(a.x-b.x) && !sgn(a.y-b.y);}
45 //向量偏序关系（先x轴再y轴）
46 bool operator < (const VECTOR a, const VECTOR b) {return !sgn(a.x-b.x) ? a.y < b.y : a.x < b.x;}
47 //取下端点
48 POINT min(POINT p1, POINT p2) {return p1.y<p2.y ? p1:p2;}
49 //取上端点
50 POINT max(POINT p1, POINT p2) {return p1.y<p2.y ? p2:p1;}
51 //点积
52 ld dot(POINT p1, POINT p2, POINT p) {return (p.x-p1.x)*(p2.x-p1.x)+(p.y-p1.y)*(p2.y-p1.y);} //（三点式：共p1端点）
53 ld dot(VECTOR a, VECTOR b) {return a.x*b.x+a.y*b.y;} //向量式
54 //叉积
55 ld cross(POINT p1, POINT p2, POINT p) {return (p.x-p1.x)*(p2.y-p1.y)-(p.y-p1.y)*(p2.x-p1.x);} //（三点式：共p1端点）
56 ld cross(VECTOR a, VECTOR b) {return a.x*b.y-a.y*b.x;} //向量式（aXb）
57 //向量长度
58 ld vlen(VECTOR a) {return sqrt(dot(a,a));}
59 //向量旋转（逆时针）
60 VECTOR rot(VECTOR a, ld sita) {return {a.x*cos(sita)-a.y*sin(sita),a.x*sin(sita)+a.y*cos(sita)};}
61 //取单位向量
62 VECTOR unit(VECTOR a) {return a/vlen(a);}
63 //取正交单位向量
64 VECTOR norm(VECTOR a) {return VECTOR(-a.y,a.x)/vlen(a);}
65 //向量夹角余弦值
66 ld vcos(VECTOR a, VECTOR b) {return dot(a,b)/(vlen(a)*vlen(b));}
67 //向量夹角正弦值
68 ld vsin(VECTOR a, VECTOR b) {return fabs(cross(a,b))/(vlen(a)*vlen(b));}
69 //向量极角
70 //【注意】使用系统反三角函数返回nan是因为超出了定义域（大多数由于精度造成）
71 //【注意】故使用系统反三角函数需处理边界
72 ld vagl(VECTOR a, VECTOR b)
73 {
74 ld d = dot(a,b)/(vlen(a)*vlen(b));
75 if(sgn(d-1) == 0) return 0;
76 else if(sgn(d+1) == 0) return PI;
77 else return acos(d);
78 } //向量间[0,PI]
79 ld vagl(VECTOR a) {return atan2(a.y, a.x);} //单向量[-PI,PI]
80
81 //直线
82 struct LINE
83 {
84 POINT p; //直线上定点
85 VECTOR v; //直线方向向量
86 LINE() {}
87 //点向式
88 LINE(POINT _p, VECTOR _v):p(_p), v(_v) {}
89 //点角式
90 LINE(POINT _p, ld a):p(_p), v(VECTOR(cos(a), sin(a))) {}
91 //通过参数t唯一确定直线上点坐标（点向式：P = p + t*v）
92 POINT point(ld t) {return p + t*v;}
93 //平移向量r
94 LINE trans(VECTOR r) {return LINE(p+r,v);}
95 };
96
97 //点到直线距离
98 ld PLdist(POINT p, LINE L) {return fabs(sin(vagl(p-L.p,L.v))*vlen(p-L.p));}
99
100 //圆
101 struct CIRCLE
102 {
103 POINT c; //圆心坐标
104 ld r; //半径
105 CIRCLE() {}
106 CIRCLE(POINT _c, ld _r):c(_c),r(_r) {}
107 //通过圆心角唯一确定圆上点坐标
108 POINT point(ld a) //a为圆心角
109 {
110 return POINT(c.x+cos(a)*r, c.y+sin(a)*r);
111 }
112 };
113
114 //两点中垂线
115 //返回中垂线的个数
116 ll PerpendicularBisector(POINT p1, POINT p2, LINE &L)
117 {
118 if(p1 == p2) return 0;
119 L.p = (p1+p2)/2;
120 L.v = rot(p2-p1, PI/2);
121 return 1;
122 }
123
124 //两直线交点
125 POINT LLitesct;
126
127 //直线与直线的交点
128 //返回交点个数
129 ll LineLineIntersection(LINE L1, LINE L2, POINT &p)
130 {
131 if(sgn(cross(L1.v, L2.v)) == 0)
132 {
133 if(sgn(cross(L2.p-L1.p,L2.v) == 0)) return -1; //两直线重合
134 else return 0; //两直线平行
135 }
136 //以L1为准
137 ld t1 = -cross(L2.p-L1.p, L2.v)/cross(L2.v, L1.v);
138 p = L1.p + t1*L1.v;
139 //以L2为准
140 //ld t2 = -cross(L1.p-L2.p, L1.v)/cross(L1.v, L2.v);
141 //p = L2.p + t2*L2.v;
142 return 1;
143 }
144
145 //直线与圆交点数组
146 vector <POINT> LCitesct;
147
148 //直线与圆的交点
149 //解方程：(at+b)^2+(ct+d)^2 = r^2;
150 //化简得：et^2+ft+g = 0
151 //返回交点个数
152 ll LineCircleIntersection(LINE L, CIRCLE C, ld &t1, ld &t2, vector <POINT>& sol)
153 {
154 ld a = L.v.x, b = L.p.x-C.c.x, c = L.v.y, d = L.p.y-C.c.y;
155 ld e = a*a + c*c, f = 2*(a*b + c*d), g = b*b + d*d - C.r*C.r;
156 ld delta = f*f - 4*e*g;
157 if(sgn(delta) < 0) //相离
158 return 0;
159 else if(sgn(delta) == 0) //相切
160 {
161 t1 = t2 = -f/(2*e);
162 sol.push_back(L.point(t1));
163 return 1;
164 }
165 else
166 {
167 t1 = (-f - sqrt(delta))/(2*e);
168 t2 = (-f + sqrt(delta))/(2*e);
169 sol.push_back(L.point(t1));
170 sol.push_back(L.point(t2));
171 return 2;
172 }
173 }
174 ll LineCircleIntersection(LINE L, CIRCLE C, POINT *p)
175 {
176 ld t1, t2;
177 ld a = L.v.x, b = L.p.x-C.c.x, c = L.v.y, d = L.p.y-C.c.y;
178 ld e = a*a + c*c, f = 2*(a*b + c*d), g = b*b + d*d - C.r*C.r;
179 ld delta = f*f - 4*e*g;
180 if(sgn(delta) < 0) //相离
181 return 0;
182 else if(sgn(delta) == 0) //相切
183 {
184 t1 = t2 = -f/(2*e);
185 p[0] = L.point(t1);
186 return 1;
187 }
188 else
189 {
190 t1 = (-f - sqrt(delta))/(2*e);
191 t2 = (-f + sqrt(delta))/(2*e);
192 p[0] = L.point(t1);
193 p[1] = L.point(t2);
194 return 2;
195 }
196 }
197
198 //圆与圆的交点数组
199 vector <POINT> CCitesct;
200
201 //两圆交点
202 ll CircleCircleIntersection(CIRCLE C1, CIRCLE C2, vector <POINT>& sol)
203 {
204 ld d = vlen(C1.c-C2.c);
205 if(sgn(d) == 0)
206 {
207 if(sgn(C1.r-C2.r) == 0) return -1; //两圆重合
208 else return 0; //同心相含
209 }
210 if(sgn(C1.r+C2.r-d) < 0) return 0; //相离
211 if(sgn(fabs(C1.r-C2.r)-d) > 0) return 0; //异心相含
212
213 //设两圆交点为p1,p2
214 ld a = vagl(C2.c - C1.c); //向量C1C2的极角
215 ld b = acos((C1.r*C1.r + d*d - C2.r*C2.r)/(2*C1.r*d)); //向量C1C2到C1P1（或C1P2）的夹角
216 POINT p1 = C1.point(a-b);
217 POINT p2 = C1.point(a+b);
218 sol.push_back(p1);
219 if(p1 == p2) return 1; //相切
220 sol.push_back(p2); //相交
221 return 2;
222 }
223
224 //圆切线方向向量
225 VECTOR LCtgt[2];
226
227 //过定点作圆的切线
228 ll Tangents(POINT p, CIRCLE C, VECTOR *v)
229 {
230 VECTOR u = C.c - p;
231 ld dist = vlen(u);
232 if(dist < C.r) return 0; //点在圆内，没有切线
233 else if(sgn(dist-C.r) == 0) //点在圆上，一条切线
234 {
235 v[0] = rot(u, PI/2);
236 return 1;
237 }
238 else //点在圆外，两条切线
239 {
240 ld a = asin(C.r/dist);
241 v[0] = rot(u, -a);
242 v[1] = rot(u, +a);
243 return 2;
244 }
245 }
246
247 //两圆公切点（分别在两圆上）
248 POINT CCtgt1[4], CCtgt2[4];
249
250 //两圆公切线
251 ll Tangents(CIRCLE C1, CIRCLE C2, POINT *a, POINT *b)
252 {
253 ll cnt = 0;
254 if(C1.r < C2.r) {swap(C1, C2); swap(a, b);}
255 ld d = dot(C1.c-C2.c, C1.c-C2.c);
256 ld rdiff = C1.r - C2.r;
257 ld rsum = C1.r + C2.r;
258 if(d < rdiff*rsum) return 0; //相含
259 ld base = atan2(C2.c.y-C1.c.y, C2.c.x-C1.c.x); //向量C1C2极角[-PI,PI]
260 if(sgn(d) == 0 && sgn(C1.r-C2.r) == 0) return -1; //重合
261 if(sgn(d-rdiff*rdiff) == 0) //内切
262 {
263 a[cnt] = C1.point(base);
264 b[cnt++] = C2.point(base);
265 return 1;
266 }
267 ld ang = acos((C1.r-C2.r)/sqrt(d)); //夹角[0,PI]
268 a[cnt] = C1.point(base+ang);
269 b[cnt++] = C2.point(base+ang);
270 a[cnt] = C1.point(base-ang);
271 b[cnt++] = C2.point(base-ang);
272 if(sgn(d-rsum*rsum)) //外切
273 {
274 a[cnt] = C1.point(base);
275 b[cnt++] = C2.point(base+PI);
276 }
277 else if(d > rsum*rsum) //相离
278 {
279 ld ang = acos((C1.r+C2.r)/sqrt(d));
280 a[cnt] = C1.point(base+ang);
281 b[cnt++] = C2.point(base+ang);
282 a[cnt] = C1.point(base-ang);
283 b[cnt++] = C2.point(base-ang);
284 }
285 return cnt; //切点个数
286 }
287
288 //三点外接圆
289 //返回外接圆的个数
290 ll CircumscribedCircle(POINT a, POINT b, POINT c, CIRCLE &C)
291 {
292 if(sgn(cross(a-b, b-c)) == 0) return 0; //三点共线
293 LINE L1, L2;
294 PerpendicularBisector(a, b, L1);
295 PerpendicularBisector(b, c, L2);
296 POINT p;
297 LineLineIntersection(L1, L2, p);
298 C.c = p;
299 C.r = vlen(p-a);
300 return 1;
301 }
302
303 //三点内接圆
304 //返回内接圆的个数
305 ll InscribedCircle(POINT a, POINT b, POINT c, CIRCLE &C)
306 {
307 if(sgn(cross(a-b, b-c)) == 0) return 0; //三点共线
308 LINE L1 = LINE(a, (vagl(b-a)+vagl(c-a))/2), L2 = LINE(b, (vagl(a-b)+vagl(c-b))/2);
309 POINT p;
310 LineLineIntersection(L1, L2, p);
311 C.c = p;
312 C.r = PLdist(p,LINE(a,b-a));
313 return 1;
314 }
315
316
317 int main()
318 {
319 //UVA-12304
320 //ios::sync_with_stdio(false);
321 string str;
322 while(cin >> str)
323 {
324 if(str == "CircumscribedCircle")
325 {
326 POINT a, b, c;
327 cin >> a.x >> a.y >> b.x >> b.y >> c.x >> c.y;
328 CIRCLE C;
329 CircumscribedCircle(a,b,c,C);
330 cout << fixed << "(" << C.c.x << "," << C.c.y << "," << C.r << ")" << endl;
331 }
332 else if(str == "InscribedCircle")
333 {
334 POINT a, b, c;
335 cin >> a.x >> a.y >> b.x >> b.y >> c.x >> c.y;
336 CIRCLE C;
337 InscribedCircle(a,b,c,C);
338 cout << fixed << "(" << C.c.x << "," << C.c.y << "," << C.r << ")" << endl;
339 }
340 //求已知圆的切线的极角
341 else if(str == "TangentLineThroughPoint")
342 {
343 CIRCLE C;
344 POINT p;
345 cin >> C.c.x >> C.c.y >> C.r >> p.x >> p.y;
346 VECTOR v[2];
347 ll cnt = Tangents(p,C,v);
348 cout << "[";
349 if(cnt == 1)
350 {
351 ld r = vagl(v[0]);
352 ld a = (sgn(r) < 0 ? PI+r : sgn(r-PI) == 0 ? 0 : r)*180/PI;
353 cout << fixed << a;
354 }
355 else if(cnt == 2)
356 {
357 ld mxr = vagl(v[0]);
358 ld mir = vagl(v[1]);
359 ld mxa = (sgn(mxr) < 0 ? PI+mxr : sgn(mxr-PI) == 0 ? 0 : mxr)*180/PI;
360 ld mia = (sgn(mir) < 0 ? PI+mir : sgn(mir-PI) == 0 ? 0 : mir)*180/PI;
361 if(mxa < mia) swap(mxa,mia);
362 cout << fixed << mia << "," << mxa;
363 }
364 cout << "]" << endl;
365 }
366 //求过定点且与已知直线相切的圆的圆心
367 else if(str == "CircleThroughAPointAndTangentToALineWithRadius")
368 {
369 CIRCLE P;
370 POINT p1, p2;
371 cin >> P.c.x >> P.c.y >> p1.x >> p1.y >> p2.x >> p2.y >> P.r;
372 VECTOR v = rot(p2-p1,PI/2);
373 ld t = P.r/vlen(v);
374 LINE L1 = LINE(p1+t*v,p2-p1), L2 = LINE(p1-t*v,p2-p1);
375 ll cnt = 0;
376 ld t1, t2;
377 vector <POINT> C;
378 cnt += LineCircleIntersection(L1,P,t1,t2,C);
379 cnt += LineCircleIntersection(L2,P,t1,t2,C);
380 sort(C.begin(),C.end());
381 cout << "[";
382 if(cnt == 1) cout << fixed << "(" << C[0].x << "," << C[0].y << ")";
383 else if(cnt == 2) cout << fixed << "(" << C[0].x << "," << C[0].y << ")" << "," << "(" << C[1].x << "," << C[1].y << ")";
384 cout << "]" << endl;
385 }
386 //求给定半径且与两条不相交直线相切的圆的圆心
387 else if(str == "CircleTangentToTwoLinesWithRadius")
388 {
389 POINT p1, p2, p3, p4;
390 ld r;
391 cin >> p1.x >> p1.y >> p2.x >> p2.y >> p3.x >> p3.y >> p4.x >> p4.y >> r;
392 LINE L1 = LINE(p1,p2-p1), L2 = LINE(p3,p4-p3);
393 VECTOR v1 = norm(p2-p1), v2 = norm(p4-p3);
394 LINE L11 = L1.trans(r*v1), L12 = L1.trans(-r*v1);
395 LINE L21 = L2.trans(r*v2), L22 = L2.trans(-r*v2);
396 vector <POINT> C;
397 POINT p;
398 LineLineIntersection(L11,L21,p);
399 C.push_back(p);
400 LineLineIntersection(L11,L22,p);
401 C.push_back(p);
402 LineLineIntersection(L12,L21,p);
403 C.push_back(p);
404 LineLineIntersection(L12,L22,p);
405 C.push_back(p);
406 sort(C.begin(),C.end());
407 cout << "[";
408 cout << fixed << "(" << C[0].x << "," << C[0].y << ")" << "," << "(" << C[1].x << "," << C[1].y << ")" << ",";
409 cout << fixed << "(" << C[2].x << "," << C[2].y << ")" << "," << "(" << C[3].x << "," << C[3].y << ")";
410 cout << "]" << endl;
411
412 }
413 else if(str == "CircleTangentToTwoDisjointCirclesWithRadius")
414 {
415 CIRCLE C1, C2;
416 ld r;
417 cin >> C1.c.x >> C1.c.y >> C1.r >> C2.c.x >> C2.c.y >> C2.r >> r;
418 C1.r += r, C2.r += r;
419 vector <POINT> p;
420 ll cnt = CircleCircleIntersection(C1,C2,p);
421 sort(p.begin(),p.end());
422 cout << "[";
423 if(cnt == 1)
424 {
425 cout << fixed << "(" << p[0].x << "," << p[0].y << ")";
426 }
427 else if(cnt == 2)
428 {
429 cout << fixed << "(" << p[0].x << "," << p[0].y << ")" << "," << "(" << p[1].x << "," << p[1].y << ")";
430 }
431 cout << "]" << endl;
432 }
433 }
434 return 0;
435 }
3. 半平面
1 /*****************************************************************
2 半平面交
3 【注意】三角与反三角最多用一次（损失精度大）。
4 【注意】数组下标一律从1开始。
5 *****************************************************************/
6
7 //#include <bits/stdc++.h>
8 #include <iostream>
9 #include <cstring>
10 #include <cstdlib>
11 #include <iomanip>
12 #include <stack>
13 #include <cstdio>
14 #include <cmath>
15 #include <algorithm>
16 #define re register
17 #define il inline
18 #define ll int
19 #define ld double
20 using namespace std;
21 const ld PI = 3.14159265358979323846;
22 const ll MAXN = 5e5+5;
23 const ld INF = 1e9;
24 const ld EPS = 1e-8;
25
26 //点坐标
27 struct POINT
28 {
29 ld x, y;
30 POINT() : x(0), y(0) {}
31 POINT(ld _x, ld _y) : x(_x), y(_y) {}
32 };
33 //向量
34 typedef POINT VECTOR;
35
36 //符号函数
37 ll sgn(ld x) {return x < -EPS ? -1 : x > EPS;}
38 //向量+向量=向量，点+向量=点
39 VECTOR operator + (VECTOR a, VECTOR b) {return {a.x+b.x,a.y+b.y};}
40 //向量-向量=向量，点-点=向量
41 VECTOR operator - (VECTOR a, VECTOR b) {return {a.x-b.x,a.y-b.y};}
42 //向量*数=向量
43 VECTOR operator * (VECTOR a, ld k) {return {a.x*k,a.y*k};}
44 VECTOR operator * (ld k, VECTOR a) {return {a.x*k,a.y*k};}
45 //向量/数=向量
46 VECTOR operator / (VECTOR a, ld k) {return {a.x/k,a.y/k};}
47 //向量==向量，点==点
48 bool operator == (const VECTOR a, const VECTOR b) {return !sgn(a.x-b.x) && !sgn(a.y-b.y);}
49 //向量偏序关系（先x轴再y轴）
50 bool operator < (const VECTOR a, const VECTOR b) {return !sgn(a.x-b.x) ? a.y < b.y : a.x < b.x;}
51 //取下端点
52 POINT min(POINT p1, POINT p2) {return p1.y<p2.y ? p1:p2;}
53 //取上端点
54 POINT max(POINT p1, POINT p2) {return p1.y<p2.y ? p2:p1;}
55 //点积
56 ld dot(POINT p1, POINT p2, POINT p) {return (p.x-p1.x)*(p2.x-p1.x)+(p.y-p1.y)*(p2.y-p1.y);} //（三点式：共p1端点）
57 ld dot(VECTOR a, VECTOR b) {return a.x*b.x+a.y*b.y;} //向量式
58 //叉积
59 ld cross(POINT p1, POINT p2, POINT p) {return (p.x-p1.x)*(p2.y-p1.y)-(p.y-p1.y)*(p2.x-p1.x);} //（三点式：共p1端点）
60 ld cross(VECTOR a, VECTOR b) {return a.x*b.y-a.y*b.x;} //向量式（aXb）
61 //向量长度
62 ld vlen(VECTOR a) {return sqrt(dot(a,a));}
63 //向量旋转（逆时针）
64 VECTOR rot(VECTOR a, ld sita) {return {a.x*cos(sita)-a.y*sin(sita),a.x*sin(sita)+a.y*cos(sita)};}
65 //取单位向量
66 VECTOR unit(VECTOR a) {return a/vlen(a);}
67 //取正交单位向量
68 VECTOR norm(VECTOR a) {return VECTOR(-a.y,a.x)/vlen(a);}
69 //向量夹角余弦值
70 ld vcos(VECTOR a, VECTOR b) {return dot(a,b)/(vlen(a)*vlen(b));}
71 //向量夹角正弦值
72 ld vsin(VECTOR a, VECTOR b) {return fabs(cross(a,b))/(vlen(a)*vlen(b));}
73 //向量极角
74 //【注意】使用系统反三角函数返回nan是因为超出了定义域（大多数由于精度造成）
75 //【注意】故使用系统反三角函数需处理边界
76 ld vagl(VECTOR a, VECTOR b)
77 {
78 ld d = dot(a,b)/(vlen(a)*vlen(b));
79 if(sgn(d-1) == 0) return 0;
80 else if(sgn(d+1) == 0) return PI;
81 else return acos(d);
82 } //向量间[0,PI]
83 ld vagl(VECTOR a) {return atan2(a.y, a.x);} //单向量[-PI,PI]
84
85 //凸包
86 ll vexcnt = 0; //凸包顶点数
87 POINT xy[MAXN], xyt[MAXN]; //多边形顶点存储，临时多边形顶点存储
88 POINT convex[MAXN]; //凸包顶点数组
89
90 //求凸包周长
91 ld convexperimeter(POINT *poly, ll n)
92 {
93 ld ans = 0;
94 for(re ll i = 1; i <= n; ++i)
95 {
96 ans += vlen(poly[i]-poly[i%n+1]);
97 }
98 return ans;
99 }
100 //计算多边形面积
101 ld polygonarea(POINT *poly, ll n)
102 {
103 ld aera = 0;
104 for(re ll i = 0; i < n; ++i)
105 {
106 aera += cross(poly[i+1], poly[(i+1)%n+1]);
107 }
108 return fabs(aera)/2; //不加绝对值为有向面积
109 }
110
111 //andrew算法求凸包（凸包数组正序为逆时针）
112 //1.严格凸包（没有重复点和三点共线）2.非严格凸包（允许重复点和三点共线）
113 void andrew(ll n)
114 {
115 vexcnt = 0; //初始化数据
116 memcpy(xyt,xy,sizeof(xy)); //备份原来顶点集
117 sort(xyt+1,xyt+n+1); //排序后xyt[1]和xyt[n]一定在凸包中
118 //求解凸包顶点集
119 //正向扫描（下凸包）
120 convex[++vexcnt] = xyt[1]; //xyt[1]一定在凸包中
121 ll j = 2;
122 while(j <= n)
123 {
124 if(vexcnt <= 1) //因为xyt[2]不一定在凸包集中
125 {
126 convex[++vexcnt] = xyt[j++];
127 }
128 else
129 {
130 //取前面两个点
131 //严格凸包：sgn(cross(convex[vexcnt]-convex[vexcnt-1],xyt[j]-convex[vexcnt]))>0
132 //非严格凸包：sgn(cross(convex[vexcnt]-convex[vexcnt-1],xyt[j]-convex[vexcnt]))>=0
133 if(sgn(cross(convex[vexcnt]-convex[vexcnt-1],xyt[j]-convex[vexcnt]))>=0) convex[++vexcnt] = xyt[j++];
134 else --vexcnt;
135 }
136 }
137 //反向扫描（上凸包）
138 //至少包含了xyt[1]和xyt[n]
139 ll record = vexcnt; //记录下凸包的结果
140 ll k = n-1;
141 while(k >= 1)
142 {
143 if(vexcnt <= record)
144 {
145 convex[++vexcnt] = xyt[k--];
146 }
147 else
148 {
149 //取前面两个点
150 //严格凸包：sgn(cross(convex[vexcnt]-convex[vexcnt-1],xyt[k]-convex[vexcnt]))>0
151 //非严格凸包：sgn(cross(convex[vexcnt]-convex[vexcnt-1],xyt[k]-convex[vexcnt]))>=0
152 if(sgn(cross(convex[vexcnt]-convex[vexcnt-1],xyt[k]-convex[vexcnt]))>0) convex[++vexcnt] = xyt[k--];
153 else --vexcnt;
154 }
155 }
156 while(convex[--vexcnt]==convex[1]); //因为和xyt[1]相等的点都在下凸包中了
157 //for(re ll i = 1; i <= vexcnt; ++i) cout << "(" << convex[i].x << "," << convex[i].y << ")" << endl;
158 return;
159 }
160
161 //直线
162 struct LINE
163 {
164 POINT p; //直线上定点
165 VECTOR v; //直线方向向量
166 LINE() {}
167 //点向式
168 LINE(POINT _p, VECTOR _v):p(_p), v(_v) {}
169 //点角式
170 LINE(POINT _p, ld a):p(_p), v(VECTOR(cos(a), sin(a))) {}
171 //通过参数t唯一确定直线上点坐标（点向式：P = p + t*v）
172 POINT point(ld t) {return p + t*v;}
173 //平移向量r
174 LINE trans(VECTOR r) {return LINE(p+r,v);}
175 };
176
177 //有向直线
178 struct DIRLINE
179 {
180 POINT p; //定点
181 VECTOR v; //方向向量
182 ld a; //极角
183 DIRLINE() {}
184 DIRLINE(POINT _p, VECTOR _v):p(_p),v(_v),a(atan2(v.y,v.x)) {}
185 bool operator < (const DIRLINE &DL) const {return a < DL.a;}
186 };
187
188 //点到直线距离
189 ld PLdist(POINT p, LINE L) {return fabs(sin(vagl(p-L.p,L.v))*vlen(p-L.p));}
190 //点在有向直线左边
191 bool OnLeft(DIRLINE DL, POINT p) {return sgn(cross(DL.v,p-DL.p)) > 0;}
192 //点在有向直线右边
193 bool OnRight(DIRLINE DL, POINT p) {return sgn(cross(DL.v,p-DL.p)) < 0;}
194
195 //两点中垂线
196 //返回中垂线的个数
197 ll PerpendicularBisector(POINT p1, POINT p2, LINE &L)
198 {
199 if(p1 == p2) return 0;
200 L.p = (p1+p2)/2;
201 L.v = rot(p2-p1, PI/2);
202 return 1;
203 }
204
205 //两直线交点
206 POINT LLitesct;
207
208 //直线与直线的交点
209 //返回交点个数
210 ll LineLineIntersection(LINE L1, LINE L2, POINT &p)
211 {
212 if(sgn(cross(L1.v, L2.v)) == 0)
213 {
214 if(sgn(cross(L2.p-L1.p,L2.v)) == 0) return -1; //两直线重合
215 else return 0; //两直线平行
216 }
217 //以L1为准
218 ld t1 = -cross(L2.p-L1.p, L2.v)/cross(L2.v, L1.v);
219 p = L1.p + t1*L1.v;
220 //以L2为准
221 //ld t2 = -cross(L1.p-L2.p, L1.v)/cross(L1.v, L2.v);
222 //p = L2.p + t2*L2.v;
223 return 1;
224 }
225 ll LineLineIntersection(DIRLINE DL1, DIRLINE DL2, POINT &p)
226 {
227 if(sgn(cross(DL1.v, DL2.v)) == 0)
228 {
229 if(sgn(cross(DL2.p-DL1.p,DL2.v)) == 0) return -1; //两直线重合
230 else return 0; //两直线平行
231 }
232 //以L1为准
233 ld t1 = -cross(DL2.p-DL1.p, DL2.v)/cross(DL2.v, DL1.v);
234 p = DL1.p + t1*DL1.v;
235 //以L2为准
236 //ld t2 = -cross(L1.p-L2.p, L1.v)/cross(L1.v, L2.v);
237 //p = L2.p + t2*L2.v;
238 return 1;
239 }
240
241 //半平面交
242 //返回凸包顶点数
243 //【注意】若半平面交可能是无界区域，则需在外面加一个坐标很大的包围框
244 ll HalfPlaneIntersection(DIRLINE *DL, ll n, POINT *poly)
245 {
246 sort(DL+1,DL+n+1);
247 ll first, last; //双端队列的首尾元素下标
248 DIRLINE *q = new DIRLINE[n+1]; //双端队列
249 POINT *p = new POINT[n+1]; //p[i]为q[i]和q[i+1]的交点
250 first = last = 1;
251 q[1] = DL[1]; //初始化只有一个半平面L[0]
252 for(re ll i = 2; i <= n; ++i)
253 {
254 while(first < last && !OnLeft(DL[i], p[last-1])) --last;
255 while(first < last && !OnLeft(DL[i], p[first])) ++first;
256 q[++last] = DL[i];
257 if(sgn(cross(q[last].v, q[last-1].v)) == 0) //两向量平行
258 {
259 --last;
260 if(OnLeft(q[last], DL[i].p)) q[last] = DL[i]; //同向取内测的一个
261 }
262 if(first < last) LineLineIntersection(q[last-1],q[last],p[last-1]);
263 }
264 while(first < last && !OnLeft(q[first],p[last-1])) --last;
265 //删除无用平面
266 if(last-first <= 1) return 0; //空集
267 LineLineIntersection(q[last],q[first],p[last]); //计算首尾两个半平面交的交点
268 ll cnt = 0;
269 for(re ll i = first; i <= last; ++i) poly[++cnt] = p[i];
270 return cnt;
271 }
272
273 int main()
274 {
275 //UVALive-3890
276 ios::sync_with_stdio(false);
277 ll n;
278 //scanf("%lld", &n);
279 while(cin >> n)
280 //while(~scanf("%d", &n))
281 {
282 if(n == 0) break;
283 POINT *p = new POINT[n+1];
284 for(re ll i = 1; i <= n; ++i)
285 {
286 cin >> p[i].x >> p[i].y;
287 //scanf("%lf%lf%lf%lf", &p1.x, &p1.y, &p2.x, &p2.y);
288 }
289 VECTOR *v = new VECTOR[n+1];
290 VECTOR *r = new VECTOR[n+1];
291 for(re ll i = 1; i <= n; ++i)
292 {
293 v[i] = p[i%n+1]-p[i];
294 r[i] = norm(v[i]); //单位正交向量
295 }
296 //二分答案
297 ld left = 0, right = 20000;
298 DIRLINE *DL = new DIRLINE[n+1];
299 POINT *poly = new POINT[n+1];
300 while(right-left > 1e-6)
301 {
302 ld mid = (left+right)/2;
303 for(re ll i = 1; i <= n; ++i)
304 {
305 DL[i] = DIRLINE(p[i]+r[i]*mid,v[i]);
306 }
307 ll cnt = HalfPlaneIntersection(DL,n,poly);
308 if(!cnt) right = mid;
309 else left = mid;
310 }
311 cout << fixed << setprecision(6) << left << endl;
312 }
313 return 0;
314 }
原文地址：https://www.cnblogs.com/BlueHeart0621/p/12238221.html
时间： 2024-10-05 03:09:42