【ACM】 计算几何

前言

因为没有拿到牌子，心情低迷，随缘更新..

计算几何

目录

1、平面几何前置知识

///非科学计数法输出 
cout.setf(ios::fixed);
///设置输出小数点 
cout << setprecision(4) << x <<endl;

2、点

2.1、点的表示

在计算几何的题目中，基本上都是实数计算，所以我们一般都使用精度更高的 double 来进行计算。double 类型的数据在读入的时候要使用 %lf ，在输出的时候要使用 %f 。

struct Point
{
    double x,y;
    Point(){}							        /// 无参构造
    Point(double _x,double _y):x(-x),y(_y){}    /// 有参构造  
}P[N];

2.2、两点之间的距离

两点距离公式： \(dis = \sqrt{(x_1-x_2)^2 + (y_1-y_2)^2}\)

double distance(Point P){
        hypot(x-P.x,y-P.y); /// 库函数，求根号下，两数平方和
}

3、向量

3.1、计算机中的向量

我们一般表示一个向量也是直接用一个点来表示，方向是从原点指向（x,y）坐标

比如说 Point (4,2) 这个点就表示的是从原点指向 B 的一个向量，因为向量是可以平移的，所以为了方便计算，我们直接用一个点来表示一个向量。

I3Kkm4.png

在计算中的，我们处理区间的运算，多数是在坐标原点处理的，可能他真实的情况是 \((A_1,C_1)\) ，\((A_2,B_2)\)​​ 这样的情况，但是我们都将其平移到原点处理

I8ZrOH.png

3.1、点积（Dot Product）

向量的基本运算是 点积 和 叉积 ，计算几何的各种操作几乎都基于这两种运算。

已知两个 \(\vec{a},\vec{b}\)​ 向量 ，它们的夹角为 \(\theta\)​ ，那么：$ a \cdot b = |a||b|\cos \theta $​

就是这两个向量的 数量积，也叫 点积 或 内积 。点积的几何意义为： \(\vec{a}\) 在 \(\vec{b}\) 上的投影乘以 \(\vec{b}\)​ 的模长

IGpXn0.png

3.1.1、向量间的投影

由点积公式：$ a \cdot b = |a||b|\cos \theta $​

double operator * (Point b){
	return x*b.x + y*b.y;
}

我们可以发现 \(|a|\cos \theta\)​​ ，\(|b|\cos \theta\)​​​​​ 都为一个向量到另一个向量的投影。

我们就可以通过点积来计算相互之间的投影:

\(|b|\cos \theta = \frac{ \vec{a} \cdot \vec{b}}{|\vec{a}|}\)​​

\(|a|\cos \theta = \frac{ \vec{a} \cdot \vec{b}}{|\vec{b}|}\)​​

3.1.2、点积与\(cos \theta\)

夹角 \(\theta\)​ 与点积大小的关系：

（补充cos的图）

  1. 若 \(\theta = 0°\)​​​​ ​，\(\vec{a} \cdot \vec{b} = |\vec{a}||\vec{b}|\)​​ ，同向共线
  2. 若 \(\theta = 180°\)​ ，\(\vec{a} \cdot \vec{b} = -|\vec{a}||\vec{b}|\)​ ，反向向共线
  3. 若 \(\theta = 90°\)​ ，\(\vec{a} \cdot \vec{b} = 0\)​ ，垂直
  4. 若 \(\theta < 90°\)​​ ，\(\vec{a} \cdot \vec{b} > 0\)​​ ，两向量的夹角为锐角
  5. 若 \(\theta > 90°\) ，\(\vec{a} \cdot \vec{b} < 0\) ，两向量的夹角为钝角

从上面我们可以发现，我们可以通过点积的结果来判断，共线的方向，以及垂直

3.2、叉积（Cross Product）

我们定义向量 \(\vec{a},\vec{b}\)​ 的向量积为一个向量，记为 \(\vec{a} \times \vec{b}\) ，其模与方向定义如下：

  1. \(|\vec{a} \times \vec{b}| = |\vec{a}||\vec{b}|\sin<a,b>\)​​
  2. \(\vec{a} \times \vec{b}\) 与 \(\vec{a},\vec{b}\) 都垂直，且 \(\vec{a} \times \vec{b} ，\vec{a},\vec{b}\)​​​ 都符合右手法则。

我的右手和官方的不太一样，平摊（正面或者反面朝上）右手，大拇指和第一个向量方向相同，四指指向另一个（如果别扭就翻过来）。

手心向上，大于零，手心向下，小于零。我觉得这样方便点。

向量积也叫外积，也叫叉积。

double operator ^ (Point b){
	return x*b.y - y*b.x;
}

3.2.1、求平行四边形面积

然而注意到向量积的模，联想到三角形面积计算公式 \(S = \frac{1}{2}ab\sin{C}\)​ ，我们可以发现向量积的几何意义是： \(|\vec{a} \times \vec{b}|\)​ 是以为 $\vec{a},\vec{b} $​​ 为邻边的平行四边形的面积。

3.2.2、判断A，B向量方向关系

若 \(\vec{a} \times \vec{b} > 0\)​​​ ，\(\vec{b}\)​​​ 在 \(\vec{a}\)​​ 的逆时针方向；

若 \(\vec{a} \times \vec{b} < 0\)​ ，\(\vec{b}\)​ 在 \(\vec{a}\)​​ 的顺时针方向；

若 \(\vec{a} \times \vec{b} = 0\)​ ， \(\vec{b}\)​ 与 \(\vec{a}\)​​​ 可能同向共线，也可能反向共线（刚才的点积就可以判断是同向共线还是反向共线），还可能平行。

3.4、向量的运算

3.4.1、向量模长

向量模长，就是原点到 (x,y) 坐标的距离（在我们的上述的建模方式下），对于向量 \(\vec{a}=(x,y)\) ，\(|\vec{a}|=\sqrt{(x-0)^2+(y-0)^2}\) ，所以我们用与求两点距离的方式来求向量的模长。

double len(){
	return hypot(x,y);
}
/// 长度的平方，便于之后的计算
double len2(){
    return x*x + y*y;
}

3.4.2、向量加、减

上图可以看到 \(\vec{AC}+\vec{AB}\) 的值就是坐标 \(D(C_x+B_x,C_y+B_y)\)

同时可以看到 \(\vec{AB}-\vec{AC}\) 的值就是坐标 \(\vec{CB}(B_x-C_x,B_y-C_y)\)

Point operator + (Point b){
	return Point( x+b.x,y+b.y );
}
Point operator - (Point b){
	return Point( x-b.x,y-b.y );
}

（ 注：这些 operator 都是直接重载的结构体运算符，如果没有这方面的前置知识，请先学习如果重载结构体的运算符 ）

3.4.3、向量数乘（放缩）

对于 \(\vec{a}=(x,y)\) ，\(\lambda \vec{a}=(\lambda x,\lambda y)\) ，除法的话等价于\(\frac{\vec{a}}{\lambda}=\frac{1}{\lambda}\vec{a}=(\frac{1}{\lambda} x,\frac{1}{\lambda} y)\)

Point operator * (double k){
	return Point(k*x,k*y);
}
Point operator / (double k){
    return Point(x/k , y/k);
}

对于数乘的具体应用是，把一个向量变为长度为 \(r\) 的向量，但是向量方向不变

Point trunc(double r){
    double l = len();
    r /= l;
    return Point(x*r,y*r);
}

3.4.4、判断向量相等

因为计算机运行精度问题，所以我们不能直接用 == 来判断 double 类型的数据是否相等，这里引入，eps 来表示题目允许的误差范围，sgn 来代替比较运算符。

const double eps = 1e-8;
int sgn(double x)
{
    if(fabs(x) < eps ) return 0;
    if( x < 0 ) return -1;
    else return 1;
}

