zoukankan      html  css  js  c++  java
  • 【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 中我们知道了 atanatan2 的区别,因为 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);
    }
    
    
  • 相关阅读:
    P4281 [AHOI2008]紧急集合 / 聚会
    P2622 关灯问题II
    CF1092F Tree with Maximum Cost
    10.28记
    威尔逊定理及证明
    CF895C Square Subsets
    洛谷 P5556 圣剑护符
    Multi-nim结论及证明
    AT2667 [AGC017D] Game on Tree
    洛谷 P4643 [国家集训队]阿狸和桃子的游戏
  • 原文地址:https://www.cnblogs.com/hoppz/p/15619756.html
Copyright © 2011-2022 走看看