计算几何基础
本文大部分内容来自b站浙师大计算几何基础以及kuangbin的模板，再次鸣谢
1.算法原理
1.1 基本符号与公式
1.1.1 基础符号表示
1.1.1.1 PI的表示
double const PI = acos(double(-1));
1.1.2 基本公式
1.1.2.1 正弦定理
在三角形 △ ABC \triangle \text{ABC} △ABC 中，若角 A , B , C A,B,C A,B,C 所对边分别为 a , b , c a,b,c a,b,c，则有：
$ \frac{a}{\sin A}=\frac{b}{\sin B}=\frac{c}{\sin C}=2R $
其中， R R R 为 △ ABC \triangle \text{ABC} △ABC 的外接圆半径
1.1.2.2 余弦定理
在三角形 △ ABC \triangle \text{ABC} △ABC 中，若角 A , B , C A,B,C A,B,C 所对边分别为 a , b , c a,b,c a,b,c ，则有：
$ \begin{aligned} a2&=b2+c^2-2bc\cos A\ b2&=a2+c^2-2ac\cos B\ c2&=a2+b^2-2ab\cos C \end{aligned} $
1.1.2.3 角度与弧度转化
弧度->角度： 弧度 * (180.0 / PI)
角度->弧度： 角度 / (180.0 / PI)
1.2 向量
设向量 v 1 = ( x 1 , y 1 ) , v 2 = ( x 2 , y 2 ) v_1=(x_1,y_1),v_2=(x_2,y_2) v1=(x1,y1),v2=(x2,y2) ，定义如下运算
1.2.1 向量加法
v 1 + v 2 = ( x 1 + x 2 , y 1 + y 2 ) v_1+v_2=(x_1+x_2,y_1+y_2) v1+v2=(x1+x2,y1+y2)
1.2.2 向量减法
v 1 − v 2 = ( x 1 − x 2 , y 1 − y 2 ) v_1-v_2=(x_1-x_2,y_1-y_2) v1−v2=(x1−x2,y1−y2)
若 P = ( x 1 , y 1 ) , Q = ( x 2 , y 2 ) P=(x_1,y_1),Q=(x_2,y_2) P=(x1,y1),Q=(x2,y2) ，则 P Q → = Q − P = ( x 2 − x 1 , y 2 − y 1 ) \overrightarrow{PQ}=Q-P=(x_2-x_1,y_2-y_1) PQ=Q−P=(x2−x1,y2−y1)
1.2.3 向量模长
∣ v 1 ∣ = x 1 2 + y 1 2 |v_1|=\sqrt{x_1^2+y_1^2} ∣v1∣=x12+y12
向量模长可以用来求两点间的距离
1.2.4 向量数乘
a v 1 = ( a x 1 , a y 1 ) , a ∈ R av_1=(ax_1,ay_1),a\in \mathbb{R} av1=(ax1,ay1),a∈R
向量数乘可以实现向量的长度伸缩
1.2.5 向量内积（点积）
v 1 ⋅ v 2 = ∣ v 1 ∣ ∣ v 2 ∣ cos < v 1 , v 2 > = x 1 x 2 + y 1 y 2 v_1\cdot v_2=|v_1||v_2|\cos<v_1,v_2>=x_1x_2+y_1y_2 v1⋅v2=∣v1∣∣v2∣cos<v1,v2>=x1x2+y1y2
v 1 ⋅ v 2 = 0 v_1\cdot v_2=0 v1⋅v2=0 当且仅当 v 1 ⊥ v 2 v_1\perp v_2 v1⊥v2
向量内积可以用来求向量间的夹角
1.2.6 向量外积（叉积）
v 1 × v 2 = ∣ x 1 y 1 x 2 y 2 ∣ = x 1 y 2 − x 2 y 1 v_1\times v_2=\begin{vmatrix} x_1 & y_1\\ x_2 & y_2 \end{vmatrix}=x_1y_2-x_2y_1 v1×v2=∣∣∣∣x1x2y1y2∣∣∣∣=x1y2−x2y1
∣ v 1 × v 2 ∣ = ∣ v 1 ∣ ∣ v 2 ∣ sin < v 1 , v 2 > |v_1\times v_2|=|v_1||v_2|\sin <v_1,v_2> ∣v1×v2∣=∣v1∣∣v2∣sin<v1,v2>
外积可以用来求面积，以 v 1 , v 2 v_1,v_2 v1,v2 为邻边的平行四边形面积为 ∣ v 1 × v 2 ∣ |v_1\times v_2| ∣v1×v2∣
v 1 × v 2 = 0 v_1\times v_2=0 v1×v2=0 当且仅当 v 1 ∥ v 2 v_1\parallel v_2 v1∥v2
外积可以用来判断向量间的位置关系，若 v 1 v_1 v1 旋转到 v 2 v_2 v2 的方向为顺时针，则 v 1 × v 2 < 0 v_1\times v_2<0 v1×v2<0 ，反之 v 1 × v 2 > 0 v_1\times v_2>0 v1×v2>0
1.2.7 向量旋转
向量
v
1
v_1
v1 逆时针旋转
θ
\theta
θ 后的坐标满足
{
x
′
=
x
1
cos
θ
−
y
1
sin
θ
y
′
=
x
1
sin
θ
+
y
1
cos
θ
\begin{cases} x'=x_1\cos \theta-y_1\sin \theta\\ y'=x_1\sin \theta+y_1\cos \theta \end{cases}
{x′=x1cosθ−y1sinθy′=x1sinθ+y1cosθ
1.3 直线与线段
1.3.1 线段相交问题
线段
A
B
AB
AB 与
C
D
CD
CD 相交（不考虑端点）的充分必要条件是
(
C
A
→
⋅
C
B
→
)
(
D
A
→
⋅
D
B
→
)
<
0
,
(
A
C
→
⋅
A
D
→
)
(
B
C
→
⋅
B
D
→
)
<
0
(\overrightarrow{CA}\cdot \overrightarrow{CB}) (\overrightarrow{DA}\cdot \overrightarrow{DB})<0, (\overrightarrow{AC}\cdot \overrightarrow{AD}) (\overrightarrow{BC}\cdot \overrightarrow{BD})<0
(CA⋅CB)(DA⋅DB)<0,(AC⋅AD)(BC⋅BD)<0
1.3.2 点到直线的距离
要计算点A到直线MN的距离，可以构建以AM，MN为邻边的平行四边形，其面积
S
=
∣
M
A
→
×
M
N
→
∣
S=|\overrightarrow{MA}\times \overrightarrow{MN}|
S=∣MA×MN∣
平行四边形的面积为底乘高，选取MN为底，高为
d
=
S
∣
M
N
→
∣
d=\frac{S}{\left|\overrightarrow{MN}\right|}
d=∣∣∣MN∣∣∣S
即为所求的A到直线MN的距离
1.3.3 两直线交点
在实际应用中，通常的已知量是直线上某一点的坐标和直线的方向向量，对于两直线
l
1
l_{1}
l1,
l
2
\ l_{2}
l2 ,设
P
(
x
1
,
y
1
)
P\left( x_{1},y_{1} \right)
P(x1,y1) ,
Q
(
x
2
,
y
2
)
\text{Q}\left( x_{2},y_{2} \right)
Q(x2,y2) 分别在
l
1
l_{1}
l1 ,
l
2
\ l_{2}
l2 上，
l
1
l_{1}
l1 ,
l
2
\ l_{2}
l2 的方向向量分别为
v
=
(
a
1
,
b
1
)
v = \left( a_{1},b_{1} \right)
v=(a1,b1) ,
w
=
(
a
2
,
b
2
)
w = \left( a_{2},b_{2} \right)
w=(a2,b2) ,由此可以得到两直线的方程
l
1
:
(
x
−
x
1
,
y
−
y
1
)
×
(
a
1
,
b
1
)
=
0
l_{1}:\left( x - x_{1},y - y_{1} \right) \times \left( a_{1},b_{1} \right) = 0
l1:(x−x1,y−y1)×(a1,b1)=0
l 2 : ( x − x 2 , y − y 2 ) × ( a 2 , b 2 ) = 0 l_{2}:\left( x - x_{2},y - y_{2} \right) \times \left( a_{2},b_{2} \right) = 0 l2:(x−x2,y−y2)×(a2,b2)=0
即
l
1
:
a
1
x
−
b
1
y
=
a
1
x
1
−
b
1
y
1
l_{1}:a_{1}x - b_{1}y = a_{1}x_{1} - b_{1}y_{1}
l1:a1x−b1y=a1x1−b1y1
l 2 : a 2 x − b 2 y = a 2 x 2 − b 2 y 2 l_{2}:a_{2}x - b_{2}y = a_{2}x_{2} - b_{2}y_{2} l2:a2x−b2y=a2x2−b2y2
联立两直线的方程，由克拉默法则得，方程组的解为
KaTeX parse error: Undefined control sequence: \ at position 467: …matrix} \right.\̲ ̲
进一步进行化简，得到
(
x
,
y
)
=
P
+
v
⋅
w
×
u
v
×
w
(x,y)=P+v\cdot \frac{w\times u}{v\times w}
(x,y)=P+v⋅v×ww×u
其中
u
=
−
P
Q
→
u=-\overrightarrow{PQ}
u=−PQ
1.4 多边形
1.4.1 点和多边形的位置关系
设有（凸）
n
(
n
≥
3
)
n(n≥3)
n(n≥3) 边形
P
0
P
2
…
P
n
−
1
P_0 P_2\dots P_{n-1}
P0P2…Pn−1，点的顺序为顺时针或逆时针，以及点A，记
θ
i
=
{
<
A
P
i
→
,
A
P
i
+
1
→
>
,
i
<
n
−
1
<
A
P
n
−
1
→
,
A
P
0
→
>
,
i
=
n
−
1
\theta_{i} = \left\{ \begin{matrix} < \overrightarrow{AP_{i}},\overrightarrow{AP_{i + 1}} > ,i < n - 1 \\ < \overrightarrow{AP_{n - 1}},\overrightarrow{AP_{0}} > ,i = n - 1 \\ \end{matrix} \right.\
θi={<APi,APi+1>,i<n−1<APn−1,AP0>,i=n−1
点在多边形内等价于
∑
i
=
0
n
−
1
θ
i
=
2
π
\sum_{i = 0}^{n - 1}\theta_{i} = 2\pi
i=0∑n−1θi=2π
1.4.2 多边形的面积
设有（凸） n ( n ≥ 3 ) n(n≥3) n(n≥3) 边形 P 0 P 2 … P n − 1 P_0 P_2\dots P_{n-1} P0P2…Pn−1 ，点的顺序为顺时针或逆时针，以及多边形内一点A，把多边形切割成如下所示n个三角形
多边形的面积等于所有三角形（有向）面积之和，代入坐标
P
i
(
x
i
,
y
i
)
,
i
=
0
,
1
,
…
,
n
−
1
P_i (x_i,y_i ),i=0,1,\dots,n-1
Pi(xi,yi),i=0,1,…,n−1 计算得
S
=
∣
1
2
∑
i
=
0
n
−
2
(
x
i
y
i
+
1
−
x
i
+
1
y
i
)
+
1
2
(
x
n
−
1
y
0
−
x
0
y
n
−
1
)
∣
S = \left| \frac{1}{2}\sum_{i = 0}^{n - 2}\left( x_{i}y_{i + 1} - x_{i + 1}y_{i} \right) + \frac{1}{2}\left( x_{n - 1}y_{0} - x_{0}y_{n - 1} \right) \right|
S=∣∣∣∣∣21i=0∑n−2(xiyi+1−xi+1yi)+21(xn−1y0−x0yn−1)∣∣∣∣∣
与A的坐标无关，因此A可任取，甚至可取在多边形外，通常为计算方便，取A为坐标原点
1.4.3 凸包
在一个实向量空间
V
V
V 中，对于给定集合
X
X
X ，所有包含
X
X
X 的凸集的交集
S
S
S 称为
X
X
X 的凸包
S
=
∩
X
⊂
K
⊂
V
,
K
is convex
K
S=\cap_{X\subset K\subset V,K\text{ is convex}}K
S=∩X⊂K⊂V,K is convexK
1.4.3.1 Graham’s scan算法
第一步：找到最下边的点，如果有多个点纵坐标相同的点都在最下方，则选取最左边的，记为点A。这一步只需要扫描一遍所有的点即可，时间复杂度为 O ( n ) O(n) O(n)
第二步：将所有的点按照 A P i AP_i APi 的极角大小进行排序，极角相同的按照到点A的距离排序。时间复杂度为 O ( n l o g n ) O(nlogn) O(nlogn)
第三步：维护一个栈，以保存当前的凸包。按第二步中排序得到的结果，依次将点加入到栈中，如果当前点与栈顶的两个点不是“向左转”的，就表明当前栈顶的点并不在凸包上，而我们需要将其弹出栈，重复这一个过程直到当前点与栈顶的两个点是“向左转”的。这一步的时间复杂度为 O ( n ) O(n) O(n)
1.4.3.2 Andrew’s monotone chain 算法
原理与Graham’s scan算法相似，但上下凸包是分开维护的
1.5 圆
1.5.1 圆的参数方程
以
(
x
0
,
y
0
)
(x_{0},y_{0})
(x0,y0)为圆心，
r
r
r为半径的圆的参数方程为
KaTeX parse error: Undefined control sequence: \ at position 99: …matrix} \right.\̲ ̲
根据圆上一点和圆心连线与
x
x
x轴正向的夹角可求得该点的坐标
1.5.2 两圆交点
设两圆 C 1 , C 2 C_{1},C_{2} C1,C2，其半径为 r 1 , r 2 ( r 1 ≥ r 2 ) r_{1},r_{2}(r_{1} \geq r_{2}) r1,r2(r1≥r2)，圆心距为 d d d，则有
①两圆重合 ⟺ d = 0 r 1 = r 2 \Longleftrightarrow d = 0\ \ r_{1} = r_{2} ⟺d=0 r1=r2
②两圆外离 ⟺ d > r 1 + r 2 \Longleftrightarrow d > r_{1} + r_{2} ⟺d>r1+r2
③两圆外切 ⟺ d = r 1 + r 2 \Longleftrightarrow d = r_{1} + r_{2} ⟺d=r1+r2
④两圆相交 ⟺ r 1 − r 2 < d < r 1 + r 2 \Longleftrightarrow r_{1} - r_{2} < d < r_{1} + r_{2} ⟺r1−r2<d<r1+r2
⑤两圆内切 ⟺ d = r 1 − r 2 \Longleftrightarrow d = r_{1} - r_{2} ⟺d=r1−r2
⑥两圆内含 ⟺ d < r 1 − r 2 \Longleftrightarrow d < r_{1} - r_{2} ⟺d<r1−r2
对于情形④，如下图所示，要求A与B的坐标，只需求
∠
A
C
1
D
\angle AC_{1}D
∠AC1D与
∠
B
C
1
D
\angle BC_{1}D
∠BC1D，进而通过圆的参数方程即可求得
∠
A
C
1
D
=
∠
C
2
C
1
D
+
∠
A
C
1
C
2
\angle AC_{1}D = \angle C_{2}C_{1}D + \angle AC_{1}C_{2}
∠AC1D=∠C2C1D+∠AC1C2
∠ B C 1 D = ∠ C 2 C 1 D − ∠ A C 1 C 2 \angle BC_{1}D = \angle C_{2}C_{1}D - \angle AC_{1}C_{2} ∠BC1D=∠C2C1D−∠AC1C2
∠ C 2 C 1 D \angle C_{2}C_{1}D ∠C2C1D可以通过 C 1 , C 2 C_{1},C_{2} C1,C2的坐标求得，而 ∠ A C 1 C 2 \angle AC_{1}C_{2} ∠AC1C2可以通过 Δ A C 1 C 2 \Delta AC_{1}C_{2} ΔAC1C2上的余弦定理求得
对于情形③和情形⑤，上述方法求得的两点坐标是相同的，即为切点的坐标
2.模板
2.1 二维计算几何
#include <bits/stdc++.h>
using namespace std;
// 计算几何模板
const double eps = 1e-8;
const double inf = 1e20;
const double pi = acos(-1.0);
const int maxp = 1010;
// 和0做比较
int sgn(double x) {
if (fabs(x) < eps) return 0; // =0
if (x < 0)
return -1; // < 0
else
return 1; // > 0
}
// 计算x的平方
inline double sqr(double x) { return x * x; }
/*
* Point
* Point() - 空构造
* Point(double _x,double _y) - 普通构造
* input() - 浮点输入
* output() - %.2f 输出
* operator == - 判断是否相等
* operator < - 判断是否更小，先比较x，然后比较y
* operator - - 返回点p减去一个点的
* operator ^ - 叉乘
* operator * - 点乘
* len() - 得到向量p的长度
* len2() - 得到向量p长度的平方
* distance(Point p) - 得到点p和其他点之间的距离
* operator + Point b - 返回向量加的结果
* operator * double k - 把点p的x和y都放大k倍
* operator / double k - 把点p的x和y都缩小k倍
* rad(Point a,Point b)- 得到点p去看点a和点b的夹角的弧度
* trunc(double r) - 把向量p进行放缩，使得其长度为r
* rotleft() - 把向量p逆时针旋转90度
* rotright() - 把向量p顺时针旋转90度
* rotate(Point p,double angle) - 绕着点p逆时针旋转90度
*/
struct Point {
double x, y;
Point() {}
Point(double _x, double _y) {
x = _x;
y = _y;
}
void input() { scanf("%lf%lf", &x, &y); }
void output() { printf("%.2f %.2f\n", x, y); }
bool operator==(Point b) const {
return sgn(x - b.x) == 0 && sgn(y - b.y) == 0;
}
bool operator<(Point b) const {
return sgn(x - b.x) == 0 ? sgn(y - b.y) < 0 : x < b.x;
}
Point operator-(const Point &b) const { return Point(x - b.x, y - b.y); }
//叉积
double operator^(const Point &b) const { return x * b.y - y * b.x; }
//点积
double operator*(const Point &b) const { return x * b.x + y * b.y; }
//返回向量长度
double len() {
// hypot(x, y), 即sqrt(x * x + y * y)
return hypot(x, y); //库函数
}
//返回长度的平方
double len2() { return x * x + y * y; }
//返回两点的距离
double distance(Point p) { return hypot(x - p.x, y - p.y); }
Point operator+(const Point &b) const { return Point(x + b.x, y + b.y); }
Point operator*(const double &k) const { return Point(x * k, y * k); }
Point operator/(const double &k) const { return Point(x / k, y / k); }
//计算pa和pb的夹角， 就是求这个点看a,b 所成的夹角的弧度
double rad(Point a, Point b) {
Point p = *this;
return fabs(atan2(fabs((a - p) ^ (b - p)), (a - p) * (b - p)));
}
//把点p进行放缩，使得点p化为长度为r的向量
Point trunc(double r) {
double l = len();
if (!sgn(l)) return *this;
r /= l;
return Point(x * r, y * r);
}
//逆时针旋转90度
Point rotleft() { return Point(-y, x); }
//顺时针旋转90度
Point rotright() { return Point(y, -x); }
//绕着p点逆时针旋转angle
Point rotate(Point p, double angle) {
Point v = (*this) - p;
double c = cos(angle), s = sin(angle);
return Point(p.x + v.x * c - v.y * s, p.y + v.x * s + v.y * c);
}
};
/*
* Stores two points
* Line() - 空构造
* Line(Point _s,Point _e) - 两点确定一条线段
* operator == - 判断两条线段是否相同
* Line(Point p,double angle) - 一个点和一个角度确定一条直线
* Line(double a,double b,double c) - 直线方程确定一条直线：ax+by+c=0
* input() - 输入线段的2个端点
* adjust() - 调整使得线段的起点s < 终点e
* length() - 求线段的长度
* angle() - 以弧度形式返回一个直线的倾斜角 0 <= angle <
* pi relation(Point p) - 返回点和直线的关系
* - 3 点在直线上
* - 1 点在直线左侧
* - 2 点在直线右侧
*
* pointonseg(double p) - 判断点p是否在线段上
* parallel(Line v) - 判断直线v是否和直线p平行
* segcrossseg(Line v) - 判断线段是否相交
* returns 0 线段不交
* returns 1 非规范相交：交点在端点
* returns 2 规范相交
* linecrossseg(Line v) - 判断直线和线段是否相交
* linecrossline(Line v) - 判断直线和直线是否平行
* 0 if parallel
* 1 if coincides
* 2 if intersects
* crosspoint(Line v) - 返回两条直线的交点（要先保证相交）
* dispointtoline(Point p) - 返回p点到直线的距离
* dispointtoseg(Point p) - 返回p点到线段的距离
* dissegtoseg(Line v) - 返回线段到线段的距离
* lineprog(Point p) - 返回p点在直线上的投影
* symmetrypoint(Point p) -
* 返回p点在直线上的反射（p点绕直线180度旋转后的点）
*
*/
struct Line {
Point s, e;
Line() {}
// 两点确定一条线段
Line(Point _s, Point _e) {
s = _s;
e = _e;
}
// 判断直线是否相等
bool operator==(Line v) { return (s == v.s) && (e == v.e); }
//一个点和倾斜角angle确定直线,0<=angle<pi
Line(Point p, double angle) {
s = p;
if (sgn(angle - pi / 2) == 0) {
e = (s + Point(0, 1));
} else {
e = (s + Point(1, tan(angle)));
}
}
//直线ax+by+c=0
Line(double a, double b, double c) {
if (sgn(a) == 0) {
s = Point(0, -c / b);
e = Point(1, -c / b);
} else if (sgn(b) == 0) {
s = Point(-c / a, 0);
e = Point(-c / a, 1);
} else {
s = Point(0, -c / b);
e = Point(1, (-c - a) / b);
}
}
void input() {
s.input();
e.input();
}
// 保证直线的两个端点满足：s < e
void adjust() {
if (e < s) swap(s, e);
}
// 求线段长度
double length() { return s.distance(e); }
//返回直线倾斜角 0<=angle<pi
double angle() {
double k = atan2(e.y - s.y, e.x - s.x);
if (sgn(k) < 0) k += pi;
if (sgn(k - pi) == 0) k -= pi;
return k;
}
//返回点和直线关系
// 1 在左侧; 2 在右侧; 3 在直线上
int relation(Point p) {
int c = sgn((p - s) ^ (e - s));
if (c < 0)
return 1;
else if (c > 0)
return 2;
else
return 3;
}
// 点在线段上的判断
bool pointonseg(Point p) {
return sgn((p - s) ^ (e - s)) == 0 && sgn((p - s) * (p - e)) <= 0;
}
// 两向量平行(对应直线平行或重合)
bool parallel(Line v) { return sgn((e - s) ^ (v.e - v.s)) == 0; }
//两线段相交判断:规范相交：交点不在端点
// 2 规范相交；1 非规范相交；0 不相交
int segcrossseg(Line v) {
int d1 = sgn((e - s) ^ (v.s - s));
int d2 = sgn((e - s) ^ (v.e - s));
int d3 = sgn((v.e - v.s) ^ (s - v.s));
int d4 = sgn((v.e - v.s) ^ (e - v.s));
if ((d1 ^ d2) == -2 && (d3 ^ d4) == -2) return 2;
return (d1 == 0 && sgn((v.s - s) * (v.s - e)) <= 0) ||
(d2 == 0 && sgn((v.e - s) * (v.e - e)) <= 0) ||
(d3 == 0 && sgn((s - v.s) * (s - v.e)) <= 0) ||
(d4 == 0 && sgn((e - v.s) * (e - v.e)) <= 0);
}
//直线和线段相交判断
//-*this line -v seg
// 2 规范相交
// 1 非规范相交
// 0 不相交
int linecrossseg(Line v) {
int d1 = sgn((e - s) ^ (v.s - s));
int d2 = sgn((e - s) ^ (v.e - s));
if ((d1 ^ d2) == -2) return 2;
return (d1 == 0 || d2 == 0);
}
//两直线关系
// 0 平行；1 重合；2 相交
int linecrossline(Line v) {
if ((*this).parallel(v)) return v.relation(s) == 3;
return 2;
}
//求两直线的交点
//要保证两直线不平行或重合
Point crosspoint(Line v) {
double a1 = (v.e - v.s) ^ (s - v.s);
double a2 = (v.e - v.s) ^ (e - v.s);
return Point((s.x * a2 - e.x * a1) / (a2 - a1),
(s.y * a2 - e.y * a1) / (a2 - a1));
}
//点到直线的距离
double dispointtoline(Point p) {
return fabs((p - s) ^ (e - s)) / length();
}
//点到线段的距离
double dispointtoseg(Point p) {
if (sgn((p - s) * (e - s)) < 0 || sgn((p - e) * (s - e)) < 0)
return min(p.distance(s), p.distance(e));
return dispointtoline(p);
}
//返回线段到线段的距离
//前提是两线段不相交，相交距离就是0了
double dissegtoseg(Line v) {
return min(min(dispointtoseg(v.s), dispointtoseg(v.e)),
min(v.dispointtoseg(s), v.dispointtoseg(e)));
}
//返回点p在直线上的投影
Point lineprog(Point p) {
return s + (((e - s) * ((e - s) * (p - s))) / ((e - s).len2()));
}
//返回点p关于直线的对称点（求反射）
Point symmetrypoint(Point p) {
Point q = lineprog(p);
return Point(2 * q.x - p.x, 2 * q.y - p.y);
}
};
//圆
/*
* 存储圆心和半径
* circle() - 空构造
* circle(Point _p, double _r) - 圆心和半径构造一个圆
* circle(double x, double y, double _r) -圆心(x, y)和半径r构造一个圆
* circle(Point a, Point b, Point c)) - 求三角形a、b、c确定外接圆
* circle(Point a, Point b, Point c, bool t)) - 求三角形的内接圆
* input() - 圆心+半径
* output() - 输出圆：圆心+半径
* operator==(circle v) - 判断和圆v是否相同
* operator<(circle v) -
* 比较2个圆的关系，先比较圆心，然后比较圆的半径 area() - 返回圆的面积
* circumference() - 返回圆的周长
* relation(Point b) - 点和圆的关系：0 圆外；1 圆上；2 圆内
* relationseg(Line v) - 线段和圆的关系：2 相交；1 相切；0 不交
* relationline(Line v)) - 直线和圆的关系：2 相交；1 相切；0 不交
* relationcircle(circle v) - 两圆的关系：5 相离；4 外切；3 相交；2
* 内切；1 内含 pointcrosscircle(circle v, Point &p1, Point &p2) -
* 求两个圆的交点，返回0表示没有交点，返回1是一个交点，2是两个交点,
* 交点存储在p1和p2 pointcrossline(Line v, Point &p1, Point &p2) -
* 求直线和圆的交点，返回交点个数，交点存储在p1和p2 返回p点到线段的距离
* gercircle(Point a, Point b, double r1, circle &c1, circle &c2) -
* 得到过a,b两点，半径为r1的两个圆，返回t为交点数，两个圆分别为r1和r2
* getcircle(Line u, Point q, double r1, circle &c1, circle &c2) -
* 得到与直线u相切，过点q,半径为r1的圆 getcircle(Line u, Line v, double r1,
* circle &c1, circle &c2, circle &c3, circle &c4) -
* 同时与直线u,v相切，半径为r1的圆 getcircle(circle cx, circle cy, double r1,
* circle &c1, circle &c2) - 同时与不相交圆cx,cy相切，半径为r1的圆
* tangentline(Point q, Line &u, Line &v) - 过一点作圆的切线(先判断点和圆的关系)
* areacircle(circle v) - 求两圆相交的面积
* areatriangle(Point a, Point b) - 求圆和三角形pab的相交面积：p为圆心
*/
struct circle {
Point p; //圆心
double r; //半径
circle() {}
// 圆的构造
circle(Point _p, double _r) {
p = _p;
r = _r;
}
circle(double x, double y, double _r) {
p = Point(x, y);
r = _r;
}
// 三角形的外接圆
// 需要Point的+ / rotate() 以及Line的crosspoint()
// 利用两条边的中垂线得到圆心
circle(Point a, Point b, Point c) {
Line u = Line((a + b) / 2, ((a + b) / 2) + ((b - a).rotleft()));
Line v = Line((b + c) / 2, ((b + c) / 2) + ((c - b).rotleft()));
p = u.crosspoint(v);
r = p.distance(a);
}
// 三角形的内切圆
// 参数bool t没有作用，只是为了和上面外接圆函数区别
circle(Point a, Point b, Point c, bool t) {
Line u, v;
double m = atan2(b.y - a.y, b.x - a.x), n = atan2(c.y - a.y, c.x - a.x);
u.s = a;
u.e = u.s + Point(cos((n + m) / 2), sin((n + m) / 2));
v.s = b;
m = atan2(a.y - b.y, a.x - b.x), n = atan2(c.y - b.y, c.x - b.x);
v.e = v.s + Point(cos((n + m) / 2), sin((n + m) / 2));
p = u.crosspoint(v);
r = Line(a, b).dispointtoseg(p);
}
//输入: 圆心+半径
void input() {
p.input();
scanf("%lf", &r);
}
//输出：圆心的x、圆心的y、圆的半径
void output() { printf("%.2lf %.2lf %.2lf\n", p.x, p.y, r); }
// 判断两个圆是否相同
bool operator==(circle v) { return (p == v.p) && sgn(r - v.r) == 0; }
// 比较2个圆的关系，先比较圆心，然后比较圆的半径
bool operator<(circle v) const {
return ((p < v.p) || ((p == v.p) && sgn(r - v.r) < 0));
}
// 返回圆的面积
double area() { return pi * r * r; }
// 返回圆的周长
double circumference() { return 2 * pi * r; }
//点和圆的关系
// 0 圆外
// 1 圆上
// 2 圆内
int relation(Point b) {
double dst = b.distance(p);
if (sgn(dst - r) < 0)
return 2;
else if (sgn(dst - r) == 0)
return 1;
return 0;
}
//线段和圆的关系
//比较的是圆心到线段的距离和半径的关系
// 2 相交；1 相切；0 不交
int relationseg(Line v) {
double dst = v.dispointtoseg(p);
if (sgn(dst - r) < 0)
return 2;
else if (sgn(dst - r) == 0)
return 1;
return 0;
}
//直线和圆的关系
//比较的是圆心到直线的距离和半径的关系
// 2 相交；1 相切；0 不交
int relationline(Line v) {
double dst = v.dispointtoline(p);
if (sgn(dst - r) < 0)
return 2;
else if (sgn(dst - r) == 0)
return 1;
return 0;
}
//两圆的关系
// 5 相离；4 外切；3 相交；2 内切；1 内含
int relationcircle(circle v) {
double d = p.distance(v.p);
if (sgn(d - r - v.r) > 0) return 5;
if (sgn(d - r - v.r) == 0) return 4;
double l = fabs(r - v.r);
if (sgn(d - r - v.r) < 0 && sgn(d - l) > 0) return 3;
if (sgn(d - l) == 0) return 2;
if (sgn(d - l) < 0) return 1;
}
//求两个圆的交点，返回0表示没有交点，返回1是一个交点，2是两个交点,
//交点存储在p1和p2 需要relationcircle()
int pointcrosscircle(circle v, Point &p1, Point &p2) {
int rel = relationcircle(v);
if (rel == 1 || rel == 5) return 0;
double d = p.distance(v.p);
double l = (d * d + r * r - v.r * v.r) / (2 * d);
double h = sqrt(r * r - l * l);
Point tmp = p + (v.p - p).trunc(l);
p1 = tmp + ((v.p - p).rotleft().trunc(h));
p2 = tmp + ((v.p - p).rotright().trunc(h));
if (rel == 2 || rel == 4) return 1;
return 2;
}
// 求直线和圆的交点，返回交点个数，交点存储在p1和p2
int pointcrossline(Line v, Point &p1, Point &p2) {
if (!(*this).relationline(v)) return 0;
Point a = v.lineprog(p);
double d = v.dispointtoline(p);
d = sqrt(r * r - d * d);
if (sgn(d) == 0) {
p1 = a;
p2 = a;
return 1;
}
p1 = a + (v.e - v.s).trunc(d);
p2 = a - (v.e - v.s).trunc(d);
return 2;
}
//得到过a,b两点，半径为r1的两个圆，返回t为交点数，两个圆分别为r1和r2
int gercircle(Point a, Point b, double r1, circle &c1, circle &c2) {
circle x(a, r1), y(b, r1);
int t = x.pointcrosscircle(y, c1.p, c2.p);
if (!t) return 0;
c1.r = c2.r = r;
return t;
}
//得到与直线u相切，过点q,半径为r1的圆
int getcircle(Line u, Point q, double r1, circle &c1, circle &c2) {
double dis = u.dispointtoline(q);
if (sgn(dis - r1 * 2) > 0) return 0;
if (sgn(dis) == 0) {
c1.p = q + ((u.e - u.s).rotleft().trunc(r1));
c2.p = q + ((u.e - u.s).rotright().trunc(r1));
c1.r = c2.r = r1;
return 2;
}
Line u1 = Line((u.s + (u.e - u.s).rotleft().trunc(r1)),
(u.e + (u.e - u.s).rotleft().trunc(r1)));
Line u2 = Line((u.s + (u.e - u.s).rotright().trunc(r1)),
(u.e + (u.e - u.s).rotright().trunc(r1)));
circle cc = circle(q, r1);
Point p1, p2;
if (!cc.pointcrossline(u1, p1, p2)) cc.pointcrossline(u2, p1, p2);
c1 = circle(p1, r1);
if (p1 == p2) {
c2 = c1;
return 1;
}
c2 = circle(p2, r1);
return 2;
}
// 同时与直线u,v相切，半径为r1的圆
int getcircle(Line u, Line v, double r1, circle &c1, circle &c2, circle &c3,
circle &c4) {
if (u.parallel(v)) return 0; //两直线平行
Line u1 = Line(u.s + (u.e - u.s).rotleft().trunc(r1),
u.e + (u.e - u.s).rotleft().trunc(r1));
Line u2 = Line(u.s + (u.e - u.s).rotright().trunc(r1),
u.e + (u.e - u.s).rotright().trunc(r1));
Line v1 = Line(v.s + (v.e - v.s).rotleft().trunc(r1),
v.e + (v.e - v.s).rotleft().trunc(r1));
Line v2 = Line(v.s + (v.e - v.s).rotright().trunc(r1),
v.e + (v.e - v.s).rotright().trunc(r1));
c1.r = c2.r = c3.r = c4.r = r1;
c1.p = u1.crosspoint(v1);
c2.p = u1.crosspoint(v2);
c3.p = u2.crosspoint(v1);
c4.p = u2.crosspoint(v2);
return 4;
}
// 同时与不相交圆cx,cy相切，半径为r1的圆
int getcircle(circle cx, circle cy, double r1, circle &c1, circle &c2) {
circle x(cx.p, r1 + cx.r), y(cy.p, r1 + cy.r);
int t = x.pointcrosscircle(y, c1.p, c2.p);
if (!t) return 0;
c1.r = c2.r = r1;
return t;
}
// 过一点作圆的切线(先判断点和圆的关系)
int tangentline(Point q, Line &u, Line &v) {
int x = relation(q);
if (x == 2) return 0;
if (x == 1) {
u = Line(q, q + (q - p).rotleft());
v = u;
return 1;
}
double d = p.distance(q);
double l = r * r / d;
double h = sqrt(r * r - l * l);
u = Line(q, p + ((q - p).trunc(l) + (q - p).rotleft().trunc(h)));
v = Line(q, p + ((q - p).trunc(l) + (q - p).rotright().trunc(h)));
return 2;
}
// 求两圆相交的面积
double areacircle(circle v) {
int rel = relationcircle(v);
if (rel >= 4) return 0.0;
if (rel <= 2) return min(area(), v.area());
double d = p.distance(v.p);
double hf = (r + v.r + d) / 2.0;
double ss = 2 * sqrt(hf * (hf - r) * (hf - v.r) * (hf - d));
double a1 = acos((r * r + d * d - v.r * v.r) / (2.0 * r * d));
a1 = a1 * r * r;
double a2 = acos((v.r * v.r + d * d - r * r) / (2.0 * v.r * d));
a2 = a2 * v.r * v.r;
return a1 + a2 - ss;
}
// 求圆和三角形pab的相交面积：p为圆心
double areatriangle(Point a, Point b) {
if (sgn((p - a) ^ (p - b)) == 0) return 0.0;
Point q[5];
int len = 0;
q[len++] = a;
Line l(a, b);
Point p1, p2;
if (pointcrossline(l, q[1], q[2]) == 2) {
if (sgn((a - q[1]) * (b - q[1])) < 0) q[len++] = q[1];
if (sgn((a - q[2]) * (b - q[2])) < 0) q[len++] = q[2];
}
q[len++] = b;
if (len == 4 && sgn((q[0] - q[1]) * (q[2] - q[1])) > 0)
swap(q[1], q[2]);
double res = 0;
for (int i = 0; i < len - 1; i++) {
if (relation(q[i]) == 0 || relation(q[i + 1]) == 0) {
double arg = p.rad(q[i], q[i + 1]);
res += r * r * arg / 2.0;
} else {
res += fabs((q[i] - p) ^ (q[i + 1] - p)) / 2.0;
}
}
return res;
}
};
/*
* n,p Line l : n为凸包上的点数，p为凸包上的所有点，l为凸包上的所有边
* input(int _n) - 输入一个大小为n的多边
* add(Point q) - 加入一个点
* getline() - 得到多边形所有边的线段
* cmp - 定义极角排序的顺序
* norm() - 进行极角排序
* getconvex(polygon &convex) - 得到凸包
* Graham(polygon &convex) - 得到凸包的另外一种方法
* isconvex() - 判断是不是凸的
* relationpoint(Point q) - 判断点和任意多边形的关系: 3 点上；2
* 边上；1 内部；0 外部 convexcut(Line u,polygon &po) -
* 直线u切割凸多边形左侧, 注意直线方向 gercircumference() -
* 得到多边形的周长 getarea() - 得到多边形的面积
* getdir() - 得到多边形的顶点的方向：1
* 表示逆时针，0表示顺时针 getbarycentre() -
* 得到多边形的重心 areacircle(circle c) - 多边形和圆交的面积
* relationcircle(circle c) - 多边形和圆关系：2 圆完全在多边形内；1
* 圆在多边形里面，碰到了多边形边界；0 其它
*/
struct polygon {
int n;
Point p[maxp];
Line l[maxp];
// 输入多边形
void input(int _n) {
n = _n;
for (int i = 0; i < n; i++) p[i].input();
}
// 加入一个点
void add(Point q) { p[n++] = q; }
// 得到多边形所有边的线段
void getline() {
for (int i = 0; i < n; i++) {
l[i] = Line(p[i], p[(i + 1) % n]);
}
}
// 极角排序
struct cmp {
Point p;
cmp(const Point &p0) { p = p0; }
bool operator()(const Point &aa, const Point &bb) {
Point a = aa, b = bb;
int d = sgn((a - p) ^ (b - p));
if (d == 0) {
return sgn(a.distance(p) - b.distance(p)) < 0;
}
return d > 0;
}
};
// 进行极角排序
// 首先需要找到最左下角的点
// 需要重载号好Point的 < 操作符(min函数要用)
void norm() {
Point mi = p[0];
for (int i = 1; i < n; i++)
mi = min(mi, p[i]); // min函数用Point的<重载
sort(p, p + n, cmp(mi));
}
// 得到凸包
// 得到的凸包里面的点编号是0 ~ n-1的
// 两种凸包的方法
// 注意如果有影响，要特判下所有点共点，或者共线的特殊情况
void getconvex(polygon &convex) {
sort(p, p + n);
convex.n = n;
for (int i = 0; i < min(n, 2); i++) {
convex.p[i] = p[i];
}
if (convex.n == 2 && (convex.p[0] == convex.p[1])) convex.n--; //特判
if (n <= 2) return;
int &top = convex.n;
top = 1;
for (int i = 2; i < n; i++) {
while (top && sgn((convex.p[top] - p[i]) ^
(convex.p[top - 1] - p[i])) <= 0)
top--;
convex.p[++top] = p[i];
}
int temp = top;
convex.p[++top] = p[n - 2];
for (int i = n - 3; i >= 0; i--) {
while (top != temp && sgn((convex.p[top] - p[i]) ^
(convex.p[top - 1] - p[i])) <= 0)
top--;
convex.p[++top] = p[i];
}
if (convex.n == 2 && (convex.p[0] == convex.p[1])) convex.n--; //特判
convex.norm(); //`原来得到的是顺时针的点，排序后逆时针`
}
// 得到凸包的另外一种方法
void Graham(polygon &convex) {
norm();
int &top = convex.n;
top = 0;
if (n == 1) {
top = 1;
convex.p[0] = p[0];
return;
}
if (n == 2) {
top = 2;
convex.p[0] = p[0];
convex.p[1] = p[1];
if (convex.p[0] == convex.p[1]) top--;
return;
}
convex.p[0] = p[0];
convex.p[1] = p[1];
top = 2;
for (int i = 2; i < n; i++) {
while (top > 1 && sgn((convex.p[top - 1] - convex.p[top - 2]) ^
(p[i] - convex.p[top - 2])) <= 0)
top--;
convex.p[top++] = p[i];
}
if (convex.n == 2 && (convex.p[0] == convex.p[1])) convex.n--; //特判
}
// 判断是不是凸的
bool isconvex() {
bool s[3];
memset(s, false, sizeof(s));
for (int i = 0; i < n; i++) {
int j = (i + 1) % n;
int k = (j + 1) % n;
s[sgn((p[j] - p[i]) ^ (p[k] - p[i])) + 1] = true;
if (s[0] && s[2]) return false;
}
return true;
}
// 判断点和任意多边形的关系
// 3 点上；2 边上；1 内部；0 外部`
int relationpoint(Point q) {
for (int i = 0; i < n; i++) {
if (p[i] == q) return 3;
}
getline();
for (int i = 0; i < n; i++) {
if (l[i].pointonseg(q)) return 2;
}
int cnt = 0;
for (int i = 0; i < n; i++) {
int j = (i + 1) % n;
int k = sgn((q - p[j]) ^ (p[i] - p[j]));
int u = sgn(p[i].y - q.y);
int v = sgn(p[j].y - q.y);
if (k > 0 && u < 0 && v >= 0) cnt++;
if (k < 0 && v < 0 && u >= 0) cnt--;
}
return cnt != 0;
}
// 直线u切割凸多边形左侧, 注意直线方向
void convexcut(Line u, polygon &po) {
int &top = po.n; //注意引用
top = 0;
for (int i = 0; i < n; i++) {
int d1 = sgn((u.e - u.s) ^ (p[i] - u.s));
int d2 = sgn((u.e - u.s) ^ (p[(i + 1) % n] - u.s));
if (d1 >= 0) po.p[top++] = p[i];
if (d1 * d2 < 0)
po.p[top++] = u.crosspoint(Line(p[i], p[(i + 1) % n]));
}
}
// 得到多边形的周长
double getcircumference() {
double sum = 0;
for (int i = 0; i < n; i++) {
sum += p[i].distance(p[(i + 1) % n]);
}
return sum;
}
// 得到多边形的面积
double getarea() {
double sum = 0;
for (int i = 0; i < n; i++) {
sum += (p[i] ^ p[(i + 1) % n]);
}
return fabs(sum) / 2;
}
// 得到多边形的顶点的方向：1 表示逆时针，0表示顺时针
bool getdir() {
double sum = 0;
for (int i = 0; i < n; i++) sum += (p[i] ^ p[(i + 1) % n]);
if (sgn(sum) > 0) return 1;
return 0;
}
// 得到多边形的重心
Point getbarycentre() {
Point ret(0, 0);
double area = 0;
for (int i = 1; i < n - 1; i++) {
double tmp = (p[i] - p[0]) ^ (p[i + 1] - p[0]);
if (sgn(tmp) == 0) continue;
area += tmp;
ret.x += (p[0].x + p[i].x + p[i + 1].x) / 3 * tmp;
ret.y += (p[0].y + p[i].y + p[i + 1].y) / 3 * tmp;
}
if (sgn(area)) ret = ret / area;
return ret;
}
// 多边形和圆交的面积
double areacircle(circle c) {
double ans = 0;
for (int i = 0; i < n; i++) {
int j = (i + 1) % n;
if (sgn((p[j] - c.p) ^ (p[i] - c.p)) >= 0)
ans += c.areatriangle(p[i], p[j]);
else
ans -= c.areatriangle(p[i], p[j]);
}
return fabs(ans);
}
// 多边形和圆关系：2 圆完全在多边形内；1 圆在多边形里面，碰到了多边形边界；0
// 其它
int relationcircle(circle c) {
getline();
int x = 2;
if (relationpoint(c.p) != 1) return 0; //圆心不在内部
for (int i = 0; i < n; i++) {
if (c.relationseg(l[i]) == 2) return 0;
if (c.relationseg(l[i]) == 1) x = 1;
}
return x;
}
};
// AB X AC
double cross(Point A, Point B, Point C) { return (B - A) ^ (C - A); }
// AB*AC
double dot(Point A, Point B, Point C) { return (B - A) * (C - A); }
// 最小矩形面积覆盖：用最小的矩形去覆盖一个凸包，返回矩形的面积
// A 必须是凸包(而且是逆时针顺序)
double minRectangleCover(polygon A) {
//要特判A.n < 3的情况
if (A.n < 3) return 0.0;
A.p[A.n] = A.p[0];
double ans = -1;
int r = 1, p = 1, q;
for (int i = 0; i < A.n; i++) {
//卡出离边A.p[i] - A.p[i+1]最远的点
while (sgn(cross(A.p[i], A.p[i + 1], A.p[r + 1]) -
cross(A.p[i], A.p[i + 1], A.p[r])) >= 0)
r = (r + 1) % A.n;
//卡出A.p[i] - A.p[i+1]方向上正向n最远的点
while (sgn(dot(A.p[i], A.p[i + 1], A.p[p + 1]) -
dot(A.p[i], A.p[i + 1], A.p[p])) >= 0)
p = (p + 1) % A.n;
if (i == 0) q = p;
//卡出A.p[i] - A.p[i+1]方向上负向最远的点
while (sgn(dot(A.p[i], A.p[i + 1], A.p[q + 1]) -
dot(A.p[i], A.p[i + 1], A.p[q])) <= 0)
q = (q + 1) % A.n;
double d = (A.p[i] - A.p[i + 1]).len2();
double tmp = cross(A.p[i], A.p[i + 1], A.p[r]) *
(dot(A.p[i], A.p[i + 1], A.p[p]) -
dot(A.p[i], A.p[i + 1], A.p[q])) /
d;
if (ans < 0 || ans > tmp) ans = tmp;
}
return ans;
}
// 直线切凸多边形，返回被切割的左侧
// 多边形是逆时针的，在q1q2的左侧
vector<Point> convexCut(const vector<Point> &ps, Point q1, Point q2) {
vector<Point> qs;
int n = ps.size();
for (int i = 0; i < n; i++) {
Point p1 = ps[i], p2 = ps[(i + 1) % n];
int d1 = sgn((q2 - q1) ^ (p1 - q1)), d2 = sgn((q2 - q1) ^ (p2 - q1));
if (d1 >= 0) qs.push_back(p1);
if (d1 * d2 < 0) qs.push_back(Line(p1, p2).crosspoint(Line(q1, q2)));
}
return qs;
}
// 半平面交
//***************************
// 定义一个半平面
struct halfplane : public Line {
double angle;
halfplane() {}
//`表示向量s->e逆时针(左侧)的半平面`
halfplane(Point _s, Point _e) {
s = _s;
e = _e;
}
halfplane(Line v) {
s = v.s;
e = v.e;
}
void calcangle() { angle = atan2(e.y - s.y, e.x - s.x); }
bool operator<(const halfplane &b) const { return angle < b.angle; }
};
// 定义一组半平面
struct halfplanes {
int n;
halfplane hp[maxp];
Point p[maxp];
int que[maxp];
int st, ed;
void push(halfplane tmp) { hp[n++] = tmp; }
//去重
void unique() {
int m = 1;
for (int i = 1; i < n; i++) {
if (sgn(hp[i].angle - hp[i - 1].angle) != 0)
hp[m++] = hp[i];
else if (sgn((hp[m - 1].e - hp[m - 1].s) ^
(hp[i].s - hp[m - 1].s)) > 0)
hp[m - 1] = hp[i];
}
n = m;
}
bool halfplaneinsert() {
for (int i = 0; i < n; i++) hp[i].calcangle();
sort(hp, hp + n);
unique();
que[st = 0] = 0;
que[ed = 1] = 1;
p[1] = hp[0].crosspoint(hp[1]);
for (int i = 2; i < n; i++) {
while (st < ed && sgn((hp[i].e - hp[i].s) ^ (p[ed] - hp[i].s)) < 0)
ed--;
while (st < ed &&
sgn((hp[i].e - hp[i].s) ^ (p[st + 1] - hp[i].s)) < 0)
st++;
que[++ed] = i;
if (hp[i].parallel(hp[que[ed - 1]])) return false;
p[ed] = hp[i].crosspoint(hp[que[ed - 1]]);
}
while (st < ed && sgn((hp[que[st]].e - hp[que[st]].s) ^
(p[ed] - hp[que[st]].s)) < 0)
ed--;
while (st < ed && sgn((hp[que[ed]].e - hp[que[ed]].s) ^
(p[st + 1] - hp[que[ed]].s)) < 0)
st++;
if (st + 1 >= ed) return false;
return true;
}
// 得到最后半平面交得到的凸多边形
// 需要先调用halfplaneinsert() 且返回true
void getconvex(polygon &con) {
p[st] = hp[que[st]].crosspoint(hp[que[ed]]);
con.n = ed - st + 1;
for (int j = st, i = 0; j <= ed; i++, j++) con.p[i] = p[j];
}
};
//***************************
// 求多个圆面积的交和并
const int maxn = 1010;
struct circles {
circle c[maxn]; // 存储所有的圆
double ans[maxn]; // ans[i]表示被覆盖了i次的面积
double pre[maxn];
int n;
circles() {}
void add(circle cc) { c[n++] = cc; }
// 判断圆x是否包含在圆y中
bool inner(circle x, circle y) {
if (x.relationcircle(y) != 1) return 0;
return sgn(x.r - y.r) <= 0 ? 1 : 0;
}
// 求圆的面积并去掉内含的圆
void init_or() {
bool mark[maxn] = {0};
int i, j, k = 0;
for (i = 0; i < n; i++) {
for (j = 0; j < n; j++)
if (i != j && !mark[j]) {
if ((c[i] == c[j]) || inner(c[i], c[j])) break;
}
if (j < n) mark[i] = 1;
}
for (i = 0; i < n; i++)
if (!mark[i]) c[k++] = c[i];
n = k;
}
// 求圆的面积交去掉内含的圆
void init_add() {
int i, j, k;
bool mark[maxn] = {0};
for (i = 0; i < n; i++) {
for (j = 0; j < n; j++)
if (i != j && !mark[j]) {
if ((c[i] == c[j]) || inner(c[j], c[i])) break;
}
if (j < n) mark[i] = 1;
}
for (i = 0; i < n; i++)
if (!mark[i]) c[k++] = c[i];
n = k;
}
// 半径为r的圆，弧度为th对应的弓形的面积
double areaarc(double th, double r) { return 0.5 * r * r * (th - sin(th)); }
// 测试SPOJVCIRCLES SPOJCIRUT
// SPOJVCIRCLES求n个圆并的面积，需要加上init_or()去掉重复圆（否则WA）
// SPOJCIRUT 是求被覆盖k次的面积，不能加init_or()
// 对于求覆盖多少次面积的问题，不能解决相同圆，而且不能init_or()
// 求多圆面积并，需要init_or,其中一个目的就是去掉相同圆
// 最后求多个圆的面积并，就是把ans[1] + ... +
// ans[n]的和，求被覆盖k次的面积和就是ans[k]
void getarea() {
memset(ans, 0, sizeof(ans));
vector<pair<double, int>> v;
for (int i = 0; i < n; i++) {
v.clear();
v.push_back(make_pair(-pi, 1));
v.push_back(make_pair(pi, -1));
for (int j = 0; j < n; j++)
if (i != j) {
Point q = (c[j].p - c[i].p);
double ab = q.len(), ac = c[i].r, bc = c[j].r;
if (sgn(ab + ac - bc) <= 0) {
v.push_back(make_pair(-pi, 1));
v.push_back(make_pair(pi, -1));
continue;
}
if (sgn(ab + bc - ac) <= 0) continue;
if (sgn(ab - ac - bc) > 0) continue;
double th = atan2(q.y, q.x),
fai = acos((ac * ac + ab * ab - bc * bc) /
(2.0 * ac * ab));
double a0 = th - fai;
if (sgn(a0 + pi) < 0) a0 += 2 * pi;
double a1 = th + fai;
if (sgn(a1 - pi) > 0) a1 -= 2 * pi;
if (sgn(a0 - a1) > 0) {
v.push_back(make_pair(a0, 1));
v.push_back(make_pair(pi, -1));
v.push_back(make_pair(-pi, 1));
v.push_back(make_pair(a1, -1));
} else {
v.push_back(make_pair(a0, 1));
v.push_back(make_pair(a1, -1));
}
}
sort(v.begin(), v.end());
int cur = 0;
for (int j = 0; j < v.size(); j++) {
if (cur && sgn(v[j].first - pre[cur])) {
ans[cur] += areaarc(v[j].first - pre[cur], c[i].r);
ans[cur] +=
0.5 * (Point(c[i].p.x + c[i].r * cos(pre[cur]),
c[i].p.y + c[i].r * sin(pre[cur])) ^
Point(c[i].p.x + c[i].r * cos(v[j].first),
c[i].p.y + c[i].r * sin(v[j].first)));
}
cur += v[j].second;
pre[cur] = v[j].first;
}
}
for (int i = 1; i < n; i++) ans[i] -= ans[i + 1];
}
};
2.2 三维计算几何
#include <bits/stdc++.h>
using namespace std;
const double eps = 1e-8;
int sgn(double x) {
if (fabs(x) < eps) return 0;
if (x < 0)
return -1;
else
return 1;
}
// 三维的向量
/*
* Point3
* Point3() - 空构造
* Point3(double _x = 0, double _y = 0, double _z = 0) - 普通构造
* input() - 浮点输入
* output() - %.2f 输出
* operator == - 判断是否相等
* operator < - 判断是否更小，先比较x，然后比较y
* operator - - 返回点p减去一个点的
* operator ^ - 叉乘
* operator * - 点乘
* len() - 得到向量p的长度
* len2() - 得到向量p长度的平方
* distance(Point3 p) - 得到点p和其他点之间的距离
* operator + Point3 b - 返回向量加的结果
* operator * double k - 把点p的x和y都放大k倍
* operator / double k - 把点p的x和y都缩小k倍
* rad(Point3 a,Point3 b)- 得到点p去看点a和点b的夹角的弧度
* trunc(double r) - 把向量p进行放缩，使得其长度为r
*/
struct Point3 {
double x, y, z;
Point3(double _x = 0, double _y = 0, double _z = 0) {
x = _x;
y = _y;
z = _z;
}
// 输入
void input() { scanf("%lf%lf%lf", &x, &y, &z); }
// 输出
void output() { printf("%.2lf %.2lf %.2lf\n", x, y, z); }
// 判断向量是否相等
bool operator==(const Point3 &b) const {
return sgn(x - b.x) == 0 && sgn(y - b.y) == 0 && sgn(z - b.z) == 0;
}
// 判断向量大小关系
bool operator<(const Point3 &b) const {
return sgn(x - b.x) == 0
? (sgn(y - b.y) == 0 ? sgn(z - b.z) < 0 : y < b.y)
: x < b.x;
}
// 计算向量的长度
double len() { return sqrt(x * x + y * y + z * z); }
// 计算向量长度的平方
double len2() { return x * x + y * y + z * z; }
// 计算向量p到向量b的距离
double distance(const Point3 &b) const {
return sqrt((x - b.x) * (x - b.x) + (y - b.y) * (y - b.y) +
(z - b.z) * (z - b.z));
}
// 计算向量p-向量b
Point3 operator-(const Point3 &b) const {
return Point3(x - b.x, y - b.y, z - b.z);
}
// 计算向量+
Point3 operator+(const Point3 &b) const {
return Point3(x + b.x, y + b.y, z + b.z);
}
// 把向量p放大k倍
Point3 operator*(const double &k) const {
return Point3(x * k, y * k, z * k);
}
// 把向量p缩小k倍
Point3 operator/(const double &k) const {
return Point3(x / k, y / k, z / k);
}
//点乘
double operator*(const Point3 &b) const {
return x * b.x + y * b.y + z * b.z;
}
//叉乘
Point3 operator^(const Point3 &b) const {
return Point3(y * b.z - z * b.y, z * b.x - x * b.z, x * b.y - y * b.x);
}
// 计算向量pa和pb的夹角
double rad(Point3 a, Point3 b) {
Point3 p = (*this);
return acos(((a - p) * (b - p)) / (a.distance(p) * b.distance(p)));
}
//变换长度
Point3 trunc(double r) {
double l = len();
if (!sgn(l)) return *this;
r /= l;
return Point3(x * r, y * r, z * r);
}
};
// 直线
/*
* Line3
* Line3() - 空构造
* Line3(Point3 _s, Point3 _e) - 普通构造
* operator== - 判断直线是否相等
* input() - 输入
* length() - 计算线段两个端点的距离
* dispointtoline(Point3 p) - 点到直线距离
* dispointtoseg(Point3 p) - 点到线段距离
* lineprog(Point3 p) - 返回点p在直线上的投影
* rotate(Point3 p, double ang) - p绕此向量逆时针arg角度
* pointonseg(Point3 p) - 点在线段上
*/
struct Line3 {
Point3 s, e;
Line3() {}
Line3(Point3 _s, Point3 _e) {
s = _s;
e = _e;
}
// 判断直线是否相等
bool operator==(const Line3 v) { return (s == v.s) && (e == v.e); }
void input() {
s.input();
e.input();
}
// 计算线段两个端点的距离
double length() { return s.distance(e); }
//点到直线距离
double dispointtoline(Point3 p) {
return ((e - s) ^ (p - s)).len() / s.distance(e);
}
//点到线段距离
double dispointtoseg(Point3 p) {
if (sgn((p - s) * (e - s)) < 0 || sgn((p - e) * (s - e)) < 0)
return min(p.distance(s), e.distance(p));
return dispointtoline(p);
}
// 返回点p在直线上的投影
Point3 lineprog(Point3 p) {
return s + (((e - s) * ((e - s) * (p - s))) / ((e - s).len2()));
}
// p绕此向量逆时针arg角度
Point3 rotate(Point3 p, double ang) {
if (sgn(((s - p) ^ (e - p)).len()) == 0) return p;
Point3 f1 = (e - s) ^ (p - s);
Point3 f2 = (e - s) ^ (f1);
double len = ((s - p) ^ (e - p)).len() / s.distance(e);
f1 = f1.trunc(len);
f2 = f2.trunc(len);
Point3 h = p + f2;
Point3 pp = h + f1;
return h + ((p - h) * cos(ang)) + ((pp - h) * sin(ang));
}
// 点在线段上
bool pointonseg(Point3 p) {
return sgn(((s - p) ^ (e - p)).len()) == 0 &&
sgn((s - p) * (e - p)) == 0;
}
};
// 平面
/*
* Plane() - 空构造
* Plane(Point3 _a, Point3 _b, Point3 _c) - 三个点构造一个平面
* pvec() - 计算法向量
* Plane(double _a, double _b, double _c, double _d) - 平面方程ax+by+cz+d = 0构造平面
* pointonplane(Point3 p) - 点在平面上的判断
* angleplane(Plane f) - 两平面夹角
* crossline(Line3 u, Point3 &p) - 平面和直线的交点，返回值是交点个数，交点存储在p内
* pointtoplane(Point3 p) - 点到平面最近点(也就是投影)
* crossplane(Plane f, Line3 &u) - 平面和平面的交线
*/
struct Plane {
Point3 a, b, c, o; // 平面上的三个点，以及法向量
Plane() {}
Plane(Point3 _a, Point3 _b, Point3 _c) {
a = _a;
b = _b;
c = _c;
o = pvec();
}
// 计算法向量
Point3 pvec() { return (b - a) ^ (c - a); }
// ax+by+cz+d = 0
Plane(double _a, double _b, double _c, double _d) {
o = Point3(_a, _b, _c);
if (sgn(_a) != 0)
a = Point3((-_d - _c - _b) / _a, 1, 1);
else if (sgn(_b) != 0)
a = Point3(1, (-_d - _c - _a) / _b, 1);
else if (sgn(_c) != 0)
a = Point3(1, 1, (-_d - _a - _b) / _c);
}
// 点在平面上的判断
bool pointonplane(Point3 p) { return sgn((p - a) * o) == 0; }
// 两平面夹角
double angleplane(Plane f) { return acos(o * f.o) / (o.len() * f.o.len()); }
// 平面和直线的交点，返回值是交点个数，交点存储在p内
int crossline(Line3 u, Point3 &p) {
double x = o * (u.e - a);
double y = o * (u.s - a);
double d = x - y;
if (sgn(d) == 0) return 0;
p = ((u.s * x) - (u.e * y)) / d;
return 1;
}
// 点到平面最近点(也就是投影)
Point3 pointtoplane(Point3 p) {
Line3 u = Line3(p, p + o);
crossline(u, p);
return p;
}
// 平面和平面的交线
int crossplane(Plane f, Line3 &u) {
Point3 oo = o ^ f.o;
Point3 v = o ^ oo;
double d = fabs(f.o * v);
if (sgn(d) == 0) return 0;
Point3 q = a + (v * (f.o * (f.a - a)) / d);
u = Line3(q, q + oo);
return 1;
}
};
2.3 平面最近点对
#include <bits/stdc++.h>
using namespace std;
const int MAXN = 100010;
const double eps = 1e-8;
const double INF = 1e20;
struct Point {
double x, y;
void input() { scanf("%lf%lf", &x, &y); }
};
double dist(Point a, Point b) {
return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
}
Point p[MAXN];
Point tmpt[MAXN];
bool cmpx(Point a, Point b) { return a.x < b.x || (a.x == b.x && a.y < b.y); }
bool cmpy(Point a, Point b) { return a.y < b.y || (a.y == b.y && a.x < b.x); }
double Closest_Pair(int left, int right) {
double d = INF;
if (left == right) return d;
if (left + 1 == right) return dist(p[left], p[right]);
int mid = (left + right) / 2;
double d1 = Closest_Pair(left, mid);
double d2 = Closest_Pair(mid + 1, right);
d = min(d1, d2);
int cnt = 0;
for (int i = left; i <= right; i++) {
if (fabs(p[mid].x - p[i].x) <= d) tmpt[cnt++] = p[i];
}
sort(tmpt, tmpt + cnt, cmpy);
for (int i = 0; i < cnt; i++) {
for (int j = i + 1; j < cnt && tmpt[j].y - tmpt[i].y < d; j++)
d = min(d, dist(tmpt[i], tmpt[j]));
}
return d;
}
int main() {
int n;
while (scanf("%d", &n) == 1 && n) {
for (int i = 0; i < n; i++) p[i].input();
sort(p, p + n, cmpx);
printf("%.2lf\n", Closest_Pair(0, n - 1));
}
return 0;
}
2.4 三维凸包
#include <bits/stdc++.h>
using namespace std;
const double eps = 1e-8;
const int MAXN = 550;
int sgn(double x){
if(fabs(x) < eps)return 0;
if(x < 0)return -1;
else return 1;
}
struct Point3{
double x,y,z;
Point3(double _x = 0, double _y = 0, double _z = 0){
x = _x;
y = _y;
z = _z;
}
void input(){
scanf("%lf%lf%lf",&x,&y,&z);
}
bool operator ==(const Point3 &b)const{
return sgn(x-b.x) == 0 && sgn(y-b.y) == 0 && sgn(z-b.z) == 0;
}
double len(){
return sqrt(x*x+y*y+z*z);
}
double len2(){
return x*x+y*y+z*z;
}
double distance(const Point3 &b)const{
return sqrt((x-b.x)*(x-b.x)+(y-b.y)*(y-b.y)+(z-b.z)*(z-b.z));
}
Point3 operator -(const Point3 &b)const{
return Point3(x-b.x,y-b.y,z-b.z);
}
Point3 operator +(const Point3 &b)const{
return Point3(x+b.x,y+b.y,z+b.z);
}
Point3 operator *(const double &k)const{
return Point3(x*k,y*k,z*k);
}
Point3 operator /(const double &k)const{
return Point3(x/k,y/k,z/k);
}
//点乘
double operator *(const Point3 &b)const{
return x*b.x + y*b.y + z*b.z;
}
//叉乘
Point3 operator ^(const Point3 &b)const{
return Point3(y*b.z-z*b.y,z*b.x-x*b.z,x*b.y-y*b.x);
}
};
struct CH3D{
struct face{
//表示凸包一个面上的三个点的编号
int a,b,c;
//表示该面是否属于最终的凸包上的面
bool ok;
};
//初始顶点数
int n;
Point3 P[MAXN];
//凸包表面的三角形数
int num;
//凸包表面的三角形
face F[8*MAXN];
int g[MAXN][MAXN];
//叉乘
Point3 cross(const Point3 &a,const Point3 &b,const Point3 &c){
return (b-a)^(c-a);
}
// 三角形面积*2
double area(Point3 a,Point3 b,Point3 c){
return ((b-a)^(c-a)).len();
}
// 四面体有向面积*6
double volume(Point3 a,Point3 b,Point3 c,Point3 d){
return ((b-a)^(c-a))*(d-a);
}
// 正：点在面同向
double dblcmp(Point3 &p,face &f){
Point3 p1 = P[f.b] - P[f.a];
Point3 p2 = P[f.c] - P[f.a];
Point3 p3 = p - P[f.a];
return (p1^p2)*p3;
}
void deal(int p,int a,int b){
int f = g[a][b];
face add;
if(F[f].ok){
if(dblcmp(P[p],F[f]) > eps)
dfs(p,f);
else {
add.a = b;
add.b = a;
add.c = p;
add.ok = true;
g[p][b] = g[a][p] = g[b][a] = num;
F[num++] = add;
}
}
}
//递归搜索所有应该从凸包内删除的面
void dfs(int p,int now){
F[now].ok = false;
deal(p,F[now].b,F[now].a);
deal(p,F[now].c,F[now].b);
deal(p,F[now].a,F[now].c);
}
bool same(int s,int t){
Point3 &a = P[F[s].a];
Point3 &b = P[F[s].b];
Point3 &c = P[F[s].c];
return fabs(volume(a,b,c,P[F[t].a])) < eps &&
fabs(volume(a,b,c,P[F[t].b])) < eps &&
fabs(volume(a,b,c,P[F[t].c])) < eps;
}
//构建三维凸包
void create(){
num = 0;
face add;
//***********************************
//此段是为了保证前四个点不共面
bool flag = true;
for(int i = 1;i < n;i++){
if(!(P[0] == P[i])){
swap(P[1],P[i]);
flag = false;
break;
}
}
if(flag)return;
flag = true;
for(int i = 2;i < n;i++){
if( ((P[1]-P[0])^(P[i]-P[0])).len() > eps ){
swap(P[2],P[i]);
flag = false;
break;
}
}
if(flag)return;
flag = true;
for(int i = 3;i < n;i++){
if(fabs( ((P[1]-P[0])^(P[2]-P[0]))*(P[i]-P[0]) ) > eps){
swap(P[3],P[i]);
flag = false;
break;
}
}
if(flag)return;
//**********************************
for(int i = 0;i < 4;i++){
add.a = (i+1)%4;
add.b = (i+2)%4;
add.c = (i+3)%4;
add.ok = true;
if(dblcmp(P[i],add) > 0)swap(add.b,add.c);
g[add.a][add.b] = g[add.b][add.c] = g[add.c][add.a] = num;
F[num++] = add;
}
for(int i = 4;i < n;i++)
for(int j = 0;j < num;j++)
if(F[j].ok && dblcmp(P[i],F[j]) > eps){
dfs(i,j);
break;
}
int tmp = num;
num = 0;
for(int i = 0;i < tmp;i++)
if(F[i].ok)
F[num++] = F[i];
}
//表面积
double area(){
double res = 0;
if(n == 3){
Point3 p = cross(P[0],P[1],P[2]);
return p.len()/2;
}
for(int i = 0;i < num;i++)
res += area(P[F[i].a],P[F[i].b],P[F[i].c]);
return res/2.0;
}
double volume(){
double res = 0;
Point3 tmp = Point3(0,0,0);
for(int i = 0;i < num;i++)
res += volume(tmp,P[F[i].a],P[F[i].b],P[F[i].c]);
return fabs(res/6);
}
//表面三角形个数
int triangle(){
return num;
}
//表面多边形个数
int polygon(){
int res = 0;
for(int i = 0;i < num;i++){
bool flag = true;
for(int j = 0;j < i;j++)
if(same(i,j)){
flag = 0;
break;
}
res += flag;
}
return res;
}
//重心
Point3 barycenter(){
Point3 ans = Point3(0,0,0);
Point3 o = Point3(0,0,0);
double all = 0;
for(int i = 0;i < num;i++){
double vol = volume(o,P[F[i].a],P[F[i].b],P[F[i].c]);
ans = ans + (((o+P[F[i].a]+P[F[i].b]+P[F[i].c])/4.0)*vol);
all += vol;
}
ans = ans/all;
return ans;
}
//点到面的距离
double ptoface(Point3 p,int i){
double tmp1 = fabs(volume(P[F[i].a],P[F[i].b],P[F[i].c],p));
double tmp2 = ((P[F[i].b]-P[F[i].a])^(P[F[i].c]-P[F[i].a])).len();
return tmp1/tmp2;
}
};
CH3D hull;
int main()
{
while(scanf("%d",&hull.n) == 1){
for(int i = 0;i < hull.n;i++)hull.P[i].input();
hull.create();
Point3 p = hull.barycenter();
double ans = 1e20;
for(int i = 0;i < hull.num;i++)
ans = min(ans,hull.ptoface(p,i));
printf("%.3lf\n",ans);
}
return 0;
}
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