1e-8 表示的是\(1\times 10^{-8}\) 是判断两值相等所允许的误差范围，比如 \(1.000000002\) 与 \(1.000000005\) 被视为相等的两个数，均表示 \(1.00000000\) ，如果两数相减的值，小于 eps 那么我们就视为这是两个相等的数，基于此，我们可重载 < ，== 运算符。小于运算符的排序规则是，x 小的，y 小的在前，x作为第一标准，y作为第二标准。

bool operator == (Point b) {
	return sgn(x-b.x) == 0 && sgn(y-b.y) == 0;
}
bool operator < (Point b){
    return sgn(x-b.x)==0?sgn(y-b.y)<0:b.x;
}

3.4.5、两向量夹角

前缀知识：

  1. atan(x)表示求的是 x 的反正切，其返回值为\([-\frac{\pi}{2},+\frac{\pi}{2}]\)之间的一个数。

  2. atan2(y，x)求的是 y/x 的反正切，其返回值为\([-\pi,+\pi]\)之间的一个数。

kuangbin用的方法是

double rad(Point a,Point b){
    Point p = *this;
    return fabs( atan2( fabs( (a-p)^(b-p) )  , (a-p)*(b-p) ) );
}

我看了一眼没有理解，就直接推了一遍，并测试了在精度为 esp=1e-10的情况下，两则的答案相同

我自己的公式推导：

对于两个从原点出发的向量

由点积公式可得： $x_1 y_1+x_2y_2 = |A|\cdot |B|\cdot \cos \theta $​​​​ ，由模长公式和反三角函数可得 \(\theta = \acos({\frac{x_1 y_1+x_2y_2}{\sqrt{x_1^2+y_1^2} \times \sqrt{x_2^2+y_2^2}} })\)​​​​

/// 求 a,b 向量对于原点的夹角
/// 如果是三个数的夹角的话，this 为中间的那个点，a,b为两端
/// return fabs( acos( ((a-tp)*(b-tp))/( (a-tp).len() * (b-tp).len() ) ) )
double rad(Point a,Point b){
    Point tp = *this;
    return fabs(acos( (tp*b)/(tp.len() * b.len()) ));
}

哦哦哦，我明白了，kuangbin的是这样推的，还是推荐用kuangbin的因为代码更加简介，我的推导只用了点积（之后会讲解），或者也可以只用叉积，但是如果我们把叉积与点积 相除，就会得到 kuangbin 的式子，下面为推导。

\(叉积：x_1 y_2-x_2y_1 = |A|\cdot |B|\cdot \sin \theta\)​

\(点积：x_1 y_1+x_2y_2 = |A|\cdot |B|\cdot \cos \theta\)​ ，相除可得

$\frac{\sin \theta }{\cos \theta } = \frac{x_1 y_2-x_2y_1}{x_1 y_1+x_2y_2} $​​，进一步推导 \(\tan \theta = \frac{\vec{a} \cdot \vec{b} }{\vec{a} \times \vec{b} }\)​​，然后求个反三角函数就好了，之后学习了点积和叉积的计算，同时我们重载了运算符之后，就会发现用这个方法来实现求夹角，代码就会比较精简。

kuangbin yyds!

3.4.6、向量旋转

设 \(\vec{a}=(x,y)\)​ ，倾角为 $\theta $​ ，长度为 \(l = \sqrt{x^2+y^2}\)​ 。则 \(x=l\cos{\theta},y=l\sin{\theta}\)​ 。令其逆时针旋转 \(\alpha\)​ ​度角，得到向量 \(\vec{b} = ( l\cos(\theta + \alpha),l\sin(\theta + \alpha) )\)​ 。

由三角恒等变换得，\(b = (l(\cos{\theta}\cos{\alpha}-\sin{\theta}\sin{\alpha}),(\sin{\theta}\sin{\alpha}+\cos{\theta}\sin{\alpha}))\)​

将 \(x,y\) 的值带入式子得到，\(b=(x\cos{\alpha}-y\sin{\alpha},y\cos{\alpha}+x\sin{\alpha})\)

IJNGxe.png

绕着 P 点旋转

 Point rotata(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);
     /// 绕原点旋转后，加上原来P点的偏移量
 }

对于特殊的选择 90°，我们知道 \(\cos{90°}=0\)​ 所以就直接变成了

///逆时针旋转90度
Point rotleft(){
	return Point(y,-x);
}
///顺时针旋转90度
	Point rotright(){
return Point(y,-x);
}

3.5、向量板子

struct Point
{
    double x,y;
    double angle;
    int id;
    Point(){}							        /// 无参构造
    Point(double _x,double _y){
        x = _x; y = _y;
    }   /// 有参构造
    /// 向量加法
    Point operator + (const Point b) const{
        return Point( x+b.x,y+b.y );
    }
    /// 向量乘法
    Point operator - (const Point b) const{
        return Point( x-b.x,y-b.y );
    }
    /// 向量数乘
    Point operator * (const double k) const{
        return Point(k*x,k*y);
    }
    /// 向量数除
    Point operator / (const double k) const{
        return Point(x/k , y/k);
    }
    bool operator == (const Point b) const{
        return sgn(x-b.x) == 0 && sgn(y-b.y) == 0;
    }
    bool operator < (const Point b) const {
        return sgn(x-b.x)==0?sgn(y-b.y)<0:x<b.x;
    }
    /// 点积
    double operator * (const Point b) const{
        return x*b.x + y*b.y;
    }
    /// 叉积
    double operator ^ (const Point b) const{
        return x*b.y - y*b.x;
    }
    /// 两点之间的距离
    double distance(const Point P) const {
        return hypot(x-P.x,y-P.y); /// 库函数，求根号下，两数平方和
    }
    /// 向量长度
    double len(){
        return hypot(x,y);
    }
    /// 长度的平方，便于之后的计算
    double len2(){
        return x*x + y*y;
    }
    /// 化为长度为 r 的向量
    Point trunc(double r){
        double l = len();
        r /= l;
        return Point(x*r,y*r);
    }
    /// 以 p 为中心点，pa,pb的夹角大小
    double rad(Point a,Point b){
        Point p = *this;
        return fabs( atan2( fabs( (a-p)^(b-p) )  , (a-p)*(b-p) ) );
    }
    /// 绕 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);
        /// 绕原点旋转后，加上原来P点的偏移量
    }
    ///逆时针旋转90度
    Point rotleft(){
        return Point(y,-x);
    }
    ///顺时针旋转90度
    Point rotright(){
        return Point(y,-x);
    }
};

4、点与直线

4.1、直线，线段的表示

在计算几何中，直线，线段我们都用两个 Point 来表示，如果把 Point 视为点的话，那么所组成的就是线段，如果把 Point 视为向量的话，那么所组成的就是直线。

struct Line
{
    Point s,e;

    Line(){}
    Line(Point _s,Point _e):s(_s),e(_e){}
};

4.1.1、由解析式来表示直线

  1. 斜截式式 \(y=kx+b\)
  2. 标准方程 \(ax + by + c = 0\)

通常，我们采用斜截式的时候是用直线与 \(x\)​ 轴的顺时针倾斜角与直线上一点来得到两点。

我们设直线的斜倾角为 \(\theta\) ，斜率为 \(k\) ，已知两点 \(P_1(x_1,y_1),P_2(x_2,y_2)\)

\(k = \frac{y1-y2}{x1-x2}\)​ ，\(\tan{\theta}= \frac{y1-y2}{x1-x2}=k\)​​ ，那么\(\theta = atan(k)\)​。

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);
    }
}

4.1.2、线段长度

就是两点之间的距离

double length(){
    return s.distance(e);
}

4.1.3、直线基于 x 轴的斜倾角

就是用 4.1.1 中的公式来求

注：const doubel pi = acos(-1);

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;
}

4.2、点与直线

4.2.1、点与直线的位置关系

IYKbvt.png

在二维平面上，点和直线有 3 种位置关系，点在直线左侧，点在直线右侧，点在直线上。记点为 \(P\)，直线上的点为 \(P_1,P_2\) ，取 \(P_1P\)，\(P_1P_2\) 构成向量用叉积就能得到位置关系。

