精华内容
下载资源
问答
  • 计算几何板子

    2018-05-22 19:48:00
    /************************************** * Au: Hany01 * Date: Feb 8th, 2018 * Prob: Geometry Template ...**************************************/ #include<bits/stdc++.h&...
    /**************************************
     * Au: Hany01
     * Date: Feb 8th, 2018
     * Prob: Geometry Template
     * Email: hany01@foxmail.com
    **************************************/
    
    #include<bits/stdc++.h>
    
    using namespace std;
    
    typedef long long LL;
    typedef pair<int, int> PII;
    #define rep(i, j) for (register int i = 0, i##_end_ = (j); i < i##_end_; ++ i)
    #define For(i, j, k) for (register int i = (j), i##_end_ = (k); i <= i##_end_; ++ i)
    #define Fordown(i, j, k) for (register int i = (j), i##_end_ = (k); i >= i##_end_; -- i)
    #define Set(a, b) memset(a, b, sizeof(a))
    #define Cpy(a, b) memcpy(a, b, sizeof(a))
    #define fir first
    #define sec second
    #define pb(a) push_back(a)
    #define mp(a, b) make_pair(a, b)
    #define ALL(a) (a).begin(), (a).end()
    #define SZ(a) ((int)(a).size())
    #define INF (0x3f3f3f3f)
    #define INF1 (2139062143)
    #define Mod (1000000007)
    #define debug(...) fprintf(stderr, __VA_ARGS__)
    
    template <typename T> inline bool chkmax(T &a, T b) { return a < b ? a = b, 1 : 0; }
    template <typename T> inline bool chkmin(T &a, T b) { return b < a ? a = b, 1 : 0; }
    
    inline int read()
    {
        register int _, __; register char c_;
        for (_ = 0, __ = 1, c_ = getchar(); c_ < '0' || c_ > '9'; c_ = getchar()) if (c_ == '-') __ = -1;
        for ( ; c_ >= '0' && c_ <= '9'; c_ = getchar()) _ = (_ << 1) + (_ << 3) + (c_ ^ 48);
        return _ * __;
    }
    
    const double eps = 1e-10, Pi = acos(-1.0);
    
    struct Point
    {
        double x, y;
        Point(double x = 0, double y = 0): x(x), y(y) {}
    };
    typedef Point Vector;
    
    struct Circle
    {
        Point O; double r;
        Circle(Point O, double r): O(O), r(r) {}
        Point point(double a) { return Point(O.x + cos(a) * r, O.y + sin(a) * r); }
    };
    
    int dcmp(double x) { if (fabs(x) < eps) return 0; else return x < 0 ? -1 : 1; }
    Vector operator + (Vector A, Vector B) { return Vector(A.x + B.x, A.y + B.y); }
    Vector operator - (Vector A, Vector B) { return Vector(A.x - B.x, A.y - B.y); }
    Vector operator * (Vector A, double p) { return Vector(A.x * p, A.y * p); }
    Vector operator / (Vector A, double p) { return Vector(A.x / p, A.y / p); }
    bool operator < (const Point& a, const Point& b) { return a.x < b.x || (a.x == b.x && a.y < b.y); }
    bool operator == (const Point& a, const Point& b) { return !dcmp(a.x - b.x) && !dcmp(a.y - b.y); }
    double Dot(Vector A, Vector B) { return A.x * B.x + A.y * B.y; }
    double Length(Vector A) { return sqrt(Dot(A, A)); }
    double Angle(Vector A, Vector B) { return acos(Dot(A, B) / Length(A) / Length(B)); }
    double Angle(Vector A) { return atan2(A.y, A.x); }
    double Cross(Vector A, Vector B) { return A.x * B.y - A.y * B.x; } //B在A的左边为正、右边为负
    double Area2(Point A, Point B, Point C) { return Cross(B - A, C - A); }
    Vector Rotate(Vector A, double rad) { return Vector(A.x * cos(rad) - A.y * sin(rad), A.x * sin(rad) + A.y * cos(rad)); }
    Vector Normal(Vector A) { double L = Length(A); return Vector(-A.y / L, A.x / L); } //向左旋转90度且长度变成1
    Point GetLineProjection(Point P, Point A, Point B) {
        Vector v = B - A;
        return A + v * (Dot(v, P - A) / Dot(v, v));
    }
    Point GetLineIntersection(Point P, Vector v, Point Q, Vector w) {
        Vector u = P - Q;
        double t = Cross(w, u) / Cross(v, w);
        return P + v * t;
    }
    double DistanceToLine(Point P, Point A, Point B) {
        Vector v1 = B - A, v2 = P - A;
        return fabs(Cross(v1, v2)) / Length(v1);
    }
    double DistanceToSegment(Point P, Point A, Point B) {
        if (A == B) return Length(P - A);
        Vector v1 = B - A, v2 = P - A, v3 = P - B;
        if (dcmp(Dot(v1, v2)) < 0) return Length(v2);
        if (dcmp(Dot(v1, v3)) > 0) return Length(v3);
        return fabs(Cross(v1, v2)) / Length(v1);
    }
    int GetLineCircleIntersection(Point A, Point B, Circle C, Point& S1, Point& S2) {
        double d = DistanceToLine(C.O, A, B);
        register int tmp = dcmp(d - C.r);
        if (tmp > 0) return 0;
        Point P = GetLineProjection(C.O, A, B);
        if (!tmp) { S1 = S2 = P; return 1; }
        double L = sqrt(C.r * C.r - d * d);
        Vector v = (B - A) / Length(B - A);
        S1 = P - v * L, S2 = P + v * L;
        return 2;
    }
    int GetCircleIntersection(Circle C1, Circle C2, Point& P1, Point& P2) {
        double d = Length(C1.O - C2.O);
        if (!dcmp(d)) { if (!dcmp(C1.r - C2.r)) return -1; return 0; }
        if (dcmp(C1.r + C2.r - d) < 0 || dcmp(fabs(C1.r - C2.r) - d) > 0) return 0;
        double a = Angle(C2.O - C1.O), da = acos((C1.r * C1.r + d * d - C2.r * C2.r) / (2 * C1.r * d));
        P1 = C1.point(a - da), P2 = C1.point(a + da);
        if (P1 == P2) return 1;
        return 2;
    }
    int GetTangents(Point P, Circle C, Point& P1, Point& P2) {
        double d = Length(C.O - P);
        if (dcmp(d - C.r) < 0) return 0;
        double aPC = Angle(P - C.O), theta = acos(C.r / d);
        P1 = C.point(aPC + theta), P2 = C.point(aPC - theta);
        if (P1 == P2) return 1;
        return 2;
    }
    int GetTangents(Circle C1, Circle C2, Point* a, Point* b) {
        int cnt = 0;
        if (C1.r < C2.r) swap(C1, C2);
        double d = Length(C1.O - C2.O), dif = C1.r - C2.r, sum = C1.r + C2.r;
        if (d < dif) return 0;//内含
        double base = Angle(C2.O - C1.O);
        if (!dcmp(d) && !dcmp(C1.r - C2.r)) return -1;//重合
        if (!dcmp(d - dif)) {//内切
            a[++ cnt] = C1.point(base), b[cnt] = a[cnt];
            return 1;
        }
        //外公切线
        double ang = acos((C1.r - C2.r) / d);
        a[++ cnt] = C1.point(base + ang), b[cnt] = C2.point(base + ang);
        a[++ cnt] = C1.point(base - ang), b[cnt] = C2.point(base - ang);
        if (!dcmp(d - sum)) a[++ cnt] = C1.point(base), b[cnt] = a[cnt];//一条内公切线
        else if (dcmp(d - sum) > 0) {
            ang = acos((C1.r + C2.r) / d);
            a[++ cnt] = C1.point(base + ang), b[cnt] = C2.point(Pi + base + ang);
            a[++ cnt] = C1.point(base - ang), b[cnt] = C2.point(Pi + base - ang);
        }
        return cnt;
    }
    bool IsPointOnSegment(Point P, Point A, Point B) {
        if (A.x > B.x) swap(A, B);
        if (P.x < A.x || P.x > B.x || dcmp((B.y - A.y) / (B.x - A.x) - (B.y - P.y) / (B.x - P.x))) return 0;
        return 1;
    }
    int IsPointInPolygon(Point P, Point* p, int n) {
        int wn = 0;
        Point A, B;
        For(i, 1, n) {
            A = p[i], B = p[i % n + 1];
            if (!IsPointOnSegment(P, A, B)) return -1;
            int k = dcmp(Cross(B - A, P - A)), d1 = dcmp(A.y - P.y), d2 = dcmp(B.y - P.y);
            if (k > 0 && d1 <= 0 && d2 > 0) ++ wn;
            if (k < 0 && d2 <= 0 && d1 > 0) -- wn;
        }
        if (wn) return 1;
        return 0;
    }
    
    int main()
    {
        return 0;
    }
    展开全文
  • 吉如一几何板子

    2019-02-01 18:15:12
    2018 world final 金牌 吉老师几何板子 struct point{ db x,y; point operator + (const point &k1;) const{return (point){k1.x+x,k1.y+y};} point operator - (const point &k1;) const{return (point){x-k1.x...
  • #include <iostream> #include <stdio.h> #include <string.h> #include <algorithm> #include <queue> #include <map> #include <vector>...set>...
    #include <iostream>
    #include <stdio.h>
    #include <string.h>
    #include <algorithm>
    #include <queue>
    #include <map>
    #include <vector>
    #include <set>
    #include <string>
    #include <math.h>
    #define ll long long
    using namespace std;
    const int maxn=1e6+10;
    const double pi=acos(-1.0);
    const double D_MAX=1e100;
    const double D_MIN=-1e100;
    const double eps=1e-9;
    int sgn(double x){ if(fabs(x) < eps)return 0;  if(x >0) return 1; return -1; }
    int dcmp(double x, double y){ if(fabs(x - y) < eps) return 0; if(x > y) return 1;return -1;}
    void usehanshu(){double x;}//floor(x)向下取整函数ceil(x)向上取整函数round(x)四舍五入函数
    struct Point  { double x,y; Point(double x=0,double y=0) { x=x;y=y; };   };
    struct Segment{ Point a,b;  Segment(Point x,Point y)     { a=x;b=y; };   };
    struct Line   { Point a,b;  Line(Point x,Point y)        { a=x;b=y; };   };
    typedef Point Vector;
    Vector operator + (Vector A, Vector B){ return Vector(A.x+B.x, A.y+B.y); } // 向量相加
    Vector operator - (Point  A, Point  B){ return Vector(A.x-B.x, A.y-B.y); } // 向量生成
    double operator * (Vector A, Vector B){ return A.x*B.x-A.y*B.y;          } // 点积
    double operator ^ (Vector A, Vector B){ return A.x*B.y-A.y*B.x;          } // 叉积
    double Dot(Vector A, Vector B)   { return A.x*B.x + A.y*B.y; }  // 点积
    double Cross(Vector A, Vector B) { return A.x*B.y-A.y*B.x;   }  // 叉积
    double Length(Vector A) { return sqrt(Dot(A, A));}           // 向量长度
    double Angle(Vector A, Vector B){ return acos(Dot(A, B)/Length(A)/Length(B));}  // 角度
    double Area2(Point A, Point B, Point C) { return Cross(B-A, C-A); }  // 四边形面积
    Vector Rotate(Vector A, double rad){ return Vector(A.x*cos(rad)-A.y*sin(rad), A.x*sin(rad)+A.y*cos(rad));}//rad为弧度 且为逆时针旋转的角
    Vector Normal(Vector A)            {double L = Length(A);return Vector(-A.y/L, A.x/L);}//向量A左转90°的单位法向量
    bool ToLeftTest(Point a, Point b, Point c){return Cross(b - a, c - b) > 0;}
    bool xx(Line A,point a){if( sgn( Cross(A.a-a,A.b-a) )==0 ) return 1;return 0;}  // 直线和点
    void XX(Line A,Line B)  // 直线和直线的情况
    {
        Point a=A.a; Point b=A.b;Point c=B.a; Point d=B.b;
        double A1=b.y-a.y,B1=-(b.x-a.x),C1=b.y*a.x-b.x*a.y;
        double A2=d.y-c.y,B2=-(d.x-c.x),C2=d.y*c.x-d.x*c.y;
        double k=A1*B2-A2*B1;
        if(fabs(k)<eps)
        {
            if( fabs( C2*A1-C1*A2)<eps && fabs(B1*C2-C1*B2)<eps ) printf("LINE\n");
            else                   printf("NONE\n");
        }
        else
        {
                double x=-(B1*C2-C1*B2)*1.000000000/k;
                double y=(A1*C2-C1*A2)*1.00000000/k;
                printf("POINT %.2f %.2f\n",x,y);
        }
    }
    
    int main()
    {
        int T;scanf("%d",&T);
        printf("INTERSECTING LINES OUTPUT\n");
        while(T--)
        {
            Point a,b,c,d;
            scanf("%lf %lf %lf %lf",&a.x,&a.y,&b.x,&b.y);
            scanf("%lf %lf %lf %lf",&c.x,&c.y,&d.x,&d.y);
            XX(Line(a,b),Line(c,d));
        }
        printf("END OF OUTPUT\n");
        return 0;
    }

     

    转载于:https://www.cnblogs.com/Andromeda-Galaxy/p/10451643.html

    展开全文
  • 计算几何板子

    2021-09-12 17:46:43
    板子

    板子题链接

    1_A:点积求投影

    1_B:点积求投影后得到所需向量

    1_C:叉积判一下

    2_A:点积和叉积

    2_B:判点在线段上和线段相交

    2_C:判一下线段是否相交,然后用直线相交方法求交点

    2_D:端点之间的距离,点到直线距离,都判一下。

    3_A:叉积算一下

    3_B:叉积判一下

    3_C:点向平行于x轴正方向搞一个射线看一下和多边形的交点奇偶性,对于平行于x轴的边,需要一些细节判一下,详见jls板子

    4_A:维护上下凸壳

    4_B:利用凸多边形的凸函数性质

    5_A:分治

    6_A:离散化+扫描线+pbds

    7_A:根据圆心之间的距离及两个圆的半径判断

    7_B:角平分线求交点,点到直线距离

    7_C:垂直平分线求交点,点到直线距离

    7_D:投影后勾股定理

    7_E:余弦定理

    7_F:相似三角形

    7_G:相似三角形

    7_H:转成三角形和圆求有向面积交

    7_I:扇形减去四边形面积

    #define vec point
    typedef long double db ; //卡精度就换long double
    const db eps = 1e-6 ; //调参
    const db pi = acos(-1.0) ;
    int sgn(db x)
    {
    	if(x < -eps) return -1 ;
    	else if (x > eps)  return 1 ;
    	else return 0 ;
    }
    int cmp(db x , db y)
    {
    	return sgn(x - y) ;
    }
    void print(int num , db x)
    {
    	cout << fixed << setprecision(num) << x << '\n' ;
    }
    struct point
    {
    	db x , y ;
    	point(){}
    	point(db x2 , db y2)
    	{
    		x = x2 , y = y2 ;
    	}
    	point operator + (const point& s)const{return (point){x + s.x , y + s.y} ;}
    	point operator - (const point& s)const{return (point){x - s.x , y - s.y} ;}
    	point operator * (const db& k)const{return (point){x * k , y * k} ;}
    	point operator / (const db& k)const{return (point){x / k , y / k} ;}
    	bool operator < (point b) const
    	{
    		return sgn(x - b.x) == 0 ? sgn(y - b.y) < 0 : x < b.x ;
    	}
    	bool equal(point p2)
    	{
    		return cmp(x , p2.x) == 0 && cmp(y , p2.y) == 0 ;
    	}
    	db get_angle(){return atan2(y , x) ;} //极角
    	db sq(db x) {return x * x ;}
    	db dis(point p) //很吃精度,可能需要1e-10或更小的eps
    	{
    		return sqrtl(sq(x - p.x) + sq(y - p.y)) ; //check sqrt or sqrtl
    	}
    	db len()
    	{
    		return sqrtl(sq(x) + sq(y)) ;
    	}
    	db len2()
    	{
    		return sq(x) + sq(y) ;
    	}
        point unit()
        {
            db w = len() ;
            return (point){x / w , y / w} ;
        }
    	vec rotate_left()  //向量逆时针旋转90度
    	{
    		return vec(-y , x) ;
    	}
    	vec rotate_right()  //向量顺时针旋转90度
    	{
    		return vec(y , -x) ;
    	}
    	point move(db r) //原点沿该点代表的向量方向移动r
    	{
    		db l = len() ;
    		if(sgn(l) == 0)  return *this ;
    		else  return point(x * r / l , y * r / l) ;
    	}
    } ;
    db cross(vec s , vec t){return s.x * t.y - s.y * t.x ;} // 叉积
    db dot(vec s , vec t){return s.x * t.x + s.y * t.y ;} // 点积
    db rad(vec p1 , vec p2){ return atan2(cross(p1 , p2) , dot(p1 , p2)) ;} // 有向弧度
    int in_mid(db k1 , db k2 , db k3) // k3 在 [k1,k2] 内
    {
    	return cmp(k1 , k3) * cmp(k2 , k3) <= 0 ;
    }
    int in_mid(point p1 , point p2 , point p3)
    {
    	return in_mid(p1.x , p2.x , p3.x) && in_mid(p1.y , p2.y , p3.y) ;
    }
    int counter_clockwise(point p1 , point p2 , point p3) // p1 p2 p3 逆时针 1 顺时针 -1 否则 0  
    {
        return sgn(cross(p2 - p1 , p3 - p1)) ;
    }
    struct line
    {
    	point s , e ; // 点s在直线上,方向向量是e ,半平面在e左侧
    	db ang ; //在运算时注意设置ang的值!
    	line() {}
    	line(point s2 , point e2)
    	{
    		s = s2 , e = e2 ;
    		ang = e.get_angle() ;
    	}
    	bool operator < (const line& l1)const{return ang < l1.ang ;}
    	int onRight(point p){return sgn(cross(e , p - s)) < 0 ;} //点p在直线右侧
    	point projection(point p) //p在该直线上的投影点
    	{
    		db t = dot(e , p - s) / e.len() ;
    		return s + e.move(t) ;
    	}
    	point reflection(point p) //p关于直线的对称点
    	{
    		point p2 = projection(p) ;
    		vec v = p2 - p ;
    		return p + v * 2 ;
    	}
    	db dis_point(point p) //点到直线的距离
    	{
    		db t = dot(e , p - s) / e.len() ;
    		point pp = s + e.move(t) ;
    		return p.dis(pp) ;
    	}
    	db dis_segment_point(point p) //点到线段的距离
    	{
    		point p1 = s ; //线段p1->p2
    		point p2 = s + e ;
    		point p3 = projection(p) ;
    		if(in_mid(p1 , p2 , p3))  return p.dis(p3) ;
    		else  return min(p.dis(p1) , p.dis(p2)) ;
    	}
    	bool on_segment(point p) //p在线段s->s+e上
    	{
    		return sgn((p.dis(s) + p.dis(s + e)) - s.dis(s + e)) == 0 ;//
    	}
    	int direction(point p2) //p0 = s , p1 = s + e 判断p2与向量p0->p1的位置关系
    	{
    		point p0 = s ;
    		point p1 = s + e ;
    		db res = cross(p1 - p0 , p2 - p0) ;
    		if(sgn(res) > 0)  return 1 ; //p0->p1逆时针旋转得到p0->p2的方向
    		else if(sgn(res) < 0)  return 2 ; //p0->p1顺时针旋转得到p0->p2的方向
    		else
    		{
    			db res2 = dot(p1 - p0 , p2 - p0) ;
    			if(on_segment(p2))  return 5 ; //p2在线段p0->p1上
    			else if(sgn(res2) < 0)  return 3 ; //p0->p1与p0->p2方向相反
    			else  return 4 ; //p0->p1与p0->p2方向相同
    		}
    	}
    	int Parallel_Orthogonal(line l2) //判断两条直线平行、正交
    	{
    		if(sgn(cross(e , l2.e)) == 0)  return 2 ; //平行
    		else if(sgn(dot(e , l2.e)) == 0)  return 1 ; //正交
    		else  return 0 ;
    	}
    	point line_intersect(point s1 , point t1 , point s2 , point t2) // 直线交点(点+向量表示直线)
    	{
    		db x = cross(t2 , s1 - s2) / cross(t1 , t2) ;
    		return s1 + t1 * x ;
    	} 
    	point line_intersect(line l2) //直线交点(点+向量表示直线)
    	{ 
    		assert(sgn(cross(e , l2.e)) != 0) ;
    		return line_intersect(s , e , l2.s , l2.e) ;
    	}
    	point segment_intersect(line l2) //线段s->s+e与线段l2.s->l2.s+l2.e的交点
    	{
    		assert(sgn(cross(e , l2.e)) != 0) ;
    		return line_intersect(s , e , l2.s , l2.e) ;
    	}
    	bool can_segment_intersect(line l2) //线段s->s+e与线段l2.s->l2.s+l2.e
    	{
    		int t1 = sgn(cross(e , l2.s - s)) ;
    		int t2 = sgn(cross(e , (l2.s + l2.e) - s)) ;
    		int t3 = sgn(cross(l2.e , s - l2.s)) ;
    		int t4 = sgn(cross(l2.e , (s + e) - l2.s)) ;
    		if(on_segment(l2.s) || on_segment(l2.s + l2.e) 
    		|| l2.on_segment(s) || l2.on_segment(s + e))  //端点在另一条线段上的情况
    		return true ;
    		else if(t1 * t2 < 0 && t3 * t4 < 0)  return true ;
    		else  return false ;
    	}
    	db segment_segment_dis(line l2) //线段s->s+e与线段l2.s->l2.s+l2.e
    	{
    		if(can_segment_intersect(l2))  return 0.0 ;
    		else
    		{
    			db res = min({s.dis(l2.s) , s.dis(l2.s + l2.e) , (s + e).dis(l2.s) , (s + e).dis(l2.s + l2.e)}) ;
    			point t = l2.projection(s) ;
    			if(l2.on_segment(t))  res = min(res , s.dis(t)) ;
    			point t2 = l2.projection(s + e) ;
    			if(l2.on_segment(t2))  res = min(res , (s + e).dis(t2)) ;
    			point t3 = projection(l2.s) ;
    			if(on_segment(t3))  res = min(res , l2.s.dis(t3)) ;
    			point t4 = projection(l2.s + l2.e) ;
    			if(on_segment(t4))  res = min(res , (l2.s + l2.e).dis(t4)) ;
    			return res ;
    		}
    	}
    	line trans(db r) //半平面u向内平移r
    	{
    		line v ;
    		v.s = s + e.rotate_left().move(r) ;
    		v.e = e ;
    		v.ang = e.get_angle() ; //注意设置ang的值
    		return v ;
    	}
    } ;
    struct Polygon
    {
    	vector<point> pp ;
    	Polygon() {}
    	Polygon(int n)
    	{
    		pp.resize(n) ;
    	}
    	int size() const
    	{
    		return pp.size() ;
    	}
    	void resize(int n)
    	{
    		pp.resize(n) ;
    	}
    	point operator [](int id) const
    	{
    		if(id < 0 || id >= size())  assert(false) ;
    		else  return pp[id] ;
    	}
    	point &operator [](int id) 
    	{
    		return pp[id] ;
    	}
    	// 凸包
    	// Andrew算法
    	// 按照x优先的顺序排序(坐标从小到大)
    	// 从第一个点开始遍历,如果下一个点在栈顶的两个元素所连成直线的左边,那么就入栈
    	// 否则如果在右边,说明凸包有更优的方案,上次的点出栈,并直到新点在栈顶两个点所在的直线的左边为止
    	// 这样求得下凸包,类似方法求得上凸包。
    	void Convex_Hull(vector<point> P , vector<point> &res , int flag = 1) // flag = 0 不严格(存在三点共线) flag = 1 严格(不存在三点共线) 
    	{
    		int n = P.size() ;
    		res.resize(0) ;
    		res.resize(n * 2) ;
    		sort(P.begin() , P.end()) ; 
    		int now = -1 ;
    		for(int i = 0 ; i < (int)P.size() ; i ++)
    		{
    			while(now > 0 && sgn(cross(res[now] - res[now - 1] , P[i] - res[now - 1])) < flag)  now -- ;
    			res[++ now] = P[i] ;
    		}
    		int pre = now ;
    		for(int i = n - 2 ; i >= 0 ; i --)
    		{
    			while(now > pre && sgn(cross(res[now] - res[now - 1] , P[i] - res[now - 1])) < flag)  now -- ;
    			res[++ now] = P[i] ;
    		}
    		res.resize(now) ;
    	}
    	// 半平面交
    	// 排序增量算法
    	// 1.以逆时针为正方向,建边。(输入方向不确定时,可用叉乘求面积看正负得知输入的顺逆方向。)
    	// 2.对线段根据极角排序。
    	// 3.去除极角相同的情况下,位置在右边的边。
    	// 4.用双端队列储存线段集合L,遍历所有线段。
    	// 5.判断该线段加入后对半平面交的影响,(对双端队列的头部和尾部进行判断,因为线段加入是有序的。)。
    	// 6.如果某条线段对于新的半平面交没有影响,则从队列中剔除掉。
    	// 7.最后剩下的线段集合L,即使最后要求的半平面交。
    	int Half_Plane_Intersection(vector<line>L , vector<point> &res) 	//L.push_back({p[i] , p[(i + 1) % m] - p[i]}) ;
    	{  
    		int n = L.size() ;
    		sort(L.begin() , L.end()) ;
    		int l = 0 , r = 0 ;
    		vector<point> p(n) ;
    		vector<line> q(n) ;
    		q[0] = L[0] ;
    		for(int i = 1 ; i <= n - 1 ; i ++)
    		{
    			while(l < r && L[i].onRight(p[r - 1]))  r -- ;
    			while(l < r && L[i].onRight(p[l]))  l ++ ;
    			q[++ r] = L[i] ;
    			if (sgn(cross(q[r].e , q[r - 1].e)) == 0)
    			{
    				r -- ;
    				if(!q[r].onRight(L[i].s))  q[r] = L[i] ;
    			}
    			if(l < r)  p[r - 1] = q[r - 1].line_intersect(q[r]) ;
    		}
    		while(l < r && q[l].onRight(p[r - 1]))  r -- ;
    		if(r - l <= 1)  return 0 ; // 交不存在
    		res.clear() ;
    		p[r] = q[l].line_intersect(q[r]) ;
    		for(int i = l ; i <= r ; i ++)  res.push_back(p[i]) ;
    		return 1 ;
    	}
    	db perimeter(vector<point> A) // 多边形用 vector<point> 表示 , 逆时针 
    	{
    		db ans = 0 ;
    		for(int i = 0 ; i < (int)A.size() ; i ++)  ans += A[i].dis(A[(i + 1) % (int)A.size()]) ;
    		return ans ;
    	} 
    	db area() // 多边形用 vector<point> 表示 , 逆时针 
    	{
    		db ans = 0 ;
    		for(int i = 0 ; i < (int)pp.size() ; i ++)  ans += cross(pp[i] , pp[(i + 1) % (int)pp.size()]) ;
    		return fabs(ans) / 2 ;
    	}
    	bool is_convex(vector<point> A) // 多边形用 vector<point> 表示 , 逆时针 
    	{
    		int n = A.size() ;
    		assert(n >= 3) ;
    		bool flag = true ;
    		for(int i = 0 ; i < n ; i ++)
    		{
    			int lst = (i - 1 + n) % n ;
    			int nxt = (i + 1) % n ;
    			if(sgn(cross(A[lst] - A[i] , A[nxt] - A[i])) > 0)  flag = false ;
    		}
    		return flag ;
    	}
    	int contain_point(point p0)  // 不一定是凸的 2 内部 1 边界 0 外部
    	{
    		int n = pp.size() ;
    		int res = 0 ;
    		for(int i = 0 ; i < n ; i ++)
    		{
    			point u = pp[(i - 1 + n) % n] ;
    			point v = pp[i] ;
    			line l ;
    			l.s = u ;
    			l.e = v - u ;
    			if(l.on_segment(p0))  return 1 ;
    			if(cmp(u.y , v.y) > 0)  swap(u , v) ;
    			if(cmp(u.y , p0.y) >= 0 || cmp(v.y , p0.y) < 0)  continue ;
    			if(sgn(cross(u - v , p0 - v)) < 0)  res ^= 1 ;
    		}
    		return res << 1 ;
    	}
    	db Convex_Diameter()
    	{
    		//原理是凸多边形的凸函数性质
    		int now = 0 ;
    		int n = pp.size() ;
    		db res = 0 ;
    		for(int i = 0 ; i < n ; i ++)
    		{
    			now = max(now , i) ;
    			while(true)
    			{
    				db k1 = pp[i].dis(pp[now % n]) ;
    				db k2 = pp[i].dis(pp[(now + 1) % n]) ;
    				res = max({res , k1 , k2}) ;
    				if(cmp(k2 , k1) > 0)  now ++ ;
    				else  break ;
    			}
    		}
    		return res ;
    	}
        Polygon line_cut_convex(point p1 , point p2)
        {
            //在已有凸多边形基础上新加一个半平面p1->p2
            //问之后的凸多边形
            int n = pp.size() ;
            line l ;
            l.s = p1 ;
            l.e = p2 - p1 ;
            Polygon res ;
            for(int i = 0 ; i < n ; i ++)
            {
                int w1 = counter_clockwise(p1 , p2 , pp[i]) ;
                int w2 = counter_clockwise(p1 , p2 , pp[(i + 1) % n]) ;
                if(w1 >= 0)  res.pp.push_back(pp[i]) ;
                line l2 ;
                l2.s = pp[i] ;
                l2.e = pp[(i + 1) % n] - pp[i] ;
                if(w1 * w2 < 0)  res.pp.push_back(l.line_intersect(l2)) ;
            }
            return res ;
        }
    } ;
    struct circle
    {
        point o ;
        db r ;
    	db area()
    	{
    		return pi * r * r ;
    	}
    	int inside(point p2)
    	{
    		return cmp(r , o.dis(p2)) ;
    	}
        circle get_circle(point p1 , point p2 , point p3) //三点确定一个圆
        {
            db a1 = p2.x - p1.x , b1 = p2.y - p1.y , c1 = (a1 * a1 + b1 * b1) / 2 ;
            db a2 = p3.x - p1.x , b2 = p3.y - p1.y , c2 = (a2 * a2 + b2 * b2) / 2 ;
            db d = a1 * b2 - a2 * b1 ;
            point o = (point){p1.x + (c1 * b2 - c2 * b1) / d , p1.y + (a1 * c2 - a2 * c1) / d} ;
            return (circle){o , p1.dis(o)} ;
        }
        vector<point> Tangent_Point_Circle(point p) //沿圆的逆时针方向给出切点,点在圆上也给出两个
        {
            //相似三角形
            assert(cmp(p.dis(o) , r) >= 0) ;
            db a = (p - o).len() ;
            assert(sgn(a) >= 0) ;
            db b = r * r / a ;
            db c = sqrt(max(db(0.0) , r * r - b * b)) ;
            point k = (p - o).unit() ;
            point m = o + k * b ;
            point v = k.rotate_left() * c ;
            return {m - v , m + v} ;
        }
    	int check_pos_circle_circle(circle c1 , circle c2) //返回两个圆的公切线数量,即位置关系
    	{
    		//根据圆心之间的距离及两个圆的半径判断
    		if(cmp(c1.r , c2.r) < 0)  swap(c1 , c2) ;
    		db d = c1.o.dis(c2.o) ;
    		int w1 = cmp(d , c1.r + c2.r) ;
    		int w2 = cmp(d , c1.r - c2.r) ;
    		if(w1 > 0)  return 4 ; //相离
    		else if(w1 == 0)  return 3 ; //外切
    		else if(w2 > 0)  return 2 ; //相交
    		else if(w2 == 0)  return 1 ; //内切
    		else  return 0 ; //内含
    	}
    	vector<point> get_circle_line(line l) //圆与直线求交点,沿l.e方向的顺序给出交点,相切给出两个
    	{
    		//投影后勾股定理
    		point k = l.projection(o) ;
    		db d = r * r - (k - o).len2() ;
    		if(sgn(d) < 0)  return {} ;
    		vec v = l.e.unit() * sqrtl(max((db)0.0 , d)) ;
    		return {k - v , k + v} ;
    	}
    	vector<point> get_circle_circle(circle c1 , circle c2) //两个圆的交点,沿圆c1逆时针给出,相切给出两个
    	{
    		//余弦定理
    		int pos = check_pos_circle_circle(c1 , c2) ;
    		if(pos == 0 || pos == 4)  return {} ;
    		db a = (c2.o - c1.o).len2() ;
    		db cosA = (c1.r * c1.r + a - c2.r * c2.r) / (2 * c1.r * sqrtl(max(a , (db)0.0))) ;
    		db b = c1.r * cosA ;
    		db c = sqrtl(max((db)0.0 , c1.r * c1.r - b * b)) ;
    		point k = (c2.o - c1.o).unit() ;
    		point m = c1.o + k * b ;
    		vec v = k.rotate_left() * c ;
    		return {m - v , m + v} ;
    	}
    	circle Incircle_of_a_Triangle(point p1 , point p2 , point p3) //三角形p1p2p3内接圆
    	{
    		//角平分线求交点,点到直线距离
    		line l1 ;
    		l1.s = p1 ;
    		l1.e = (p2 - p1).unit() + (p3 - p1).unit() ;
    		line l2 ;
    		l2.s = p2 ;
    		l2.e = (p1 - p2).unit() + (p3 - p2).unit() ;
    		point p = l1.line_intersect(l2) ;
    		line l3 ;
    		l3.s = p1 ;
    		l3.e = p2 - p1 ;
    		circle cc ;
    		cc.o = p ;
    		cc.r = l3.dis_point(p) ;
    		return cc ;
    	}
    	circle Circumscribed_Circle_of_a_Triangle(point p1 , point p2 , point p3) //三角形p1p2p3外接圆
    	{
    		//垂直平分线求交点,点到直线距离
    		line l1 ;
    		l1.s = p1 + (p2 - p1) / 2 ;
    		l1.e = (p2 - p1).rotate_left() ; //顺时针或逆时针旋转90度均可,因为是直线求交
    		line l2 ;
    		l2.s = p2 + (p3 - p2) / 2 ;
    		l2.e = (p3 - p2).rotate_left() ;
    		point p = l1.line_intersect(l2) ;
    		circle cc ;
    		cc.o = p ;
    		cc.r = p.dis(p1) ;
    		return cc ;
    	}
    	vector<line> Tangent_out_Circle_Circle(circle c1 , circle c2) //外切线
    	{
    		//相似三角形
    		int pos = check_pos_circle_circle(c1 , c2) ;
    		if(pos == 0)  return {} ;
    		if(pos == 1)
    		{
    			point k = get_circle_circle(c1 , c2)[0] ;
    			vec v = k - c1.o ;
    			v = v.rotate_left() ;
    			line l ;
    			l.s = k ;
    			l.e = v ;
    			return {l} ;
    		}
    		if(cmp(c1.r , c2.r) == 0)
    		{
    			vec v = (c2.o - c1.o).rotate_left().unit() * c1.r ;
    			vec p1 = c1.o + v ;
    			line l1 ;
    			l1.s = p1 ;
    			l1.e = c2.o - c1.o ;
    			vec p2 = c1.o - v ;
    			line l2 ;
    			l2.s = p2 ;
    			l2.e = c2.o - c1.o ;
    			return {l1 , l2} ;
    		}
    		else
    		{
    			vec v = c2.o - c1.o ;
    			v = v.unit() * (c1.o.dis(c2.o) * c1.r / (c1.r - c2.r)) ;
    			point H = c1.o + v ;
    			vector<point> pp1 = c1.Tangent_Point_Circle(H) ;
    			line l1 ;
    			l1.s = H ;
    			l1.e = pp1[0] - H ;
    			line l2 ;
    			l2.s = H ;
    			l2.e = pp1[1] - H ;
    			return {l1 , l2} ;
    		}
    	}
    	vector<line> Tangent_in_Circle_Circle(circle c1 , circle c2) //内切线
    	{
    		int pos = check_pos_circle_circle(c1 , c2) ;
    		if(pos <= 2)  return {} ;
    		if(pos == 3)
    		{
    			point k = get_circle_circle(c1 , c2)[0] ;
    			vec v = k - c1.o ;
    			v = v.rotate_left() ;
    			line l ;
    			l.s = k ;
    			l.e = v ;
    			return {l} ;
    		}
    		db t = c1.r / (c1.r + c2.r) * c1.o.dis(c2.o) ;
    		point p = c1.o + (c2.o - c1.o).unit() * t ;
    		vector<point> pp1 = c1.Tangent_Point_Circle(p) ;
    		line l1 ;
    		l1.s = p ;
    		l1.e = pp1[0] - p ;
    		line l2 ;
    		l2.s = p ;
    		l2.e = pp1[1] - p ;
    		return {l1 , l2} ;
    	}
    	vector<line> Tangent_Circle_Circle(circle c1 , circle c2) //两个圆的公切线
    	{
    		//外切线和内切线一起考虑
    		if(cmp(c1.r , c2.r) < 0)  swap(c1 , c2) ;
    		vector<line> A = Tangent_out_Circle_Circle(c1 , c2) ;
    		vector<line> B = Tangent_in_Circle_Circle(c1 , c2) ;
    		for(auto u : B)  A.push_back(u) ;
    		return A ;
    	}
    	db get_area_of_Circle_Triangle(circle c1 , point p2 , point p3)
    	{
    		//圆c1与三角形c1 p2 p3的交
    		point k = c1.o ;
    		c1.o = c1.o - k ;
    		p2 = p2 - k ;
    		p3 = p3 - k ;
    		int pos1 = c1.inside(p2) ;
    		int pos2 = c1.inside(p3) ;
    		line l ;
    		l.s = p2 ;
    		l.e = p3 - p2 ;
    		vector<point> pp = c1.get_circle_line(l) ;
    		if(pos1 >= 0)
    		{
    			if(pos2 >= 0)  return cross(p2 , p3) / 2 ;
    			return c1.r * c1.r * rad(pp[1] , p3) / 2 + cross(p2 , pp[1]) / 2 ;
    		}
    		else if(pos2 >= 0)
    		{
    			return c1.r * c1.r * rad(p2 , pp[0]) / 2 + cross(pp[0] , p3) / 2 ;
    		}
    		else
    		{
    			line l ;
    			l.s = p2 ;
    			l.e = p3 - p2 ;
    			int pos = cmp(c1.r , l.dis_segment_point(c1.o)) ;
    			if(pos <= 0)  return c1.r * c1.r * rad(p2 , p3) / 2 ;
    			return cross(pp[0] , pp[1]) / 2 + c1.r * c1.r * (rad(p2 , pp[0]) + rad(pp[1] , p3)) / 2 ;
    		}
    	}
    	db get_area_of_Circle_Polygon(circle c1 , vector<point> pp) //pp按照逆时针给出
    	{
    		//转成三角形和圆求有向面积交
    		int n = pp.size() ;
    		db res = 0 ;
    		for(int i = 0 ; i < n ; i ++)
    		{
    			point now = pp[i] ; //三角形o pp[i] pp[(i + 1) % n] 需要按照顺序给出
    			point nxt = pp[(i + 1) % n] ;
    			res += get_area_of_Circle_Triangle(c1 , now , nxt) ; //已经除2
    		}
    		return fabs(res) ;
    	}
    	db get_area_of_Circle_Circle(circle c1 , circle c2) //两个圆的相交面积
    	{
    		//扇形减去四边形面积
    		if(cmp(c1.r , c2.r) < 0)  swap(c1 , c2) ;
    		int pos = check_pos_circle_circle(c1 , c2) ;
    		if(pos >= 3)  return 0.0 ;
    		if(pos <= 1)  return c2.area() ;
    		vector<point> pp = get_circle_circle(c1 , c2) ;
    		db S1 = fabs(2.0 * rad(pp[0] - c1.o , c2.o - c1.o)) * c1.r * c1.r / 2 ;
    		db S2 = fabs(2.0 * rad(pp[0] - c2.o , c1.o - c2.o)) * c2.r * c2.r / 2 ;
    		db S3 = fabs(cross(pp[0] - c1.o , c2.o - c1.o)) ;
    		return S1 + S2 - S3 ;
    	}
    } c ;
    db closest_point(vector<point> &A , int l , int r)
    {
    	if(r - l <= 5)
    	{
    		db res = 1e20 ;
    		for(int i = l ; i <= r ; i ++)
    			for(int j = i + 1 ; j <= r ; j ++)
    				res = min(res , A[i].dis(A[j])) ;	
    		return res ;
    	}
    	int mid = (l + r) / 2 ;
    	db res = min(closest_point(A , l , mid) , closest_point(A , mid + 1 , r)) ;
    	vector<point> B ;
    	for(int i = l ; i <= r ; i ++)
    		if(cmp(fabs(A[i].x - A[mid].x) , res) <= 0)
    			B.push_back(A[i]) ;
    	sort(B.begin() , B.end() , [&](point p1 , point p2){ return cmp(p1.y , p2.y) < 0 ; }) ;
    	for(int i = 0 ; i < B.size() ; i ++)
    		for(int j = i + 1 ; j < B.size() && B[j].y - B[i].y < res ; j ++)
    			res = min(res , B[i].dis(B[j])) ;
    	return res ;
    }

    展开全文
  • 包含点、线、多边形、凸多边形、圆、半平面交等算法的板子
  • 计算几何模板

    2018-09-16 04:30:00
    计算几何模板 if( (a.y2-a.y1)*(b.x2-b.x1)==(b.y2-b.y1)*(a.x2-a.x1) ) //向量判别是否平行 x1y2=x2y1 { if( (a.y1-b.y2)*(a.x2-b.x1)==(a.y2-b.y1)*(a.x1-b.x2) && (b.y1-a.y1)*(a.x2-b.x2)==(a.y2-b.y2)*(b.x1-a...
  • 转自:... 二维几何 // `计算几何模板` const double eps = 1e-8; const double inf = 1e20; const double pi = acos(-1.0); const int maxp = 1010; //`Compares a double ...

    转自:https://kuangbin.github.io/2019/04/28/20190428/#more

    二维几何

    // `计算几何模板`
    const double eps = 1e-8;
    const double pi = acos(-1.0);
    const int maxp = 1010;
    //`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;}
    /*
     * Point
     * Point()               - Empty constructor
     * Point(double _x,double _y)  - constructor
     * input()             - double input
     * output()            - %.2f output
     * operator ==         - compares x and y
     * operator <          - compares first by x, then by y
     * operator -          - return new Point after subtracting curresponging x and y
     * operator ^          - cross product of 2d points
     * operator *          - dot product
     * len()               - gives length from origin
     * len2()              - gives square of length from origin
     * distance(Point p)   - gives distance from p
     * operator + Point b  - returns new Point after adding curresponging x and y
     * operator * double k - returns new Point after multiplieing x and y by k
     * operator / double k - returns new Point after divideing x and y by k
     * rad(Point a,Point b)- returns the angle of Point a and Point b from this Point
     * trunc(double r)     - return Point that if truncated the distance from center to r
     * rotleft()           - returns 90 degree ccw rotated point
     * rotright()          - returns 90 degree cw rotated point
     * rotate(Point p,double angle) - returns Point after rotateing the Point centering at p by angle radian ccw
     */
    struct Point{
    	double x,y;
    	Point(){}
    	Point(double _x,double _y){
    		x = _x;
    		y = _y;
    	}
    	void input(){
    		scanf("%lf%lf",&x,&y);
    	}
    	void output(){
    		printf("%.2f %.2f\n",x,y);
    	}
    	bool operator == (Point b)const{
    		return sgn(x-b.x) == 0 && sgn(y-b.y) == 0;
    	}
    	bool operator < (Point b)const{
    		return sgn(x-b.x)== 0?sgn(y-b.y)<0:x<b.x;
    	}
    	Point operator -(const Point &b)const{
    		return Point(x-b.x,y-b.y);
    	}
    	//叉积
    	double operator ^(const Point &b)const{
    		return x*b.y - y*b.x;
    	}
    	//点积
    	double operator *(const Point &b)const{
    		return x*b.x + y*b.y;
    	}
    	//返回长度
    	double len(){
    		return hypot(x,y);//库函数
    	}
    	//返回长度的平方
    	double len2(){
    		return x*x + y*y;
    	}
    	//返回两点的距离
    	double distance(Point p){
    		return hypot(x-p.x,y-p.y);
    	}
    	Point operator +(const Point &b)const{
    		return Point(x+b.x,y+b.y);
    	}
    	Point operator *(const double &k)const{
    		return Point(x*k,y*k);
    	}
    	Point operator /(const double &k)const{
    		return Point(x/k,y/k);
    	}
    	//`计算pa  和  pb 的夹角`
    	//`就是求这个点看a,b 所成的夹角`
    	//`测试 LightOJ1203`
    	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(x*r,y*r);
    	}
    	//`逆时针旋转90度`
    	Point rotleft(){
    		return Point(-y,x);
    	}
    	//`顺时针旋转90度`
    	Point rotright(){
    		return Point(y,-x);
    	}
    	//`绕着p点逆时针旋转angle`
    	Point rotate(Point p,double angle){
    		Point v = (*this) - p;
    		double c = cos(angle), s = sin(angle);
    		return Point(p.x + v.x*c - v.y*s,p.y + v.x*s + v.y*c);
    	}
    };
    /*
     * Stores two points
     * Line()                         - Empty constructor
     * Line(Point _s,Point _e)        - Line through _s and _e
     * operator ==                    - checks if two points are same
     * Line(Point p,double angle)     - one end p , another end at angle degree
     * Line(double a,double b,double c) - Line of equation ax + by + c = 0
     * input()                        - inputs s and e
     * adjust()                       - orders in such a way that s < e
     * length()                       - distance of se
     * angle()                        - return 0 <= angle < pi
     * relation(Point p)              - 3 if point is on line
     *                                  1 if point on the left of line
     *                                  2 if point on the right of line
     * pointonseg(double p)           - return true if point on segment
     * parallel(Line v)               - return true if they are parallel
     * segcrossseg(Line v)            - returns 0 if does not intersect
     *                                  returns 1 if non-standard intersection
     *                                  returns 2 if intersects
     * linecrossseg(Line v)           - line and seg
     * linecrossline(Line v)          - 0 if parallel
     *                                  1 if coincides
     *                                  2 if intersects
     * crosspoint(Line v)             - returns intersection point
     * dispointtoline(Point p)        - distance from point p to the line
     * dispointtoseg(Point p)         - distance from p to the segment
     * dissegtoseg(Line v)            - distance of two segment
     * lineprog(Point p)              - returns projected point p on se line
     * symmetrypoint(Point p)         - returns reflection point of p over se
     *
     */
    struct Line{
    	Point s,e;
    	Line(){}
    	Line(Point _s,Point _e){
    		s = _s;
    		e = _e;
    	}
    	bool operator ==(Line v){
    		return (s == v.s)&&(e == v.e);
    	}
    	//`根据一个点和倾斜角angle确定直线,0<=angle<pi`
    	Line(Point p,double angle){
    		s = p;
    		if(sgn(angle-pi/2) == 0){
    			e = (s + Point(0,1));
    		}
    		else{
    			e = (s + Point(1,tan(angle)));
    		}
    	}
    	//ax+by+c=0
    	Line(double a,double b,double c){
    		if(sgn(a) == 0){
    			s = Point(0,-c/b);
    			e = Point(1,-c/b);
    		}
    		else if(sgn(b) == 0){
    			s = Point(-c/a,0);
    			e = Point(-c/a,1);
    		}
    		else{
    			s = Point(0,-c/b);
    			e = Point(1,(-c-a)/b);
    		}
    	}
    	void input(){
    		s.input();
    		e.input();
    	}
    	void adjust(){
    		if(e < s)swap(s,e);
    	}
    	//求线段长度
    	double length(){
    		return s.distance(e);
    	}
    	//`返回直线倾斜角 0<=angle<pi`
    	double angle(){
    		double k = atan2(e.y-s.y,e.x-s.x);
    		if(sgn(k) < 0)k += pi;
    		if(sgn(k-pi) == 0)k -= pi;
    		return k;
    	}
    	//`点和直线关系`
    	//`1  在左侧`
    	//`2  在右侧`
    	//`3  在直线上`
    	int relation(Point p){
    		int c = sgn((p-s)^(e-s));
    		if(c < 0)return 1;
    		else if(c > 0)return 2;
    		else return 3;
    	}
    	// 点在线段上的判断
    	bool pointonseg(Point p){
    		return sgn((p-s)^(e-s)) == 0 && sgn((p-s)*(p-e)) <= 0;
    	}
    	//`两向量平行(对应直线平行或重合)`
    	bool parallel(Line v){
    		return sgn((e-s)^(v.e-v.s)) == 0;
    	}
    	//`两线段相交判断`
    	//`2 规范相交`
    	//`1 非规范相交`
    	//`0 不相交`
    	int segcrossseg(Line v){
    		int d1 = sgn((e-s)^(v.s-s));
    		int d2 = sgn((e-s)^(v.e-s));
    		int d3 = sgn((v.e-v.s)^(s-v.s));
    		int d4 = sgn((v.e-v.s)^(e-v.s));
    		if( (d1^d2)==-2 && (d3^d4)==-2 )return 2;
    		return (d1==0 && sgn((v.s-s)*(v.s-e))<=0) ||
    			(d2==0 && sgn((v.e-s)*(v.e-e))<=0) ||
    			(d3==0 && sgn((s-v.s)*(s-v.e))<=0) ||
    			(d4==0 && sgn((e-v.s)*(e-v.e))<=0);
    	}
    	//`直线和线段相交判断`
    	//`-*this line   -v seg`
    	//`2 规范相交`
    	//`1 非规范相交`
    	//`0 不相交`
    	int linecrossseg(Line v){
    		int d1 = sgn((e-s)^(v.s-s));
    		int d2 = sgn((e-s)^(v.e-s));
    		if((d1^d2)==-2) return 2;
    		return (d1==0||d2==0);
    	}
    	//`两直线关系`
    	//`0 平行`
    	//`1 重合`
    	//`2 相交`
    	int linecrossline(Line v){
    		if((*this).parallel(v))
    			return v.relation(s)==3;
    		return 2;
    	}
    	//`求两直线的交点`
    	//`要保证两直线不平行或重合`
    	Point crosspoint(Line v){
    		double a1 = (v.e-v.s)^(s-v.s);
    		double a2 = (v.e-v.s)^(e-v.s);
    		return Point((s.x*a2-e.x*a1)/(a2-a1),(s.y*a2-e.y*a1)/(a2-a1));
    	}
    	//点到直线的距离
    	double dispointtoline(Point p){
    		return fabs((p-s)^(e-s))/length();
    	}
    	//点到线段的距离
    	double dispointtoseg(Point p){
    		if(sgn((p-s)*(e-s))<0 || sgn((p-e)*(s-e))<0)
    			return min(p.distance(s),p.distance(e));
    		return dispointtoline(p);
    	}
    	//`返回线段到线段的距离`
    	//`前提是两线段不相交,相交距离就是0了`
    	double dissegtoseg(Line v){
    		return min(min(dispointtoseg(v.s),dispointtoseg(v.e)),min(v.dispointtoseg(s),v.dispointtoseg(e)));
    	}
    	//`返回点p在直线上的投影`
    	Point lineprog(Point p){
    		return s + ( ((e-s)*((e-s)*(p-s)))/((e-s).len2()) );
    	}
    	//`返回点p关于直线的对称点`
    	Point symmetrypoint(Point p){
    		Point q = lineprog(p);
    		return Point(2*q.x-p.x,2*q.y-p.y);
    	}
    };
    //圆
    struct circle{
    	Point p;//圆心
    	double r;//半径
    	circle(){}
    	circle(Point _p,double _r){
    		p = _p;
    		r = _r;
    	}
    	circle(double x,double y,double _r){
    		p = Point(x,y);
    		r = _r;
    	}
    	//`三角形的外接圆`
    	//`需要Point的+ /  rotate()  以及Line的crosspoint()`
    	//`利用两条边的中垂线得到圆心`
    	//`测试:UVA12304`
    	circle(Point a,Point b,Point c){
    		Line u = Line((a+b)/2,((a+b)/2)+((b-a).rotleft()));
    		Line v = Line((b+c)/2,((b+c)/2)+((c-b).rotleft()));
    		p = u.crosspoint(v);
    		r = p.distance(a);
    	}
    	//`三角形的内切圆`
    	//`参数bool t没有作用,只是为了和上面外接圆函数区别`
    	//`测试:UVA12304`
    	circle(Point a,Point b,Point c,bool t){
    		Line u,v;
    		double m = atan2(b.y-a.y,b.x-a.x), n = atan2(c.y-a.y,c.x-a.x);
    		u.s = a;
    		u.e = u.s + Point(cos((n+m)/2),sin((n+m)/2));
    		v.s = b;
    		m = atan2(a.y-b.y,a.x-b.x) , n = atan2(c.y-b.y,c.x-b.x);
    		v.e = v.s + Point(cos((n+m)/2),sin((n+m)/2));
    		p = u.crosspoint(v);
    		r = Line(a,b).dispointtoseg(p);
    	}
    	//输入
    	void input(){
    		p.input();
    		scanf("%lf",&r);
    	}
    	//输出
    	void output(){
    		printf("%.2lf %.2lf %.2lf\n",p.x,p.y,r);
    	}
    	bool operator == (circle v){
    		return (p==v.p) && sgn(r-v.r)==0;
    	}
    	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;
    	}
    	//`两圆的关系`
    	//`5 相离`
    	//`4 外切`
    	//`3 相交`
    	//`2 内切`
    	//`1 内含`
    	//`需要Point的distance`
    	//`测试:UVA12304`
    	int relationcircle(circle v){
    		double d = p.distance(v.p);
    		if(sgn(d-r-v.r) > 0)return 5;
    		if(sgn(d-r-v.r) == 0)return 4;
    		double l = fabs(r-v.r);
    		if(sgn(d-r-v.r)<0 && sgn(d-l)>0)return 3;
    		if(sgn(d-l)==0)return 2;
    		if(sgn(d-l)<0)return 1;
    	}
    	//`求两个圆的交点,返回0表示没有交点,返回1是一个交点,2是两个交点`
    	//`需要relationcircle`
    	//`测试:UVA12304`
    	int pointcrosscircle(circle v,Point &p1,Point &p2){
    		int rel = relationcircle(v);
    		if(rel == 1 || rel == 5)return 0;
    		double d = p.distance(v.p);
    		double l = (d*d+r*r-v.r*v.r)/(2*d);
    		double h = sqrt(r*r-l*l);
    		Point tmp = p + (v.p-p).trunc(l);
    		p1 = tmp + ((v.p-p).rotleft().trunc(h));
    		p2 = tmp + ((v.p-p).rotright().trunc(h));
    		if(rel == 2 || rel == 4)
    			return 1;
    		return 2;
    	}
    	//`求直线和圆的交点,返回交点个数`
    	int pointcrossline(Line v,Point &p1,Point &p2){
    		if(!(*this).relationline(v))return 0;
    		Point a = v.lineprog(p);
    		double d = v.dispointtoline(p);
    		d = sqrt(r*r-d*d);
    		if(sgn(d) == 0){
    			p1 = a;
    			p2 = a;
    			return 1;
    		}
    		p1 = a + (v.e-v.s).trunc(d);
    		p2 = a - (v.e-v.s).trunc(d);
    		return 2;
    	}
    	//`得到过a,b两点,半径为r1的两个圆`
    	int gercircle(Point a,Point b,double r1,circle &c1,circle &c2){
    		circle x(a,r1),y(b,r1);
    		int t = x.pointcrosscircle(y,c1.p,c2.p);
    		if(!t)return 0;
    		c1.r = c2.r = r1;
    		return t;
    	}
    	//`得到与直线u相切,过点q,半径为r1的圆`
    	//`测试:UVA12304`
    	int getcircle(Line u,Point q,double r1,circle &c1,circle &c2){
    		double dis = u.dispointtoline(q);
    		if(sgn(dis-r1*2)>0)return 0;
    		if(sgn(dis) == 0){
    			c1.p = q + ((u.e-u.s).rotleft().trunc(r1));
    			c2.p = q + ((u.e-u.s).rotright().trunc(r1));
    			c1.r = c2.r = r1;
    			return 2;
    		}
    		Line u1 = Line((u.s + (u.e-u.s).rotleft().trunc(r1)),(u.e + (u.e-u.s).rotleft().trunc(r1)));
    		Line u2 = Line((u.s + (u.e-u.s).rotright().trunc(r1)),(u.e + (u.e-u.s).rotright().trunc(r1)));
    		circle cc = circle(q,r1);
    		Point p1,p2;
    		if(!cc.pointcrossline(u1,p1,p2))cc.pointcrossline(u2,p1,p2);
    		c1 = circle(p1,r1);
    		if(p1 == p2){
    			c2 = c1;
    			return 1;
    		}
    		c2 = circle(p2,r1);
    		return 2;
    	}
    	//`同时与直线u,v相切,半径为r1的圆`
    	//`测试:UVA12304`
    	int getcircle(Line u,Line v,double r1,circle &c1,circle &c2,circle &c3,circle &c4){
    		if(u.parallel(v))return 0;//两直线平行
    		Line u1 = Line(u.s + (u.e-u.s).rotleft().trunc(r1),u.e + (u.e-u.s).rotleft().trunc(r1));
    		Line u2 = Line(u.s + (u.e-u.s).rotright().trunc(r1),u.e + (u.e-u.s).rotright().trunc(r1));
    		Line v1 = Line(v.s + (v.e-v.s).rotleft().trunc(r1),v.e + (v.e-v.s).rotleft().trunc(r1));
    		Line v2 = Line(v.s + (v.e-v.s).rotright().trunc(r1),v.e + (v.e-v.s).rotright().trunc(r1));
    		c1.r = c2.r = c3.r = c4.r = r1;
    		c1.p = u1.crosspoint(v1);
    		c2.p = u1.crosspoint(v2);
    		c3.p = u2.crosspoint(v1);
    		c4.p = u2.crosspoint(v2);
    		return 4;
    	}
    	//`同时与不相交圆cx,cy相切,半径为r1的圆`
    	//`测试:UVA12304`
    	int getcircle(circle cx,circle cy,double r1,circle &c1,circle &c2){
    		circle x(cx.p,r1+cx.r),y(cy.p,r1+cy.r);
    		int t = x.pointcrosscircle(y,c1.p,c2.p);
    		if(!t)return 0;
    		c1.r = c2.r = r1;
    		return t;
    	}
    
    	//`过一点作圆的切线(先判断点和圆的关系)`
    	//`测试:UVA12304`
    	int tangentline(Point q,Line &u,Line &v){
    		int x = relation(q);
    		if(x == 2)return 0;
    		if(x == 1){
    			u = Line(q,q + (q-p).rotleft());
    			v = u;
    			return 1;
    		}
    		double d = p.distance(q);
    		double l = r*r/d;
    		double h = sqrt(r*r-l*l);
    		u = Line(q,p + ((q-p).trunc(l) + (q-p).rotleft().trunc(h)));
    		v = Line(q,p + ((q-p).trunc(l) + (q-p).rotright().trunc(h)));
    		return 2;
    	}
    	//`求两圆相交的面积`
    	double areacircle(circle v){
    		int rel = relationcircle(v);
    		if(rel >= 4)return 0.0;
    		if(rel <= 2)return min(area(),v.area());
    		double d = p.distance(v.p);
    		double hf = (r+v.r+d)/2.0;
    		double ss = 2*sqrt(hf*(hf-r)*(hf-v.r)*(hf-d));
    		double a1 = acos((r*r+d*d-v.r*v.r)/(2.0*r*d));
    		a1 = a1*r*r;
    		double a2 = acos((v.r*v.r+d*d-r*r)/(2.0*v.r*d));
    		a2 = a2*v.r*v.r;
    		return a1+a2-ss;
    	}
    	//`求两圆相交的面积(精度更高)(需要long double)`
    	double areacircle2(circle v)
    	{
    		double a=hypot(p.x-v.p.x,p.y-v.p.y),b=r,c=v.r;
    		double s1=pi*r*r,s2=pi*v.r*v.r;
    		if(sgn(a-b-c)>=0)
    			return 0;
    		if(sgn(a+min(b,c)-max(b,c))<=0)
    			return min(s1,s2);
    		else
    		{
    			double cta1=2*acos((a*a+b*b-c*c)/(2*a*b));
    			double cta2=2*acos((a*a+c*c-b*b)/(2*a*c));
    			return cta1/(2*pi)*s1-0.5*sin(cta1)*b*b+cta2/(2*pi)*s2-0.5*sin(cta2)*c*c;
    		}
    	}
    	//`求圆和三角形pab的相交面积`
    	//`测试:POJ3675 HDU3982 HDU2892`
    	double areatriangle(Point a,Point b){
    		if(sgn((p-a)^(p-b)) == 0)return 0.0;
    		Point q[5];
    		int len = 0;
    		q[len++] = a;
    		Line l(a,b);
    		Point p1,p2;
    		if(pointcrossline(l,q[1],q[2])==2){
    			if(sgn((a-q[1])*(b-q[1]))<0)q[len++] = q[1];
    			if(sgn((a-q[2])*(b-q[2]))<0)q[len++] = q[2];
    		}
    		q[len++] = b;
    		if(len == 4 && sgn((q[0]-q[1])*(q[2]-q[1]))>0)swap(q[1],q[2]);
    		double res = 0;
    		for(int i = 0;i < len-1;i++){
    			if(relation(q[i])==0||relation(q[i+1])==0){
    				double arg = p.rad(q[i],q[i+1]);
    				res += r*r*arg/2.0;
    			}
    			else{
    				res += fabs((q[i]-p)^(q[i+1]-p))/2.0;
    			}
    		}
    		return res;
    	}
    };
    
    /*
     * n,p  Line l for each side
     * input(int _n)                        - inputs _n size polygon
     * add(Point q)                         - adds a point at end of the list
     * getline()                            - populates line array
     * cmp                                  - comparision in convex_hull order
     * norm()                               - sorting in convex_hull order
     * getconvex(polygon &convex)           - returns convex hull in convex
     * Graham(polygon &convex)              - returns convex hull in convex
     * isconvex()                           - checks if convex
     * relationpoint(Point q)               - returns 3 if q is a vertex
     *                                                2 if on a side
     *                                                1 if inside
     *                                                0 if outside
     * convexcut(Line u,polygon &po)        - left side of u in po
     * gercircumference()                   - returns side length
     * getarea()                            - returns area
     * getdir()                             - returns 0 for cw, 1 for ccw
     * getbarycentre()                      - returns barycenter
     *
     */
    struct polygon{
    	int n;
    	Point p[maxp];
    	Line l[maxp];
    	void input(int _n){
    		n = _n;
    		for(int i = 0;i < n;i++)
    			p[i].input();
    	}
    	void add(Point q){
    		p[n++] = q;
    	}
    	void getline(){
    		for(int i = 0;i < n;i++){
    			l[i] = Line(p[i],p[(i+1)%n]);
    		}
    	}
    	struct cmp{
    		Point p;
    		cmp(const Point &p0){p = p0;}
    		bool operator()(const Point &aa,const Point &bb){
    			Point a = aa, b = bb;
    			int d = sgn((a-p)^(b-p));
    			if(d == 0){
    				return sgn(a.distance(p)-b.distance(p)) < 0;
    			}
    			return d > 0;
    		}
    	};
    	//`进行极角排序`
    	//`首先需要找到最左下角的点`
    	//`需要重载号好Point的 < 操作符(min函数要用) `
    	void norm(){
    		Point mi = p[0];
    		for(int i = 1;i < n;i++)mi = min(mi,p[i]);
    		sort(p,p+n,cmp(mi));
    	}
    	//`得到凸包`
    	//`得到的凸包里面的点编号是0$\sim$n-1的`
    	//`两种凸包的方法`
    	//`注意如果有影响,要特判下所有点共点,或者共线的特殊情况`
    	//`测试 LightOJ1203  LightOJ1239`
    	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();//`原来得到的是顺时针的点,排序后逆时针`
    	}
    	//`得到凸包的另外一种方法`
    	//`测试 LightOJ1203  LightOJ1239`
    	void Graham(polygon &convex){
    		norm();
    		int &top = convex.n;
    		top = 0;
    		if(n == 1){
    			top = 1;
    			convex.p[0] = p[0];
    			return;
    		}
    		if(n == 2){
    			top = 2;
    			convex.p[0] = p[0];
    			convex.p[1] = p[1];
    			if(convex.p[0] == convex.p[1])top--;
    			return;
    		}
    		convex.p[0] = p[0];
    		convex.p[1] = p[1];
    		top = 2;
    		for(int i = 2;i < n;i++){
    			while( top > 1 && sgn((convex.p[top-1]-convex.p[top-2])^(p[i]-convex.p[top-2])) <= 0 )
    				top--;
    			convex.p[top++] = p[i];
    		}
    		if(convex.n == 2 && (convex.p[0] == convex.p[1]))convex.n--;//特判
    	}
    	//`判断是不是凸的`
    	bool isconvex(){
    		bool s[3];
    		memset(s,false,sizeof(s));
    		for(int i = 0;i < n;i++){
    			int j = (i+1)%n;
    			int k = (j+1)%n;
    			s[sgn((p[j]-p[i])^(p[k]-p[i]))+1] = true;
    			if(s[0] && s[2])return false;
    		}
    		return true;
    	}
    	//`判断点和任意多边形的关系`
    	//` 3 点上`
    	//` 2 边上`
    	//` 1 内部`
    	//` 0 外部`
    	int relationpoint(Point q){
    		for(int i = 0;i < n;i++){
    			if(p[i] == q)return 3;
    		}
    		getline();
    		for(int i = 0;i < n;i++){
    			if(l[i].pointonseg(q))return 2;
    		}
    		int cnt = 0;
    		for(int i = 0;i < n;i++){
    			int j = (i+1)%n;
    			int k = sgn((q-p[j])^(p[i]-p[j]));
    			int u = sgn(p[i].y-q.y);
    			int v = sgn(p[j].y-q.y);
    			if(k > 0 && u < 0 && v >= 0)cnt++;
    			if(k < 0 && v < 0 && u >= 0)cnt--;
    		}
    		return cnt != 0;
    	}
    	//`直线u切割凸多边形左侧`
    	//`注意直线方向`
    	//`测试:HDU3982`
    	void convexcut(Line u,polygon &po){
    		int &top = po.n;//注意引用
    		top = 0;
    		for(int i = 0;i < n;i++){
    			int d1 = sgn((u.e-u.s)^(p[i]-u.s));
    			int d2 = sgn((u.e-u.s)^(p[(i+1)%n]-u.s));
    			if(d1 >= 0)po.p[top++] = p[i];
    			if(d1*d2 < 0)po.p[top++] = u.crosspoint(Line(p[i],p[(i+1)%n]));
    		}
    	}
    	//`得到周长`
    	//`测试 LightOJ1239`
    	double getcircumference(){
    		double sum = 0;
    		for(int i = 0;i < n;i++){
    			sum += p[i].distance(p[(i+1)%n]);
    		}
    		return sum;
    	}
    	//`得到面积`
    	double getarea(){
    		double sum = 0;
    		for(int i = 0;i < n;i++){
    			sum += (p[i]^p[(i+1)%n]);
    		}
    		return fabs(sum)/2;
    	}
    	//`得到方向`
    	//` 1 表示逆时针,0表示顺时针`
    	bool getdir(){
    		double sum = 0;
    		for(int i = 0;i < n;i++)
    			sum += (p[i]^p[(i+1)%n]);
    		if(sgn(sum) > 0)return 1;
    		return 0;
    	}
    	//`得到重心`
    	Point getbarycentre(){
    		Point ret(0,0);
    		double area = 0;
    		for(int i = 1;i < n-1;i++){
    			double tmp = (p[i]-p[0])^(p[i+1]-p[0]);
    			if(sgn(tmp) == 0)continue;
    			area += tmp;
    			ret.x += (p[0].x+p[i].x+p[i+1].x)/3*tmp;
    			ret.y += (p[0].y+p[i].y+p[i+1].y)/3*tmp;
    		}
    		if(sgn(area)) ret = ret/area;
    		return ret;
    	}
    	//`多边形和圆交的面积`
    	//`测试:POJ3675 HDU3982 HDU2892`
    	double areacircle(circle c){
    		double ans = 0;
    		for(int i = 0;i < n;i++){
    			int j = (i+1)%n;
    			if(sgn( (p[j]-c.p)^(p[i]-c.p) ) >= 0)
    				ans += c.areatriangle(p[i],p[j]);
    			else ans -= c.areatriangle(p[i],p[j]);
    		}
    		return fabs(ans);
    	}
    	//`多边形和圆关系`
    	//` 2 圆完全在多边形内`
    	//` 1 圆在多边形里面,碰到了多边形边界`
    	//` 0 其它`
    	int relationcircle(circle c){
    		getline();
    		int x = 2;
    		if(relationpoint(c.p) != 1)return 0;//圆心不在内部
    		for(int i = 0;i < n;i++){
    			if(c.relationseg(l[i])==2)return 0;
    			if(c.relationseg(l[i])==1)x = 1;
    		}
    		return x;
    	}
    };
    //`AB X AC`
    double cross(Point A,Point B,Point C){
    	return (B-A)^(C-A);
    }
    //`AB*AC`
    double dot(Point A,Point B,Point C){
    	return (B-A)*(C-A);
    }
    //`最小矩形面积覆盖`
    //` A 必须是凸包(而且是逆时针顺序)`
    //` 测试 UVA 10173`
    double minRectangleCover(polygon A){
    	//`要特判A.n < 3的情况`
    	if(A.n < 3)return 0.0;
    	A.p[A.n] = A.p[0];
    	double ans = -1;
    	int r = 1, p = 1, q;
    	for(int i = 0;i < A.n;i++){
    		//`卡出离边A.p[i] - A.p[i+1]最远的点`
    		while( sgn( cross(A.p[i],A.p[i+1],A.p[r+1]) - cross(A.p[i],A.p[i+1],A.p[r]) ) >= 0 )
    			r = (r+1)%A.n;
    		//`卡出A.p[i] - A.p[i+1]方向上正向n最远的点`
    		while(sgn( dot(A.p[i],A.p[i+1],A.p[p+1]) - dot(A.p[i],A.p[i+1],A.p[p]) ) >= 0 )
    			p = (p+1)%A.n;
    		if(i == 0)q = p;
    		//`卡出A.p[i] - A.p[i+1]方向上负向最远的点`
    		while(sgn(dot(A.p[i],A.p[i+1],A.p[q+1]) - dot(A.p[i],A.p[i+1],A.p[q])) <= 0)
    			q = (q+1)%A.n;
    		double d = (A.p[i] - A.p[i+1]).len2();
    		double tmp = cross(A.p[i],A.p[i+1],A.p[r]) *
    			(dot(A.p[i],A.p[i+1],A.p[p]) - dot(A.p[i],A.p[i+1],A.p[q]))/d;
    		if(ans < 0 || ans > tmp)ans = tmp;
    	}
    	return ans;
    }
    
    //`直线切凸多边形`
    //`多边形是逆时针的,在q1q2的左侧`
    //`测试:HDU3982`
    vector<Point> convexCut(const vector<Point> &ps,Point q1,Point q2){
    	vector<Point>qs;
    	int n = ps.size();
    	for(int i = 0;i < n;i++){
    		Point p1 = ps[i], p2 = ps[(i+1)%n];
    		int d1 = sgn((q2-q1)^(p1-q1)), d2 = sgn((q2-q1)^(p2-q1));
    		if(d1 >= 0)
    			qs.push_back(p1);
    		if(d1 * d2 < 0)
    			qs.push_back(Line(p1,p2).crosspoint(Line(q1,q2)));
    	}
    	return qs;
    }
    //`半平面交`
    //`测试 POJ3335 POJ1474 POJ1279`
    //***************************
    struct halfplane:public Line{
    	double angle;
    	halfplane(){}
    	//`表示向量s->e逆时针(左侧)的半平面`
    	halfplane(Point _s,Point _e){
    		s = _s;
    		e = _e;
    	}
    	halfplane(Line v){
    		s = v.s;
    		e = v.e;
    	}
    	void calcangle(){
    		angle = atan2(e.y-s.y,e.x-s.x);
    	}
    	bool operator <(const halfplane &b)const{
    		return angle < b.angle;
    	}
    };
    struct halfplanes{
    	int n;
    	halfplane hp[maxp];
    	Point p[maxp];
    	int que[maxp];
    	int st,ed;
    	void push(halfplane tmp){
    		hp[n++] = tmp;
    	}
    	//去重
    	void unique(){
    		int m = 1;
    		for(int i = 1;i < n;i++){
    			if(sgn(hp[i].angle-hp[i-1].angle) != 0)
    				hp[m++] = hp[i];
    			else if(sgn( (hp[m-1].e-hp[m-1].s)^(hp[i].s-hp[m-1].s) ) > 0)
    				hp[m-1] = hp[i];
    		}
    		n = m;
    	}
    	bool halfplaneinsert(){
    		for(int i = 0;i < n;i++)hp[i].calcangle();
    		sort(hp,hp+n);
    		unique();
    		que[st=0] = 0;
    		que[ed=1] = 1;
    		p[1] = hp[0].crosspoint(hp[1]);
    		for(int i = 2;i < n;i++){
    			while(st<ed && sgn((hp[i].e-hp[i].s)^(p[ed]-hp[i].s))<0)ed--;
    			while(st<ed && sgn((hp[i].e-hp[i].s)^(p[st+1]-hp[i].s))<0)st++;
    			que[++ed] = i;
    			if(hp[i].parallel(hp[que[ed-1]]))return false;
    			p[ed]=hp[i].crosspoint(hp[que[ed-1]]);
    		}
    		while(st<ed && sgn((hp[que[st]].e-hp[que[st]].s)^(p[ed]-hp[que[st]].s))<0)ed--;
    		while(st<ed && sgn((hp[que[ed]].e-hp[que[ed]].s)^(p[st+1]-hp[que[ed]].s))<0)st++;
    		if(st+1>=ed)return false;
    		return true;
    	}
    	//`得到最后半平面交得到的凸多边形`
    	//`需要先调用halfplaneinsert() 且返回true`
    	void getconvex(polygon &con){
    		p[st] = hp[que[st]].crosspoint(hp[que[ed]]);
    		con.n = ed-st+1;
    		for(int j = st,i = 0;j <= ed;i++,j++)
    			con.p[i] = p[j];
    	}
    };
    //***************************
    
    const int maxn = 1010;
    struct circles{
    	circle c[maxn];
    	double ans[maxn];//`ans[i]表示被覆盖了i次的面积`
    	double pre[maxn];
    	int n;
    	circles(){}
    	void add(circle cc){
    		c[n++] = cc;
    	}
    	//`x包含在y中`
    	bool inner(circle x,circle y){
    		if(x.relationcircle(y) != 1)return 0;
    		return sgn(x.r-y.r)<=0?1:0;
    	}
    	//圆的面积并去掉内含的圆
    	void init_or(){
    		bool mark[maxn] = {0};
    		int i,j,k=0;
    		for(i = 0;i < n;i++){
    			for(j = 0;j < n;j++)
    				if(i != j && !mark[j]){
    					if( (c[i]==c[j])||inner(c[i],c[j]) )break;
    				}
    			if(j < n)mark[i] = 1;
    		}
    		for(i = 0;i < n;i++)
    			if(!mark[i])
    				c[k++] = c[i];
    		n = k;
    	}
    	//`圆的面积交去掉内含的圆`
    	void init_add(){
    		int i,j,k;
    		bool mark[maxn] = {0};
    		for(i = 0;i < n;i++){
    			for(j = 0;j < n;j++)
    				if(i != j && !mark[j]){
    					if( (c[i]==c[j])||inner(c[j],c[i]) )break;
    				}
    			if(j < n)mark[i] = 1;
    		}
    		for(i = 0;i < n;i++)
    			if(!mark[i])
    				c[k++] = c[i];
    		n = k;
    	}
    	//`半径为r的圆,弧度为th对应的弓形的面积`
    	double areaarc(double th,double r){
    		return 0.5*r*r*(th-sin(th));
    	}
    	//`测试SPOJVCIRCLES SPOJCIRUT`
    	//`SPOJVCIRCLES求n个圆并的面积,需要加上init\_or()去掉重复圆(否则WA)`
    	//`SPOJCIRUT 是求被覆盖k次的面积,不能加init\_or()`
    	//`对于求覆盖多少次面积的问题,不能解决相同圆,而且不能init\_or()`
    	//`求多圆面积并,需要init\_or,其中一个目的就是去掉相同圆`
    	void getarea(){
    		memset(ans,0,sizeof(ans));
    		vector<pair<double,int> >v;
    		for(int i = 0;i < n;i++){
    			v.clear();
    			v.push_back(make_pair(-pi,1));
    			v.push_back(make_pair(pi,-1));
    			for(int j = 0;j < n;j++)
    				if(i != j){
    					Point q = (c[j].p - c[i].p);
    					double ab = q.len(),ac = c[i].r, bc = c[j].r;
    					if(sgn(ab+ac-bc)<=0){
    						v.push_back(make_pair(-pi,1));
    						v.push_back(make_pair(pi,-1));
    						continue;
    					}
    					if(sgn(ab+bc-ac)<=0)continue;
    					if(sgn(ab-ac-bc)>0)continue;
    					double th = atan2(q.y,q.x), fai = acos((ac*ac+ab*ab-bc*bc)/(2.0*ac*ab));
    					double a0 = th-fai;
    					if(sgn(a0+pi)<0)a0+=2*pi;
    					double a1 = th+fai;
    					if(sgn(a1-pi)>0)a1-=2*pi;
    					if(sgn(a0-a1)>0){
    						v.push_back(make_pair(a0,1));
    						v.push_back(make_pair(pi,-1));
    						v.push_back(make_pair(-pi,1));
    						v.push_back(make_pair(a1,-1));
    					}
    					else{
    						v.push_back(make_pair(a0,1));
    						v.push_back(make_pair(a1,-1));
    					}
    				}
    			sort(v.begin(),v.end());
    			int cur = 0;
    			for(int j = 0;j < v.size();j++){
    				if(cur && sgn(v[j].first-pre[cur])){
    					ans[cur] += areaarc(v[j].first-pre[cur],c[i].r);
    					ans[cur] += 0.5*(Point(c[i].p.x+c[i].r*cos(pre[cur]),c[i].p.y+c[i].r*sin(pre[cur]))^Point(c[i].p.x+c[i].r*cos(v[j].first),c[i].p.y+c[i].r*sin(v[j].first)));
    				}
    				cur += v[j].second;
    				pre[cur] = v[j].first;
    			}
    		}
    		for(int i = 1;i < n;i++)
    			ans[i] -= ans[i+1];
    	}
    };

    三维几何

    const double eps = 1e-8;
    int sgn(double x){
    	if(fabs(x) < eps)return 0;
    	if(x < 0)return -1;
    	else return 1;
    }
    struct Point3{
    	double x,y,z;
    	Point3(double _x = 0,double _y = 0,double _z = 0){
    		x = _x;
    		y = _y;
    		z = _z;
    	}
    	void input(){
    		scanf("%lf%lf%lf",&x,&y,&z);
    	}
    	void output(){
    		printf("%.2lf %.2lf %.2lf\n",x,y,z);
    	}
    	bool operator ==(const Point3 &b)const{
    		return sgn(x-b.x) == 0 && sgn(y-b.y) == 0 && sgn(z-b.z) == 0;
    	}
    	bool operator <(const Point3 &b)const{
    		return sgn(x-b.x)==0?(sgn(y-b.y)==0?sgn(z-b.z)<0:y<b.y):x<b.x;
    	}
    	double len(){
    		return sqrt(x*x+y*y+z*z);
    	}
    	double len2(){
    		return x*x+y*y+z*z;
    	}
    	double distance(const Point3 &b)const{
    		return sqrt((x-b.x)*(x-b.x)+(y-b.y)*(y-b.y)+(z-b.z)*(z-b.z));
    	}
    	Point3 operator -(const Point3 &b)const{
    		return Point3(x-b.x,y-b.y,z-b.z);
    	}
    	Point3 operator +(const Point3 &b)const{
    		return Point3(x+b.x,y+b.y,z+b.z);
    	}
    	Point3 operator *(const double &k)const{
    		return Point3(x*k,y*k,z*k);
    	}
    	Point3 operator /(const double &k)const{
    		return Point3(x/k,y/k,z/k);
    	}
    	//点乘
    	double operator *(const Point3 &b)const{
    		return x*b.x+y*b.y+z*b.z;
    	}
    	//叉乘
    	Point3 operator ^(const Point3 &b)const{
    		return Point3(y*b.z-z*b.y,z*b.x-x*b.z,x*b.y-y*b.x);
    	}
    	double rad(Point3 a,Point3 b){
    		Point3 p = (*this);
    		return acos( ( (a-p)*(b-p) )/ (a.distance(p)*b.distance(p)) );
    	}
    	//变换长度
    	Point3 trunc(double r){
    		double l = len();
    		if(!sgn(l))return *this;
    		r /= l;
    		return Point3(x*r,y*r,z*r);
    	}
    };
    struct Line3
    {
    	Point3 s,e;
    	Line3(){}
    	Line3(Point3 _s,Point3 _e)
    	{
    		s = _s;
    		e = _e;
    	}
    	bool operator ==(const Line3 v)
    	{
    		return (s==v.s)&&(e==v.e);
    	}
    	void input()
    	{
    		s.input();
    		e.input();
    	}
    	double length()
    	{
    		return s.distance(e);
    	}
    	//点到直线距离
    	double dispointtoline(Point3 p)
    	{
    		return ((e-s)^(p-s)).len()/s.distance(e);
    	}
    	//点到线段距离
    	double dispointtoseg(Point3 p)
    	{
    		if(sgn((p-s)*(e-s)) < 0 || sgn((p-e)*(s-e)) < 0)
    			return min(p.distance(s),e.distance(p));
    		return dispointtoline(p);
    	}
    	//`返回点p在直线上的投影`
    	Point3 lineprog(Point3 p)
    	{
    		return s + ( ((e-s)*((e-s)*(p-s)))/((e-s).len2()) );
    	}
    	//`p绕此向量逆时针arg角度`
    	Point3 rotate(Point3 p,double ang)
    	{
    		if(sgn(((s-p)^(e-p)).len()) == 0)return p;
    		Point3 f1 = (e-s)^(p-s);
    		Point3 f2 = (e-s)^(f1);
    		double len = ((s-p)^(e-p)).len()/s.distance(e);
    		f1 = f1.trunc(len); f2 = f2.trunc(len);
    		Point3 h = p+f2;
    		Point3 pp = h+f1;
    		return h + ((p-h)*cos(ang)) + ((pp-h)*sin(ang));
    	}
    	//`点在直线上`
    	bool pointonseg(Point3 p)
    	{
    		return sgn( ((s-p)^(e-p)).len() ) == 0 && sgn((s-p)*(e-p)) == 0;
    	}
    };
    struct Plane
    {
    	Point3 a,b,c,o;//`平面上的三个点,以及法向量`
    	Plane(){}
    	Plane(Point3 _a,Point3 _b,Point3 _c)
    	{
    		a = _a;
    		b = _b;
    		c = _c;
    		o = pvec();
    	}
    	Point3 pvec()
    	{
    		return (b-a)^(c-a);
    	}
    	//`ax+by+cz+d = 0`
    	Plane(double _a,double _b,double _c,double _d)
    	{
    		o = Point3(_a,_b,_c);
    		if(sgn(_a) != 0)
    			a = Point3((-_d-_c-_b)/_a,1,1);
    		else if(sgn(_b) != 0)
    			a = Point3(1,(-_d-_c-_a)/_b,1);
    		else if(sgn(_c) != 0)
    			a = Point3(1,1,(-_d-_a-_b)/_c);
    	}
    	//`点在平面上的判断`
    	bool pointonplane(Point3 p)
    	{
    		return sgn((p-a)*o) == 0;
    	}
    	//`两平面夹角`
    	double angleplane(Plane f)
    	{
    		return acos(o*f.o)/(o.len()*f.o.len());
    	}
    	//`平面和直线的交点,返回值是交点个数`
    	int crossline(Line3 u,Point3 &p)
    	{
    		double x = o*(u.e-a);
    		double y = o*(u.s-a);
    		double d = x-y;
    		if(sgn(d) == 0)return 0;
    		p = ((u.s*x)-(u.e*y))/d;
    		return 1;
    	}
    	//`点到平面最近点(也就是投影)`
    	Point3 pointtoplane(Point3 p)
    	{
    		Line3 u = Line3(p,p+o);
    		crossline(u,p);
    		return p;
    	}
    	//`平面和平面的交线`
    	int crossplane(Plane f,Line3 &u)
    	{
    		Point3 oo = o^f.o;
    		Point3 v = o^oo;
    		double d = fabs(f.o*v);
    		if(sgn(d) == 0)return 0;
    		Point3 q = a + (v*(f.o*(f.a-a))/d);
    		u = Line3(q,q+oo);
    		return 1;
    	}
    };

     平面最近点对

    const int MAXN = 100010;
    const double eps = 1e-8;
    const double INF = 1e20;
    struct Point{
    	double x,y;
    	void input(){
    		scanf("%lf%lf",&x,&y);
    	}
    };
    double dist(Point a,Point b){
    	return sqrt((a.x-b.x)*(a.x-b.x) + (a.y-b.y)*(a.y-b.y));
    }
    Point p[MAXN];
    Point tmpt[MAXN];
    bool cmpx(Point a,Point b){
    	return a.x < b.x || (a.x == b.x && a.y < b.y);
    }
    bool cmpy(Point a,Point b){
    	return a.y < b.y || (a.y == b.y && a.x < b.x);
    }
    double Closest_Pair(int left,int right){
    	double d = INF;
    	if(left == right)return d;
    	if(left+1 == right)return dist(p[left],p[right]);
    	int mid = (left+right)/2;
    	double d1 = Closest_Pair(left,mid);
    	double d2 = Closest_Pair(mid+1,right);
    	d = min(d1,d2);
    	int cnt = 0;
    	for(int i = left;i <= right;i++){
    		if(fabs(p[mid].x - p[i].x) <= d)
    			tmpt[cnt++] = p[i];
    	}
    	sort(tmpt,tmpt+cnt,cmpy);
    	for(int i = 0;i < cnt;i++){
    		for(int j = i+1;j < cnt && tmpt[j].y - tmpt[i].y < d;j++)
    			d = min(d,dist(tmpt[i],tmpt[j]));
    	}
    	return d;
    }
    int main(){
    	int n;
    	while(scanf("%d",&n) == 1 && n){
    		for(int i = 0;i < n;i++)p[i].input();
    		sort(p,p+n,cmpx);
    		printf("%.2lf\n",Closest_Pair(0,n-1));
    	}
        return 0;
    }

    三维凸包:HDU4273

    const double eps = 1e-8;
    const int MAXN = 550;
    int sgn(double x){
    	if(fabs(x) < eps)return 0;
    	if(x < 0)return -1;
    	else return 1;
    }
    struct Point3{
    	double x,y,z;
    	Point3(double _x = 0, double _y = 0, double _z = 0){
    		x = _x;
    		y = _y;
    		z = _z;
    	}
    	void input(){
    		scanf("%lf%lf%lf",&x,&y,&z);
    	}
    	bool operator ==(const Point3 &b)const{
    		return sgn(x-b.x) == 0 && sgn(y-b.y) == 0 && sgn(z-b.z) == 0;
    	}
    	double len(){
    		return sqrt(x*x+y*y+z*z);
    	}
    	double len2(){
    		return x*x+y*y+z*z;
    	}
    	double distance(const Point3 &b)const{
    		return sqrt((x-b.x)*(x-b.x)+(y-b.y)*(y-b.y)+(z-b.z)*(z-b.z));
    	}
    	Point3 operator -(const Point3 &b)const{
    		return Point3(x-b.x,y-b.y,z-b.z);
    	}
    	Point3 operator +(const Point3 &b)const{
    		return Point3(x+b.x,y+b.y,z+b.z);
    	}
    	Point3 operator *(const double &k)const{
    		return Point3(x*k,y*k,z*k);
    	}
    	Point3 operator /(const double &k)const{
    		return Point3(x/k,y/k,z/k);
    	}
    	//点乘
    	double operator *(const Point3 &b)const{
    		return x*b.x + y*b.y + z*b.z;
    	}
    	//叉乘
    	Point3 operator ^(const Point3 &b)const{
    		return Point3(y*b.z-z*b.y,z*b.x-x*b.z,x*b.y-y*b.x);
    	}
    };
    struct CH3D{
    	struct face{
    		//表示凸包一个面上的三个点的编号
    		int a,b,c;
    		//表示该面是否属于最终的凸包上的面
    		bool ok;
    	};
    	//初始顶点数
    	int n;
    	Point3 P[MAXN];
    	//凸包表面的三角形数
    	int num;
    	//凸包表面的三角形
    	face F[8*MAXN];
    	int g[MAXN][MAXN];
    	//叉乘
    	Point3 cross(const Point3 &a,const Point3 &b,const Point3 &c){
    		return (b-a)^(c-a);
    	}
    	//`三角形面积*2`
    	double area(Point3 a,Point3 b,Point3 c){
    		return ((b-a)^(c-a)).len();
    	}
    	//`四面体有向面积*6`
    	double volume(Point3 a,Point3 b,Point3 c,Point3 d){
    		return ((b-a)^(c-a))*(d-a);
    	}
    	//`正:点在面同向`
    	double dblcmp(Point3 &p,face &f){
    		Point3 p1 = P[f.b] - P[f.a];
    		Point3 p2 = P[f.c] - P[f.a];
    		Point3 p3 = p - P[f.a];
    		return (p1^p2)*p3;
    	}
    	void deal(int p,int a,int b){
    		int f = g[a][b];
    		face add;
    		if(F[f].ok){
    			if(dblcmp(P[p],F[f]) > eps)
    				dfs(p,f);
    			else {
    				add.a = b;
    				add.b = a;
    				add.c = p;
    				add.ok = true;
    				g[p][b] = g[a][p] = g[b][a] = num;
    				F[num++] = add;
    			}
    		}
    	}
    	//递归搜索所有应该从凸包内删除的面
    	void dfs(int p,int now){
    		F[now].ok = false;
    		deal(p,F[now].b,F[now].a);
    		deal(p,F[now].c,F[now].b);
    		deal(p,F[now].a,F[now].c);
    	}
    	bool same(int s,int t){
    		Point3 &a = P[F[s].a];
    		Point3 &b = P[F[s].b];
    		Point3 &c = P[F[s].c];
    		return fabs(volume(a,b,c,P[F[t].a])) < eps &&
    			fabs(volume(a,b,c,P[F[t].b])) < eps &&
    			fabs(volume(a,b,c,P[F[t].c])) < eps;
    	}
    	//构建三维凸包
    	void create(){
    		num = 0;
    		face add;
    
    		//***********************************
    		//此段是为了保证前四个点不共面
    		bool flag = true;
    		for(int i = 1;i < n;i++){
    			if(!(P[0] == P[i])){
    				swap(P[1],P[i]);
    				flag = false;
    				break;
    			}
    		}
    		if(flag)return;
    		flag = true;
    		for(int i = 2;i < n;i++){
    			if( ((P[1]-P[0])^(P[i]-P[0])).len() > eps ){
    				swap(P[2],P[i]);
    				flag = false;
    				break;
    			}
    		}
    		if(flag)return;
    		flag = true;
    		for(int i = 3;i < n;i++){
    			if(fabs( ((P[1]-P[0])^(P[2]-P[0]))*(P[i]-P[0]) ) > eps){
    				swap(P[3],P[i]);
    				flag = false;
    				break;
    			}
    		}
    		if(flag)return;
    		//**********************************
    
    		for(int i = 0;i < 4;i++){
    			add.a = (i+1)%4;
    			add.b = (i+2)%4;
    			add.c = (i+3)%4;
    			add.ok = true;
    			if(dblcmp(P[i],add) > 0)swap(add.b,add.c);
    			g[add.a][add.b] = g[add.b][add.c] = g[add.c][add.a] = num;
    			F[num++] = add;
    		}
    		for(int i = 4;i < n;i++)
    			for(int j = 0;j < num;j++)
    				if(F[j].ok && dblcmp(P[i],F[j]) > eps){
    					dfs(i,j);
    					break;
    				}
    		int tmp = num;
    		num = 0;
    		for(int i = 0;i < tmp;i++)
    			if(F[i].ok)
    				F[num++] = F[i];
    	}
    	//表面积
    	//`测试:HDU3528`
    	double area(){
    		double res = 0;
    		if(n == 3){
    			Point3 p = cross(P[0],P[1],P[2]);
    			return p.len()/2;
    		}
    		for(int i = 0;i < num;i++)
    			res += area(P[F[i].a],P[F[i].b],P[F[i].c]);
    		return res/2.0;
    	}
    	double volume(){
    		double res = 0;
    		Point3 tmp = Point3(0,0,0);
    		for(int i = 0;i < num;i++)
    			res += volume(tmp,P[F[i].a],P[F[i].b],P[F[i].c]);
    		return fabs(res/6);
    	}
    	//表面三角形个数
    	int triangle(){
    		return num;
    	}
    	//表面多边形个数
    	//`测试:HDU3662`
    	int polygon(){
    		int res = 0;
    		for(int i = 0;i < num;i++){
    			bool flag = true;
    			for(int j = 0;j < i;j++)
    				if(same(i,j)){
    					flag = 0;
    					break;
    				}
    			res += flag;
    		}
    		return res;
    	}
    	//重心
    	//`测试:HDU4273`
    	Point3 barycenter(){
    		Point3 ans = Point3(0,0,0);
    		Point3 o = Point3(0,0,0);
    		double all = 0;
    		for(int i = 0;i < num;i++){
    			double vol = volume(o,P[F[i].a],P[F[i].b],P[F[i].c]);
    			ans = ans + (((o+P[F[i].a]+P[F[i].b]+P[F[i].c])/4.0)*vol);
    			all += vol;
    		}
    		ans = ans/all;
    		return ans;
    	}
    	//点到面的距离
    	//`测试:HDU4273`
    	double ptoface(Point3 p,int i){
    		double tmp1 = fabs(volume(P[F[i].a],P[F[i].b],P[F[i].c],p));
    		double tmp2 = ((P[F[i].b]-P[F[i].a])^(P[F[i].c]-P[F[i].a])).len();
    		return tmp1/tmp2;
    	}
    };
    CH3D hull;
    int main()
    {
        while(scanf("%d",&hull.n) == 1){
    		for(int i = 0;i < hull.n;i++)hull.P[i].input();
    		hull.create();
    		Point3 p = hull.barycenter();
    		double ans = 1e20;
    		for(int i = 0;i < hull.num;i++)
    			ans = min(ans,hull.ptoface(p,i));
    		printf("%.3lf\n",ans);
    	}
        return 0;
    }

    叉乘性质:

    展开全文
  • // `计算几何模板` const double eps = 1e-8; const double inf = 1e20; const double pi = acos(-1.0); const int maxp = 1010; //`Compares a double to zero` int sgn(double x){ if(fabs(x) < eps)re...
  • acm计算几何模板

    2019-05-09 15:43:08
    不要过分相信的我抄的,我在用的时候也经常会遇到错,所以如果不对,还请自己改,或者自己找一份原kuangbin大佬的计算几何板子。 # include using namespace std ; # define me(x,y) memset(x,y,sizeof x)...
  • 这道题目就是单纯的计算几何板子,我照着板子自己敲了一遍,希望大家不熟悉板子的也可以自己敲一下。 两点: 要选择适合自己的板子 注意最后结果要排序 题目要求的六种结果: 三角形的外接圆 三角形的内切圆 过...
  • } 2.2叉积的运用: 2.2.1计算多边形的面积 2.2.2判断凹凸多边形 2.2.3折线的拐向的方向 2.2.4:判断点是否在线段上 代码模板 struct line//线段 { Point s; Point e; line(Point s, Point e){s = a, e = b;} line(){...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 1,659
精华内容 663
关键字:

计算几何板子