计算几何基础

最新推荐文章于 2024-06-23 20:47:20 发布
最新推荐文章于 2024-06-23 20:47:20 发布
阅读量324 收藏 1
点赞数
分类专栏： ACM--计算几何
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/m0_49959202/article/details/110849110
版权

文章目录

计算几何基础

本文大部分内容来自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​>=x1​x2​+y1​y2​

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​=∣∣∣∣​x1​x2​​y1​y2​​∣∣∣∣​=x1​y2​−x2​y1​

∣ 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′=x1​cosθ−y1​sinθy′=x1​sinθ+y1​cosθ​

1.3 直线与线段

1.3.1 线段相交问题
1.png

线段 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 点到直线的距离
2.png

要计算点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​:a1​x−b1​y=a1​x1​−b1​y1​

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​:a2​x−b2​y=a2​x2​−b2​y2​

联立两直线的方程，由克拉默法则得，方程组的解为
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} P0​P2​…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} P0​P2​…Pn−1​ ，点的顺序为顺时针或逆时针，以及多边形内一点A，把多边形切割成如下所示n个三角形

3.png

多边形的面积等于所有三角形（有向）面积之和，代入坐标 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=∣∣∣∣∣​21​i=0∑n−2​(xi​yi+1​−xi+1​yi​)+21​(xn−1​y0​−x0​yn−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 convex​K

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 圆的参数方程
4.png

以 ( 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 两圆交点
5.png

设两圆 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 ∠AC1​D与 ∠ B C 1 D \angle BC_{1}D ∠BC1​D，进而通过圆的参数方程即可求得
∠ 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} ∠AC1​D=∠C2​C1​D+∠AC1​C2​

∠ 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} ∠BC1​D=∠C2​C1​D−∠AC1​C2​

∠ C 2 C 1 D \angle C_{2}C_{1}D ∠C2​C1​D可以通过 C 1 , C 2 C_{1},C_{2} C1​,C2​的坐标求得，而 ∠ A C 1 C 2 \angle AC_{1}C_{2} ∠AC1​C2​可以通过 Δ A C 1 C 2 \Delta AC_{1}C_{2} ΔAC1​C2​上的余弦定理求得