///点与直线关系
///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;
}

4.2.2、点到直线的距离

IYo7Jx.png

如图，我们设点 \(P_0\)，取直线点 \(P_1,P_2\)，我们用\(P_0-P_1，P_2-P_1\) 来将 \(P_1\) 点移动到原点，我们要求的点到直线的距离的长度设为 \(l\)，\(P_0-P_1，P_2-P_1\) 两个向量的夹角设为 \(\theta\)。

通过叉积公式 \(\vec{a} \times \vec{b} = |\vec{a}||\vec{b}|\sin{\theta}\)​

我们可以得到 \(|\vec{b}|\sin{\theta}= \frac{\vec{a} \times \vec{b}}{|\vec{a}|}=l\)​​，这样我们就得到了点到直线的距离的计算公式，因为叉积是有正负的，所以最后要取个绝对值。

double dispointtoline(Point p){
    return fabs( (p-s)^(e-s) ) / length();
}

4.2.3、点在直线上的投影（垂足）

还是上面哪个图，这次我们是要求 \(|\vec{b}|*\cos|\theta|\)​​ , 我自己想的就是直接用点积来求，但是 kuangbin 和 黑书 都是用的成比例的方式，是因为求 cos 损失精度？需要 fabs ？

...

我去模拟了一下，才发现 lineProg 这个函数的返回值是 Point ！！！所以这波是数学不过关。 我之前一直理解的是标量投影。但是这里我们求的是矢量投影，见百度百科。

IYqIGn.png

设 \(P_0\)​​​ 到 \(P_1P_2\)​​​ 的垂足为 \(P_3\)​​​ （均为向量），为了计算不用cos 所以我们要求一个 \(k = \frac{|P_3-P_1|}{|P_2-P_1|}\)​​​ ，长度比例嘛， \(P_3\)​​​ 就是我们要求的 矢投影。

由点积的公式：$ (P_0-P_1) \cdot (P_2-P_1) = |P_0-P_1||P_2-P_1|\cos \theta $​ ，我们把 \(P_3\) 引入

公式变为 ： \((P_0-P_1) \cdot (P_2-P_1) = |P_2-P_1| * |P3-P_1|\)​​

带入 k 的式子中： $ k = \frac{|P_3-P_1|}{|P_2-P_1|}=\frac{(P_0-P_1) \cdot (P_2-P_1)}{|P_2-P_1| * |P_2-P_1|} $

所以，$P_3 = P_1 + k*(P_2-P_1) $​

/// 返回点p在直线上的投影
Point lineprog(Point p){
	return s + ( ( (e-s)*((e-s)*(p-s)) ) / ( (e-s).len2() ) );
}

这个也想相当于是，求点到直线的垂足向量。

4.2.4、点关于直线的对称点

求一个点对一条直线的对称点。先求点 \(P\) 在直线上的投影点 \(P_3\) （承接上一小节），再求对称点 。

就是把垂足的长度延长一倍就好了。

Point symmetrypoint(Point p){
    Point q = lineprog(p);
    return Point(2*q.x-p.x,2*q.y-p.y);
}

一组测试数据

int main()
{
    Point p = Point(4,5);
    Line l = Line( Point(0,0),Point(1,10) );

    Point te = l.lineprog(p) ;
    printf("%.10f\n%.10f\n",te.x,te.y);
    Point t = l.symmetrypoint(p);
    printf("%.10f\n%.10f",t.x,t.y);


    return 0;
}
/* 输出
0.5346534653
5.3465346535
-2.9306930693
5.6930693069
*/

4.3、点与线段

4.3.1、点与线段的位置关系

判断点 \(p\)​ 是否在线段上， 先用叉积判断是共线，再判断是否与端点组成的线段与原线段形成钝角。

在 3.1.2 节，我们知道：

若 \(\theta = 180°\) ，\(\vec{a} \cdot \vec{b} = -|\vec{a}||\vec{b}|\) ，反向向共线

同时点 \(p\)​ 还有可能在端点上，所以最后为 \(\le 0\)​

bool point_on_seg(Point p){
	return sgn((p-s)^(e-s) ) == 0 && sgn( (p-s)*(p-e) ) <= 0 ;
}

4.3.2、点到线段的距离

点到线段的距离和点到直线的距离有所不同，设点为 \(P\) ，线段\(AB\) 。

在 P到AB的垂线长度，P到A的距离，P到B的距离，三者中取最小值。It7o1P.png

这是 if 判断的情况

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);
}

5、直线与线段

5.1、直线与直线的位置关系

5.1.1、判断直线与直线平行

在3.2.2小节，我们知道：

若 \(\vec{a} \times \vec{b} = 0\) ， \(\vec{b}\) 与 \(\vec{a}\)​ 可能同向共线，也可能反向共线（刚才的点积就可以判断是同向共线还是反向共线），还可能平行。

bool parallel(Line v){
	return sgn( (e-s)^(v.e-v.s) ) == 0 ;
}

5.1.2、两直线位置关系

思路就是，先判断平行，再判断是否共线，如果都不是的话那肯定是相交了。

这里共线我们直接用点与直线的共线来判断。

/// 判断两直线的位置关系
///0 平行
///1 共线
///2 相交
int linecrossline(Line v){
    if( v.parallel(*this) ) return v.relation(s) == 3;
    return 3;
}

5.2、直线与线段的位置关系

(-1)^(1) = 2

/// 直线与线段位置关系 -*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);
}

5.3、线段与线段的位置关系

IUi9yD.png

可用上图数据进行验证

///两线段相交判断
///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);
}

5.4、线段到线段的距离

double disssegtoseg(Line v){
	return min( { dispointtoseg(v.s),dispointtoseg(v.e) ,v.dispointtoseg(s), v.dispointtoseg(e) } );
}

直线到直线的距离，只有当两直线平行的时候才有意义，在两直线平行的情况下，我们直接随便取一个直线的 s,e 然后求点到直线的距离就可以了。

5.5、两条直线的交点

两直线，两线段的求交点是一样的，但是都要先判断是否相交，再求交点

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) );
}

6、线板子

struct Line
{
    Point s,e;
    Line(){}
    Line(Point _s,Point _e){
        s =_s; e= _e;
    }

    /// 点斜
    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);
        }
    }
    /// 线段长度
    double length(){
        return s.distance(e);
    }
    /// 返回 0 <= angle <= pi 的基于 x轴 的斜倾角
    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;
    }
    /// 点到直线的距离
    double dispointtoline(Point p){
        return fabs( (p-s)^(e-s) ) /length();
    }
    /// 返回点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);
    }
    /// 点与线段的位置关系
    bool point_on_seg(Point p){
        return sgn((p-s)^(e-s) ) == 0 && sgn( (p-s)*(p-e) ) <= 0 ;
    }
    /// 点到线段的距离
    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);
    }
    /// 判断平行
    bool parallel(Line v){
        return sgn( (e-s)^(v.e-v.s) ) == 0 ;
    }
    /// 判断两直线的位置关系
    //0 平行
    //1 共线
    //2 相交
    int linecrossline(Line v){
        if( v.parallel(*this) ) return v.relation(s) == 3;
        return 3;
    }
    /// 直线与线段位置关系 -*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);
    }

    ///两线段相交判断
    ///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);
    }
    /// 线段到线段的距离
    double disssegtoseg(Line v){
        return min( min( dispointtoseg(v.s),dispointtoseg(v.e)) , min(v.dispointtoseg(s), v.dispointtoseg(e) ) );
    }

    /// 求两直线的交点
    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) );
    }

};

6、多边形

之前的多数来自于 kuangbin计算几何板子 ， 黑书 ，之后的可能就是我自己想的很多方法了

6.1、多边形的表示

maxp 为点的数量，狂兵的把求凸包，求极角排序都写在 struct 里面了，但是我觉得这个应该是在 点集上进行的操作，所以我写在了外面。

const int maxp = 1e5+10;

struct Polygon
{
    int n;
    Point p[maxp];
    Line l[maxp];
    /// 在多边形中添加点
    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]);
        }
    }
};

6.2、多边形判定与求得

6.2.1、极角排序

对于一个二维平面上的点集，我们选择点集中，最左下角的点，做 J 点到 x轴 的平行线，为直线 JA，然后逆时针旋转 直线JA ，接触其他点的先后循序就是 极角排序 的结果。

IUw4tH.png

如果有两个点与 JA 都相交，那么我们定义与 J 的距离短的优先级高。

IUw7ct.png

6.2.1.1、atan2得到极角排序

在 3.4.5 中我们知道了 atan 与 atan2 的区别，因为 atan2 的取值是在 \([-\pi,+\pi]\) 所以，我们用 atan2 来求，直线之间的夹角，因为极角排序其实就是夹角排序。

之前狂兵的板子中的 Point 是没有angle 属性的。我这里把极角排序从多边形中独立出来，多边形就只处理多边形那些关系，不管求凸包，极角排序这些，因为这些其实都是 点集 在干的事。

我调整了一下点的位置，把他们都变成了整数。

IU2CvQ.png

我先是这样写的 atan2 cmpAtan2Point 为最左下角的 J

te[i].angle = atan2(  te[i].x - cmpAtan2Point.x,te[i].y - cmpAtan2Point.y );

然后我发现我是 x,y 写反了，最后他变成了 J与y轴平行，顺时针的过去！还有一个值得注意的问题：

\(atan2 \in [-\pi,+\pi]\) 所以，刚才的与 y轴平行 的情况中 B 点的值为负，反而排到前面了，所以计算 atan2 的时候小于零的数都 \(+2π\)​ 。​

同时要注意，在 Point中，我们之前重载的 < 是 x小就小 x 的优先级比 y 高，这里我们需要提高 y 的优先级

修改的代码：

/// Point 的重载
bool operator < (Point b) const {
	return sgn(y-b.y)==0?sgn(x-b.x)<0:y<b.y;
}
/// -------------------------------------------------///
/// 用来指定排序的极点
Point cmpAtan2Point = Point(inf,inf);
/// sort cmp
bool cmpAtan2(Point a,Point b)
{
    if( sgn(a.angle - b.angle) == 0 )    return a.distance(cmpAtan2Point) < b.distance(cmpAtan2Point );
    return a.angle < b.angle;
}

Point te[N];
int main()
{
    te[0] = {6,10};
    te[1] = {8,2};
    te[2] = {10,6};
    te[3] = {10,14};
    te[4] = {12,4};
    te[5] = {12,8};
    te[6] = {16,12};
    te[7] = {18,10};
    te[8] = {18,8};
    te[9] = {18,4};
    te[10] = {6,8};	/// 我又加了一个，来验证负数的情况
    /// 取到最小的数，所以需要重载运算符
    for(int i = 0; i < 11 ; i++){
        cmpAtan2Point = min(cmpAtan2Point,te[i]);
    }
    cout << cmpAtan2Point.x << ' ' << cmpAtan2Point.y << endl;
    for(int i = 0; i < 11 ; i++){
        if( te[i] == cmpAtan2Point ){
            te[i].angle = 0;
        }
        else{
            te[i].angle = atan2(  te[i].x - cmpAtan2Point.x,te[i].y - cmpAtan2Point.y );
            if( sgn(te[i].angle) < 0 ) te[i].angle += 2*pi;
        }
    }

    sort(te,te+11,cmpAtan2);

    for(int i = 0 ; i < 11 ; i++){
        cout << te[i].angle <<endl;
        cout <<te[i].x << ' ' << te[i].y <<endl;
    }

    return 0;
}

最后的结果：（注意这是 做J垂直于Y顺时针的极角排序）

IUWmX4.png

修改为 J平行于X轴，逆时针扫

te[i].angle = atan2( te[i].y - cmpAtan2Point.y , te[i].x - cmpAtan2Point.x); ///只用修改这里

image-20211110142422341

6.2.1.2、叉积得到极角排序

这个比较老火，主要是用来排序，以左下角为极点的。

bool cmpCross(const Point a,const Point b)
{
    double d =  sgn((a-cmpPoint)^(b-cmpPoint))  ;
    if( d == 0 ) return sgn( a.distance(cmpPoint) - b.distance(cmpPoint) ) < 0;
    return d > 0;
}

6.2.2、得到凸包

Gradham 扫描法的变种，Andrew 算法

///--------------------- 凸包 --------------------///

/// 求之前的点集
Point convexP[N];
int Sizep;      /// 点集大小
void getconvex(Polygon &convex)
{
    /// 按照 x 从小到大排序，如果 x 相同，按 y 排序。
    sort(convexP,convexP+Sizep);
    convex.n = Sizep;
    for(int i = 0 ; i < min(Sizep,2) ; i ++){
        convex.p[i] = convexP[i];
    }
    /// 特判
    if( convex.n == 2 && ( convex.p[0] == convex.p[1] ) ) convex.n--;
    if( Sizep <= 2 ) return ;
    int &top = convex.n;
    top = 1;
    for(int i = 2; i < Sizep ; i++){
        while( top && sgn( (convex.p[top] - convexP[i])^(convex.p[top-1]-convexP[i]) ) <= 0  ) top--;
        convex.p[++top] = convexP[i];
    }
    int temp = top;
    convex.p[++top] = convexP[Sizep-2];

    for(int i = Sizep-3 ; i >= 0 ; i--){
        while( top != temp && sgn( (convex.p[top] -convexP[i])^(convex.p[top-1]-convexP[i])) <= 0 ) top --;
        convex.p[++top] = convexP[i];
    }
    if( convex.n == 2 && ( convex.p[0] == convex.p[1] ) ) convex.n--;
    /// 这样求出来的顺时针的凸包，如果要逆时针的话，来一次极角排序就好了。
    /// 极角排序
    cmpPoint = convex.p[0];
    sort(convex.p,convex.p+convex.n,cmpCross);
}

6.2.3、判定多边形是否为凸多边形

如果直接用点集来判断的话，要先求一波极角排序