对于情形③和情形⑤，上述方法求得的两点坐标是相同的，即为切点的坐标

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;
}
「已注销」
关注
  • 0
    点赞
  • 踩
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
专栏目录
arcgis 批量计算几何_GIS算法基础（一） 计算几何基础（上）
weixin_39997310的博客
01-12 1028
最近在学习GIS算法，在学习过程中，想把一些经典的算法或者思想记录下来，分享给大家计算几何基础本来是计算机图形学的内容，但是GIS在图像处理中是离不开计算机处理的，所以GIS算法基础第一个应该是计算几何基础。如何把空间实体的点线面以及他们之间的关系(例如，相交，包含，邻接，叠置等等)用计算机语言表达出来，应该是这里解决的第一个问题。前人已经有过这样这样的解决模型---DE-9IM模型(Dimens...
ACM 计算几何基础模板
qq_39561693的博客
03-25 1136
在ACM计算几何中，我们常常重复用到许多方法，例如求点距，求直线方程，等等。不妨利用模板总结之。这些模板为我从网上学习后加入自己的理解糅合得来qwq，看了很多的博客，查了很多书籍，其中印象最深的是林夕林夕关于计算几何模板的总结，大家也可以去参考他的应该比我的更详细qwq（不过大佬最近好像出没比较少了） 所以我这个主要还是留着自己参考吧，能帮到他人再好不过了 为了不产生歧义，方便理解，本模板使用笛卡尔坐标系，封装的多是数学中常见的函数。 前置知识：高中数学 文章目录基础声明点向量直线 基础声明 首先，打acm
建图学习一
weixin_37781153的博客
12-26 111
使用目的：定位、导航、避障、重建、交互 类型：稠密、稀疏 单目稠密建图
计算几何基础入门详解
slass-mot-odet
08-16 2172
我把计算几何分为了六大块 一、计算几何基础(点，线，面) 二、计算几何凸包 三、计算几何圆形 四、计算几何旋转卡壳 五、计算几何半平面交 六、三维计算几何(本博客内容不予介绍) 一、计算几何基础 1、最开始要提一下的还是精度问题，这在计算几何中可以说非常关键！！！ const double eps=1e-8; //我一般都用eps表示，控制在1e-8左右，当然有些题会有它的要求 int sgn(d...
计算几何基础【用图来助你理解几何算法】
qq_45249273的博客
03-28 7411
计算几何基础 1.基本概念 1.1计算几何的引入 ​ 计算几何是几何学的一个重要分支，也是计算机科学的一个分支,研究解决几何问题的算法。在现代工程与数学、计算机图形学、机器人学、VLSI设计、计算机辅助设计等学科领域中,计算几何都有重要应用。 计算几何问题的输入一般是关于一组几何物体(如点、线)的描述;输出常常是有关这些物体相关问题的回答，如直线是否相交、点围成的面积等问题。 1.2 浮点数造成的误差 ​ 在学习计算几何之前，我们来回顾一下高中知识，我们求直线方程，或者是交点时，经常出现根号多少多少，以至我
计算几何基础汇总
WHYZ2017
05-20 530
01:    [POJ3348] Cows02:    [POJ1265] Area03:    [POJ3714] Raid
计算几何基础知识
千古的博客
05-20 1590
向量的叉乘 三维叉积：两个向量 a=(x1,y1,z1),b=(x2,y2,z2) 的叉积的结果是一个向量c。记作c=a×b。 其中i,j,k是三个轴上的单位向量。展开结果得到c=(y1-z2,x2z1-x1z2,x1y2-x2y1)。展开方式有三阶行列式的展开方式得到(见下)。其中的i,j,k作为单位向量合并，并写在对应坐标处。 根据叉积的计算式子，c 的
计算几何基础工具 叉积介绍
一七得七
08-09 1102
这篇文章将介绍计算几何中一个基础而重要的工具——叉积。   在这之前，我们先要解决一些基本问题。                            点、线段的代码表示   结合结构体或者类的知识，这里其实很好理解，但是在具体的编码过程中怎么写，这个根据个人有很大的灵活性。   比如对于点，就可以用如下的语句表达：   struct point {  double x , y;
计算几何基础———— 空间向量的基本运算重载，点乘和叉乘
Harukaze
11-30 1922
定义向量 以OSG （OpenSceneGraph）为例 定于向量的数据类型 typedef float value_type 向量的维度 enum{num_components = 3; //这里以三维为例子} 三维向量的各个分量可以用一维数组表示 ; value_type _v[3]; ##重载运算符 ： 重载等于 inline bool operator == (const Vec3f& v) const { return _v[0] == v._v[0] && _v[
算法课程笔记——计算几何基础
2302_79123319的博客
06-23 832
算法课程笔记——计算几何基础。
计算几何基础知识提纲
12-20
"计算几何基础知识提纲" 计算几何是计算机科学和数学的交叉领域，涉及到点、向量、线、面、体等几何对象的计算和处理。计算几何基础知识提纲是计算几何的基础知识，包括向量、坐标系、直线与多边形等内容。 一、...
计算几何基础凸包学习笔记
12-25
计算几何基础凸包学习笔记
计算机几何基础
12-15
学习计算机语言的基础文档，对于思维的培养有很大的好处。
ACM计算几何基础模板，杭电老学长精心整理.pdf
03-25
包含点、线、多边形、凸多边形、圆、半平面交等算法的板子
1111111111111111111111
08-13
1111111111111111
javascript笔试题参考整理(答案)(可编辑修改word版).docx
08-13
javascript笔试题参考整理(答案)(可编辑修改word版)
财务收支表-公式计算查询.xlsx
08-13
工资表，财务报表，对账表，付款申请，财务报告，费用支出表，财务收支 适用人群：学习不同技术领域的小白或进阶学习者；可作为毕设项目、课程设计、大作业、工程实训或初期项目立项。
集团信息化顶层规划 .pptx
最新发布
08-13
集团信息化顶层规划 .pptx
计算几何基础：三角形的面积计算
计算几何基础——三角形 计算几何是计算机科学中的一门重要学科，它研究如何使用计算机对几何对象进行处理和分析。在计算机科学中，计算几何有很多应用，如计算机视觉、机器人学、计算机辅助设计等。今天，我们将从...
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值