bool isconvex(){
        bool s[2];
        memset(s,0,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;
}

6.3、多边形周长

double getcircumference(){
    double sum = 0;
    for(int i = 0 ; i < n ; i++){
        sum += p[i].distance(p[(i+1)%n]);
    }
    return sum;
}

6.4、多边形面积

如果是凹多边形，要按顺序排列才行。凸多边形，跑一边极角。

double getarea(){
    double sum = 0;
    for(int i = 0; i < n ; i++){
        sum += (p[i]^p[(i+1)%n]);
    }
    return fabs(sum)/2;
}

6.5、多边形的重心

对于面积切割的加权平均。三角形的顶点重量都是分均的，但是多边形就不是了，所以要面积分均才是真正的重心。下图的这个5边形的重心要向上移，而不是在中间。

I0UqAO.png

Point getbarycenter(){
    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;
}

6.6、点与多边形

6.6.1、点与多边形的位置关系

//3 点上
//2 边上
//1 内部
//0 外部
int relationpoint(Point q){
    /// 在点上
    for(int i = 0; i < n; i++){
        if( p[i] == q ) return 3;
    }
    getLine();  /// 之前getline了的话，可以注释节约时间
    /// 在边上
    for(int i = 0; i < n ; i++){
        if( l[i].point_on_seg(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;
}

6.7、直线与多边形

  1. 判断线段 \(AB\) 是否在任意多边形 \(Poly\) 以内：不相交且两端点 \(A,B\) 均在多边形以内。

  2. 判断线段 \(AB\) 是否在凸多边形 \(Poly\) 以内：两端点 \(A,B\)​ 均在多边形以内。

6.8、多边形与多边形

判断任意两个多边形是否相离：属于不同多边形的任意两边都不相交 且 一个多边形上的任意点都不被另一个多边形所包含。

int relationpolygon(Polygon Poly){
    getLine();
    Poly.getLine();
    for(int i1 = 0 ; i1 < n ; i1++){
        int j1 = (i1+1)%n;
        for(int i2 = 0 ; i2 <= Poly.n ; i2++){
            int j2 = (i2+1)%Poly.n;
            Line l1 = Line(p[i1],p[j1]);
            Line l2 = Line(Poly.p[i2],Poly.p[j2]);
            if( l1.segcrossseg(l2) ) return 0;
            if( Poly.relationpoint(p[i1]) || relationpoint(Poly.p[i2]) ) return 0;
        }
    }
    return 1;
}

6.5、多边形与圆

6.6、多边形板子

struct Polygon
{
    int n;
    Point p[maxp];
    Line l[maxp];

    /// 在多边形中添加点
    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]);
        }
    }
    /// 判断多边形是不是凸包
    /// 如果是直接对点集效验的话，要先极角排序
    bool isconvex(){
        bool s[2];
        memset(s,0,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;
    }
    /// 多边形周长
    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;
    }
    /// 重心
    Point getbarycenter(){
        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;
    }
    /// 判断点与任意多边形的关系
    //3 点上
    //2 边上
    //1 内部
    //0 外部
    int relationpoint(Point q){
        /// 在点上
        for(int i = 0; i < n; i++){
            if( p[i] == q ) return 3;
        }
//        getLine();  /// 之前getline了的话，可以注释节约时间
        /// 在边上
        for(int i = 0; i < n ; i++){
            if( l[i].point_on_seg(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;
    }
    /// 判断多边形是否与多边形相离
    int relationpolygon(Polygon Poly){
        getLine();
        Poly.getLine();
        for(int i1 = 0 ; i1 < n ; i1++){
            int j1 = (i1+1)%n;
            for(int i2 = 0 ; i2 <= Poly.n ; i2++){
                int j2 = (i2+1)%Poly.n;
                Line l1 = Line(p[i1],p[j1]);
                Line l2 = Line(Poly.p[i2],Poly.p[j2]);
                if( l1.segcrossseg(l2) ) return 0;
                if( Poly.relationpoint(p[i1]) || relationpoint(Poly.p[i2]) ) return 0;
            }
        }
        return 1;
    }

};

/// 用来指定排序的极点
Point cmpPoint = Point(inf,inf);

/// 需重载 Point 的小于符号
bool cmpAtan2(Point a,Point b)
{
    if( sgn(a.angle - b.angle) == 0 )    return a.distance(cmpPoint) < b.distance(cmpPoint );
    return a.angle < b.angle;
}

/// 慎用，第一个基本上就够了
bool cmpCross(const Point a,const Point b)
{

    double d =  sgn((a-cmpPoint)^(b-cmpPoint))  ;
    if( d == 0 ) return sgn( a.distance(cmpPoint) - b.distance(cmpPoint) ) < 0;
    return d > 0;
}

///--------------------- 凸包 --------------------///

/// 求之前的点集
Point convexP[N];
int Sizep;      /// 点集大小
void getconvex(Polygon &convex)
{
    /// 按照 x 从小到大排序，如果 x 相同，按 y 排序。
    sort(convexP,convexP+Sizep);
    convex.n = Sizep;
    for(int i = 0 ; i < min(Sizep,2) ; i ++){
        convex.p[i] = convexP[i];
    }
    /// 特判
    if( convex.n == 2 && ( convex.p[0] == convex.p[1] ) ) convex.n--;
    if( Sizep <= 2 ) return ;
    int &top = convex.n;
    top = 1;
    for(int i = 2; i < Sizep ; i++){
        while( top && sgn( (convex.p[top] - convexP[i])^(convex.p[top-1]-convexP[i]) ) <= 0  ) top--;
        convex.p[++top] = convexP[i];
    }
    int temp = top;
    convex.p[++top] = convexP[Sizep-2];

    for(int i = Sizep-3 ; i >= 0 ; i--){
        while( top != temp && sgn( (convex.p[top] -convexP[i])^(convex.p[top-1]-convexP[i])) <= 0 ) top --;
        convex.p[++top] = convexP[i];
    }
    if( convex.n == 2 && ( convex.p[0] == convex.p[1] ) ) convex.n--;
    /// 这样求出来的顺时针的凸包，如果要逆时针的话，来一次极角排序就好了。
    /// 极角排序
//    cmpPoint = convex.p[0];
//    sort(convex.p,convex.p+convex.n,cmpCross);
}

7、圆

7.1、圆的表示

struct Circle
{
    Point c;
    double r;
    Circle(){}
    Circle(Point _c,double _r):c(_c),r(_r){}
    Circle(double x,double y,double _r){    /// 点的坐标，圆心
        c = Point(x,y); r = _r;
    }
    bool operator == (Circle v){
        return (p == v.p) && sgn(r-v.r) == 0;
    }   /// 以圆心为主排序，半径为副排序
    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;
    }
};

7.2、点和圆的关系

/// 点和圆的位置关系
//0 圆外
//1 圆上
//2 圆内
int relationPoint(Point b){
    double dst = b.distance(p);
    if( sgn(dst-r)<0 ) return 2;
    if( sgn(dst-r)==0 ) return 1;
    return 0;
}

7.3、直线和圆的关系

7.3.1、直线和圆的位置关系

比较的是圆心到直线的距离和半径的关系

//0 圆外
    //1 圆上
    //2 园内
    int relationLine(Line v){
        double dst = v.dispointtoline(p);
        if( sgn(dst-r) < 0 ) return 2;
        if( sgn(dst-r) == 1 ) return 1;
        return 0;
    }

7.3.2、直线和圆的交点

7.4、线段和圆的关系

比较的是圆心到线段的距离和半径的关系

int relationSeg(Line v){
    double d = p.dispointtoSeg(p);
    if( sgn(dst-r) < 0 ) return 2;
    if( sgn(dst-r) == 0 ) return 1;
    return 0;
}

7.5、圆与圆的关系

7.5.1、圆与圆的位置关系

//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;
}

7.5.2、圆与圆的交点个数

7.5.3、圆与圆的交点

7.5.4、两圆相交面积

7.6、直线和圆的交点

7.7、圆与三角形

7.7.1、三角形内接圆

7.7.2、三角形外接圆

7.7.3、圆与三角形相交面积

7.8、最小圆覆盖

8、专题

8.1、求凸包

6.2.2

8.2、最近点对

8.3、旋转卡壳

IrDFMD.png

参考博客

8.4、最远点对

8.4、半平面交

8.5、最小矩形覆盖

9、测试题

9.1、【LightOJ 1203】凸包+最小顶点夹角

9.2、【2020ICPC 昆明】线段与直线相交判断

这个题可以说是我计算几何的梦魇。今年上半年打昆明的时候，很快就发现了这个题是铜，但是到最后都没有过，今天知道原因了，没有排序的优先级有问题....

题意

给你一个点集，让你求任意一个点，与其他点所组成的直线与一条线段的第 k 个交点。保证没有3点共线。

题解

直接我们打个表，把所有的交点已经坐标都求出来，计算几何怎么搞都行，只要是基础的那些求交点，平行，都是\(O(1)\)​ 的世界复杂度。但是一定要对交点排序，排序的标准是距离线段第一个点的距离的长度！！。

还有就是，计算几何别用 cin ！！！ 计算几何别用 cin ！！！ 计算几何别用 cin ！！！

#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;
}
const int N = 1e3+10;
const int pi = acos(-1);
//square of a double
inline double sqr(double x){return x*x;}


struct Point
{
    double x,y;
    Point(){}							        /// 无参构造
    Point(double _x,double _y):x(_x),y(_y){}    /// 有参构造
    /// 向量加法
    Point operator + (Point b){
        return Point( x+b.x,y+b.y );
    }
    /// 向量乘法
    Point operator - (Point b){
        return Point( x-b.x,y-b.y );
    }
    /// 向量数乘
    Point operator * (double k){
        return Point(k*x,k*y);
    }
    /// 向量数除
    Point operator / (double k){
        return Point(x/k , y/k);
    }
    bool operator == (Point b) {
        return sgn(x-b.x) == 0 && sgn(y-b.y) == 0;
    }
    bool operator < (Point b){
        return sgn(x-b.x)==0?sgn(y-b.y)<0:b.x;
    }
    /// 点积
    double operator * (Point b){
        return x*b.x + y*b.y;
    }
    /// 叉积
    double operator ^ (Point b){
        return x*b.y - y*b.x;
    }
    /// 两点之间的距离
    double distance(Point P){
        return hypot(x-P.x,y-P.y); /// 库函数，求根号下，两数平方和
    }

}P[N];

struct Line
{
    Point s,e;

    Line(){}
    Line(Point _s,Point _e):s(_s),e(_e){}

    /// 直线与线段位置关系 -*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);
    }

    /// 求两直线的交点
    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 t1,t2,t3,t4;
vector<Point> ans[N];
bool cmp(Point a,Point b){
    return a.distance(Point{t1,t2}) < b.distance(Point{t1,t2});
}
int main()
{

    int n,m;
    scanf("%d%d",&n,&m);
    scanf("%lf%lf%lf%lf",&t1,&t2,&t3,&t4);
    Line seg = Line(Point(t1,t2),Point(t3,t4));

    for(int i = 1 ; i <= n ; i++){
        scanf("%lf%lf",&P[i].x,&P[i].y);
    }

    for(int i = 1; i <= n ; i ++){
        for(int j = 1; j <= n ; j++){
            if(i ==j)continue;
            Line tel = Line(P[i],P[j]);

            if( tel.linecrossseg(seg) ){
                Point tep = tel.crosspoint(seg);
                ans[i].push_back(tep);
            }

        }
    }

    for(int i = 1 ; i <= n ; i++){
            sort(ans[i].begin(),ans[i].end(),cmp);
    }

    while(m--){
        int te1,te2;
        scanf("%d%d",&te1,&te2);
        if( ans[te1].size() < te2 ){
            puts("-1");
            continue;
        }
        printf("%.6f %.6f\n",ans[te1][te2-1].x,ans[te1][te2-1].y);


    }


    return 0;
}

9.3、【POJ - 2007】极角排序

#include <iostream>
#include <cmath>
#include <stdio.h>
#include <algorithm>
using namespace std;

const double eps = 1e-8;
const double inf = 1e14;
int sgn(double x)
{
    if(fabs(x) < eps ) return 0;
    if( x < 0 ) return -1;
    else return 1;
}
struct Point
{
    double x,y;
    double angle;
    Point(){}							        /// 无参构造
    Point(double _x,double _y){
        x = _x; y = _y;
    }   /// 有参构造
    /// 向量加法
    Point operator + (Point b) const{
        return Point( x+b.x,y+b.y );
    }
    /// 向量乘法
    Point operator - (Point b) const{
        return Point( x-b.x,y-b.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;
    }
    /// 叉积
    double operator ^ (Point b) const{
        return x*b.y - y*b.x;
    }


}P;

/// 用来指定排序的极点
Point cmpPoint = Point(inf,inf);

bool cmpCross(const Point &a,const Point &b)
{
    double d =  (a-cmpPoint)^(b-cmpPoint)  ;
//    if( sgn(d) == 0 ) return sgn( a.distance(cmpPoint) - b.distance(cmpPoint) ) < 0;
    return d >= 0;
}

Point te[100];

int main()
{
    int it = 0;
    while(scanf("%lf%lf",&te[it].x,&te[it].y) != EOF) ++it;

    cmpPoint = Point(0,0);

    sort(te+1,te+it,cmpCross);

    for(int i = 0 ; i < it ; i++){

        printf("(%.0f,%.0f)\n",te[i].x,te[i].y);
    }

    return 0;
}

/*
80 20
50 -60
0 0
70 -50
60 30
-30 -50
90 -20
-30 -40
-10 -60
90 10
*/

9.4、【CF 598C】极角排序

9.5、【LightOj 1203】 凸包 + 边夹角

#include <iostream>
#include <cmath>
#include <stdio.h>
#include <algorithm>
#include <vector>
using namespace std;


const double eps = 1e-18;
int sgn(double x)
{
    if(fabs(x) < eps ) return 0;
    if( x < 0 ) return -1;
    else return 1;
}
const int N = 1e5+10;
const int maxp = 1e5+10;
const double pi = acos(-1);
const double inf = 1e14;
//square of a double
inline double sqr(double x){return x*x;}

/*
 * Point ()
 * distance(Point P)    - distance from this to P
 * len()                - get the length from (0,0) to (x,y)
 * trunc(double r)      - transform to the vector which length is r
 * rad(Point a)
 */


struct Point
{
    double x,y;
    double angle;
    int id;
    Point(){}							        /// 无参构造
    Point(double _x,double _y){
        x = _x; y = _y;
    }   /// 有参构造
    Point operator + (const Point b) const{return Point( x+b.x,y+b.y );}
    Point operator - (const Point b) const{return Point( x-b.x,y-b.y );}
    /// 向量数乘
    Point operator * (const double k) const{
        return Point(k*x,k*y);
    }
    /// 向量数除
    Point operator / (const double k) const{
        return Point(x/k , y/k);
    }
    bool operator == (const Point b) const{
        return sgn(x-b.x) == 0 && sgn(y-b.y) == 0;
    }
    bool operator < (const Point b) const {
        return sgn(x-b.x)==0?sgn(y-b.y)<0:x<b.x;
    }
    /// 点积
    double operator * (const Point b) const{
        return x*b.x + y*b.y;
    }
    /// 叉积
    double operator ^ (const Point b) const{
        return x*b.y - y*b.x;
    }
    /// 两点之间的距离
    double distance(const Point P) const {
        return hypot(x-P.x,y-P.y); /// 库函数，求根号下，两数平方和
    }
    /// 向量长度
    double len(){
        return hypot(x,y);
    }
    /// 以 p 为中心点，pa,pb的夹角大小
    double rad(Point a,Point b){
        Point p = *this;
        return fabs( atan2( fabs( (a-p)^(b-p) )  , (a-p)*(b-p) ) );
    }
};


struct Polygon
{
    int n;
    Point p[maxp];
};

/// 用来指定排序的极点
Point cmpPoint = Point(inf,inf);

/// 慎用，第一个基本上就够了
bool cmpCross(const Point a,const Point b)
{

    double d =  sgn((a-cmpPoint)^(b-cmpPoint))  ;
    if( d == 0 ) return sgn( a.distance(cmpPoint) - b.distance(cmpPoint) ) < 0;
    return d > 0;
}

///--------------------- 凸包 --------------------///

/// 求之前的点集
Point convexP[N];
int Sizep;      /// 点集大小
void getconvex(Polygon &convex)
{
    /// 按照 x 从小到大排序，如果 x 相同，按 y 排序。
    sort(convexP,convexP+Sizep);
    convex.n = Sizep;
    for(int i = 0 ; i < min(Sizep,2) ; i ++){
        convex.p[i] = convexP[i];
    }
    /// 特判
    if( convex.n == 2 && ( convex.p[0] == convex.p[1] ) ) convex.n--;
    if( Sizep <= 2 ) return ;
    int &top = convex.n;
    top = 1;
    for(int i = 2; i < Sizep ; i++){
        while( top && sgn( (convex.p[top] - convexP[i])^(convex.p[top-1]-convexP[i]) ) <= 0  ) top--;
        convex.p[++top] = convexP[i];
    }
    int temp = top;
    convex.p[++top] = convexP[Sizep-2];

    for(int i = Sizep-3 ; i >= 0 ; i--){
        while( top != temp && sgn( (convex.p[top] -convexP[i])^(convex.p[top-1]-convexP[i])) <= 0 ) top --;
        convex.p[++top] = convexP[i];
    }
    if( convex.n == 2 && ( convex.p[0] == convex.p[1] ) ) convex.n--;

}

Polygon poly;

void solve()
{
    int n;
    cin >> n;
    Sizep = n;

    for(int i = 0; i < n ; i++){
        double tex,tey;
        scanf("%lf%lf",&tex,&tey);
        convexP[i] = Point(tex,tey);
    }
    if( n <= 2 ) {puts("0.000000");return ;}
    getconvex(poly);
    double ans_angle = 1e5;
    for(int i = 0;  i < poly.n ; i++){
        int j = (i+1)%poly.n;
        int k = (j+1)%poly.n;
//        printf("%.2f %.2f %.2f\n",poly.p[j].x,poly.p[i].x,poly.p[k].x);
        double te_angle = poly.p[j].rad(poly.p[i],poly.p[k]);
        ans_angle = min(te_angle,ans_angle);
    }
    printf("%.6f\n",((ans_angle)*180.0)/pi) ;
    return ;
}

int main()
{
    int t;
    scanf("%d",&t);
    for(int i = 1; i <= t ; i++){
        printf("Case %d: ",i);
        solve();
    }

}

/*
2
11
6 10
8 2
10 6
10 14
12 4
12 8
16 12
18 10
18 8
18 4
6 8
4
0 0
10 0
10 10
2 1

1
3
0 0
3 0
0 3
*/

10、kuangbin 计算几何板子

const double eps = 1e-8;
const double pi = acos( -1.0);

///Compares a double to zero
int sgn(double x)
{
    if( fabs(x) < eps ) return 0;
    if( x < 0 ) return -1;
    else return 1;
}
///square of a double
inline double sqr(double x) { return x * x; }


/////////////////////////////////////////////////
struct Point
{
    double x,y;
    Point(){}               ///no arguments constructor
    Point(double _x,double _y) {
        x = _x , y = _y;    ///arguments constructor
    }
    /*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);
    }
    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);
    }
    ///叉积
    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(){
        return hypot(x,y);  ///<cmath>
    }
    ///长度的平方
    double len2(){
        return x * x + y * y;
    }
    ///返回两点的距离
    double distance(Point p){
        return hypot( x - p.x , y - p.y );
    }
    ///计算 pa 和 pb 的夹角
    double rad(Point a,Point b){
        Point p = *this;
        return fabs( atan2( fabs( (a-p)^(b-p) )  , (a-p)*(b-p) ) );
    }
    ///化为长度为r的向量
    Point trunc(double r){
        double l = len();
        if( !sgn(l) ) return *this;
        r /= l;
        return Point(-y,x);
    }
    ///逆时针旋转90度
    Point rotleft(){
        return Point(y,-x);
    }
    ///顺时针旋转90度
    Point rotright(){
        return Point(y,-x);
    }
    ///绕着p点逆时针
    Point rotata(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);
    }
};

struct Line
{
    Point s,e;
    Line(){}
    Line( Point _s, Point _e ){ s =_s ; e=_e; }
    ///由斜倾角angle与任意直线一点确定直线 y = kx + b;
    void input( Point _s, Point _e ){ s =_s ; e=_e; }

    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);
        }
    }
    double length(){ return s.distance(e);}
    ///直线与线段相交判断
    ///-*this line -v seg
    ///2规范相交，1非规范相交,0不相交
    bool linecrossseg(Line v){
        return sgn( (v.s - e) ^ (s - e) ) * sgn(( v.e-e ) ^ (s -e) ) <= 0;
    }
		///点与直线关系
    ///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 point_on_seg(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;
    }
    ///两直线关系 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 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);
    }
    /// 点到直线的距离
    double dispointtoline(Point p){
        return fabs( (p-s)^(e-s) ) / length();
    }

    /// 返回点p在直线上的投影
    Point lineprog(Point p){
        return s + ( ( (e-s)*((e-s)*(p-s)) ) / ( (e-s).len2() ) );
    }
    ///两线段相交判断
    ///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);
    }

};

struct triangle
{
    Point A,B,C;
    Line a,b,c;

    triangle(){}
    triangle(Point _A,Point _B,Point _C){ A = _A ; B = _B ; C = _C;}

    ///求重心
    Point incenter(){
        return Point( ( A.x + B.x + C.x ) / 3, ( A.y + B.y + C.y ) / 3);
    }

};

///已知三点求圆心与半径模板

void cal(int a,int b,int c)//求外心 ，外心为三角形三边的垂直平分线交点，
{
    double a1 = p[b].x - p[a].x, b1 = p[b].y - p[a].y, c1 = (a1*a1 + b1*b1)/2;
    double a2 = p[c].x - p[a].x, b2 = p[c].y - p[a].y, c2 = (a2*a2 + b2*b2)/2;
    double d = a1 * b2 - a2 * b1;
    x = p[a].x + (c1*b2 - c2*b1)/d,y = p[a].y + (a1*c2 - a2*c1)/d;
    r  = dis2(a);
}


struct circle{
    Point p;    ///圆心
    double r;   ///半径
    circle(){}
    circle( Point _p,double _r ) { p = _p ; r = _r; }

    bool operator == (circle v){
        return (p == v.p) && sgn(r - v.r) == 0;
    }
    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;
    }
    ///线段与园的关系
    ///比较的是圆心到线段的距离和半径的的关系
    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;
    }

    /// 直线和圆的关系
    /// 比较的是圆心到直线的距离和半径的关系
    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;
    }

    /// 求直线和圆的交点个数
    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;
    }

    /// 求圆和三角形 pab 的相交面积
    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;
    }

};

const int maxp = 1100;
const int maxl = 2200;

struct polygon
{
    int n;      ///点的数量
    Point p[maxp];
    Line l[maxl];


    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;
        }
    };
    ///极角排序
    ///mi为最左下角的点
    void norm(){
        Point mi = p[0];
        for(int i = 1; i < n; i ++) mi = min(mi,p[i]);
        sort(p, p + n, cmp(mi) );
    }
    /// 判断任意点与多边形的关系
    /// 3在顶点上
    /// 2在边上
    /// 1在内部
    /// 0在外面
    int relationpoint(Point tep)
    {
        for(int i = 0 ; i < n ; i++){
            if( p[i] == tep ) return 3;
        }
        for(int i = 0 ; i < n; i++){
            if( l[i].point_on_seg(tep) ) return 2;
        }
        int tecnt = 0;
        for(int i = 0 ; i < n ; i++)
        {
            int j = (i + 1) % n;
            int c = sgn( (tep - p[j]) ^ (p[i] - p[j]) );
            int u = sgn( p[i].y - tep.y );
            int v = sgn( p[j].y - tep.y );
            if( c > 0 && u < 0 && v >=0 ) tecnt ++;
            if( c < 0 && u >= 0 && v < 0 ) tecnt --;
        }
        return tecnt != 0;
    }

    /// 得到凸包
    /// 得到的凸包里的点编号是 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();///得到的是顺时针的点,排序后逆时针
    }

    ///判断是不是凸多边形,用点集，不是得到凸包之后的多边形
    bool isconvex(){
        bool s[2];
        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;
    }

    ///得到周长
    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;
    }
    ///得到重心
    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){
        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;
    }

};

11、Hoppz 板子

#include <iostream>
#include <cmath>
#include <stdio.h>
#include <algorithm>
#include <vector>
#include <cstring>
using namespace std;


const double eps = 1e-18;
int sgn(double x)
{
    if(fabs(x) < eps ) return 0;
    if( x < 0 ) return -1;
    else return 1;
}
const int N = 1e5+10;
const int maxp = 10;
const double pi = acos(-1);
const double inf = 1e14;
//square of a double
inline double sqr(double x){return x*x;}

/*
 * Point ()
 * distance(Point P)    - distance from this to P
 * len()                - get the length from (0,0) to (x,y)
 * trunc(double r)      - transform to the vector which length is r
 * rad(Point a)
 */


struct Point
{
    double x,y;
    double angle;
    int id;
    Point(){}							        /// 无参构造
    Point(double _x,double _y){
        x = _x; y = _y;
    }   /// 有参构造
    /// 向量加法
    Point operator + (const Point b) const{
        return Point( x+b.x,y+b.y );
    }
    /// 向量乘法
    Point operator - (const Point b) const{
        return Point( x-b.x,y-b.y );
    }
    /// 向量数乘
    Point operator * (const double k) const{
        return Point(k*x,k*y);
    }
    /// 向量数除
    Point operator / (const double k) const{
        return Point(x/k , y/k);
    }
    bool operator == (const Point b) const{
        return sgn(x-b.x) == 0 && sgn(y-b.y) == 0;
    }
    bool operator < (const Point b) const {
        return sgn(x-b.x)==0?sgn(y-b.y)<0:x<b.x;
    }
    /// 点积
    double operator * (const Point b) const{
        return x*b.x + y*b.y;
    }
    /// 叉积
    double operator ^ (const Point b) const{
        return x*b.y - y*b.x;
    }
    /// 两点之间的距离
    double distance(const Point P) const {
        return hypot(x-P.x,y-P.y); /// 库函数，求根号下，两数平方和
    }
    /// 向量长度
    double len(){
        return hypot(x,y);
    }
    /// 长度的平方，便于之后的计算
    double len2(){
        return x*x + y*y;
    }
    /// 化为长度为 r 的向量
    Point trunc(double r){
        double l = len();
        r /= l;
        return Point(x*r,y*r);
    }
    /// 以 p 为中心点，pa,pb的夹角大小
    double rad(Point a,Point b){
        Point p = *this;
        return fabs( atan2( fabs( (a-p)^(b-p) )  , (a-p)*(b-p) ) );
    }
    /// 绕 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);
        /// 绕原点旋转后，加上原来P点的偏移量
    }
    ///逆时针旋转90度
    Point rotleft(){
        return Point(y,-x);
    }
    ///顺时针旋转90度
    Point rotright(){
        return Point(y,-x);
    }
};

struct Line
{
    Point s,e;
    Line(){}
    Line(Point _s,Point _e){
        s =_s; e= _e;
    }

    /// 点斜
    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);
        }
    }
    /// 线段长度
    double length(){
        return s.distance(e);
    }
    /// 返回 0 <= angle <= pi 的基于 x轴 的斜倾角
    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;
    }
    /// 点到直线的距离
    double dispointtoline(Point p){
        return fabs( (p-s)^(e-s) ) /length();
    }
    /// 返回点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);
    }
    /// 点与线段的位置关系
    bool point_on_seg(Point p){
        return sgn((p-s)^(e-s) ) == 0 && sgn( (p-s)*(p-e) ) <= 0 ;
    }
    /// 点到线段的距离
    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);
    }
    /// 判断平行
    bool parallel(Line v){
        return sgn( (e-s)^(v.e-v.s) ) == 0 ;
    }
    /// 判断两直线的位置关系
    //0 平行
    //1 共线
    //2 相交
    int linecrossline(Line v){
        if( v.parallel(*this) ) return v.relation(s) == 3;
        return 3;
    }
    /// 直线与线段位置关系 -*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);
    }

    ///两线段相交判断
    ///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);
    }
    /// 线段到线段的距离
    double disssegtoseg(Line v){
        return min( min( dispointtoseg(v.s),dispointtoseg(v.e)) , min(v.dispointtoseg(s), v.dispointtoseg(e) ) );
    }

    /// 求两直线的交点
    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) );
    }
};

struct Polygon
{
    int n;
    Point p[maxp];
    Line l[maxp];

    /// 在多边形中添加点
    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]);
        }
    }
    /// 判断多边形是不是凸包
    /// 如果是直接对点集效验的话，要先极角排序
    bool isconvex(){
        bool s[2];
        memset(s,0,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;
    }
    /// 多边形周长
    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;
    }
    /// 重心
    Point getbarycenter(){
        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;
    }
    /// 判断点与任意多边形的关系
    //3 点上
    //2 边上
    //1 内部
    //0 外部
    int relationpoint(Point q){
        /// 在点上
        for(int i = 0; i < n; i++){
            if( p[i] == q ) return 3;
        }
//        getLine();  /// 之前getline了的话，可以注释节约时间
        /// 在边上
        for(int i = 0; i < n ; i++){
            if( l[i].point_on_seg(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;
    }
    /// 判断多边形是否与多边形相离
    int relationpolygon(Polygon Poly){
        getLine();
        Poly.getLine();
        for(int i1 = 0 ; i1 < n ; i1++){
            int j1 = (i1+1)%n;
            for(int i2 = 0 ; i2 <= Poly.n ; i2++){
                int j2 = (i2+1)%Poly.n;
                Line l1 = Line(p[i1],p[j1]);
                Line l2 = Line(Poly.p[i2],Poly.p[j2]);
                if( l1.segcrossseg(l2) ) return 0;
                if( Poly.relationpoint(p[i1]) || relationpoint(Poly.p[i2]) ) return 0;
            }
        }
        return 1;
    }

};

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;
    }

    bool operator == (Circle v){
        return (p == v.p) && sgn(r-v.r) == 0;
    }   /// 以圆心为主排序，半径为副排序
    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 relationPoint(Point b){
        double dst = b.distance(p);
        if( sgn(dst-r)<0 ) return 2;
        if( sgn(dst-r)==0 ) return 1;
        return 0;
    }
    /// 直线和圆的关系
    //0 相离
    //1 相切
    //2 相交
    int relationLine(Line v){
        double dst = v.dispointtoline(p);
        if( sgn(dst-r) < 0 ) return 2;
        if( sgn(dst-r) == 0 ) return 1;
        return 0;
    }
    /// 线段和圆的关系
    int relationSeg(Line v){
        double dst = v.dispointtoseg(p);
        if( sgn(dst-r) < 0 ) return 2;
        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;
    }

};


/// 用来指定排序的极点
Point cmpPoint = Point(inf,inf);

/// 需重载 Point 的小于符号
bool cmpAtan2(Point a,Point b)
{
    if( sgn(a.angle - b.angle) == 0 )    return a.distance(cmpPoint) < b.distance(cmpPoint );
    return a.angle < b.angle;
}

/// 慎用，第一个基本上就够了
bool cmpCross(const Point a,const Point b)
{

    double d =  sgn((a-cmpPoint)^(b-cmpPoint))  ;
    if( d == 0 ) return sgn( a.distance(cmpPoint) - b.distance(cmpPoint) ) < 0;
    return d > 0;
}

///--------------------- 凸包 --------------------///

/// 求之前的点集
Point convexP[N];
int Sizep;      /// 点集大小
void getconvex(Polygon &convex)
{
    /// 按照 x 从小到大排序，如果 x 相同，按 y 排序。
    sort(convexP,convexP+Sizep);
    convex.n = Sizep;
    for(int i = 0 ; i < min(Sizep,2) ; i ++){
        convex.p[i] = convexP[i];
    }
    /// 特判
    if( convex.n == 2 && ( convex.p[0] == convex.p[1] ) ) convex.n--;
    if( Sizep <= 2 ) return ;
    int &top = convex.n;
    top = 1;
    for(int i = 2; i < Sizep ; i++){
        while( top && sgn( (convex.p[top] - convexP[i])^(convex.p[top-1]-convexP[i]) ) <= 0  ) top--;
        convex.p[++top] = convexP[i];
    }
    int temp = top;
    convex.p[++top] = convexP[Sizep-2];

    for(int i = Sizep-3 ; i >= 0 ; i--){
        while( top != temp && sgn( (convex.p[top] -convexP[i])^(convex.p[top-1]-convexP[i])) <= 0 ) top --;
        convex.p[++top] = convexP[i];
    }
    if( convex.n == 2 && ( convex.p[0] == convex.p[1] ) ) convex.n--;
    /// 这样求出来的顺时针的凸包，如果要逆时针的话，来一次极角排序就好了。
    /// 极角排序
//    cmpPoint = convex.p[0];
//    sort(convex.p,convex.p+convex.n,cmpCross);
}
【推广】 免费学中医，健康全家人
原文地址：https://www.cnblogs.com/hoppz/p/15619756.html