zoukankan      html  css  js  c++  java
  • 计算几何模板

    计算几何基础模板

    基础类型

    const double eps = 1e-8;
    const double PI=acos(-1.0);
    int sgn(double x) {
        if(fabs(x) < eps) return 0;
        if(x < 0) return -1;
        return 1;
    }
    struct Point {
        double x,y;
        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);
        }
        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;
        }
        Point operator *(double b)const {
            return Point(x*b,y*b);
        }
        Point Rotate(double rad){ //逆时针旋转
            return Point(x*cos(rad)-y*sin(rad),x*sin(rad)+y*cos(rad));
        }
        double angle(){
            return atan2(y,x);
        }
        double len(){//Vector
            return sqrt(x*x+y*y);
        }
        void stdd(){//Vector
            double le=len();
            x/=le;y/=le;
        }
        void input(){
            scanf("%lf%lf",&x,&y);
        }
        void output(){
            printf("(%.8f,%.8f) ",x,y);
        }
    };
    double xmult(Point p0,Point p1,Point p2) { //p0p1 X p0p2
        return (p1-p0)^(p2-p0);
    }
    double dist(Point a,Point b) {//a到b距离
        return sqrt( (b - a)*(b - a) );
    }
    double getArea(Point a,Point b,Point c){//三角形面积S
        return fabs(xmult(a,b,c))/2.0;
    }
    double distLP(Point a,Point b,Point p){//p到ab的距离
        return getArea(a,b,p)*2/dist(a,b);
    }
    double calArea(vector<Point> &p){//按逆时针排列的多边形面积
        double ans=0;
        int m=p.size();
        for(int i=0;i<m;i++){
            ans+=p[i]^p[(i+1)%m];
        }
        return fabs(ans/2.0);
    }
    struct Line{
    	Point s,t;
    	double ang;
    	Line(Point X=Point(),Point Y=Point()){
    		s=X,t=Y,ang=(Y-X).angle();
    	}
    	double getangle(){
            return ang=(t-s).angle();
    	}
        double getDistance(Point A){ //点到直线的距离
            return fabs((A-s)^(A-t))/dist(s,t);
        }
    };
    Point getIntersectPoint(Line a, Line b) {//两直线交,注意如果两直线重合会出错,要先判断
        double a1 = a.s.y - a.t.y, b1 = a.t.x - a.s.x, c1 = a.s.x * a.t.y - a.t.x * a.s.y;
        double a2 = b.s.y - b.t.y, b2 = b.t.x - b.s.x, c2 = b.s.x * b.t.y - b.t.x * b.s.y;
        return Point((c1*b2-c2*b1)/(a2*b1-a1*b2), (a2*c1-a1*c2)/(a1*b2-a2*b1));
    }
    struct Circle {
        Point o;
        double r;
        Circle(){}
        void input(){
            o.input();
            scanf("%lf",&r);
        }
        void output(){
            printf("%.8f %.8f %.8f
    ",o.x,o.y,r);
        }
    };
    

    判断两个线段是否相交

    bool jiao(Point a,Point b,Point c,Point d){//判断线段ab和cd是否相交
        if( max(a.x,b.x)<min(c.x,d.x)||
            max(a.y,b.y)<min(c.y,d.y)||
            max(c.x,d.x)<min(a.x,b.x)||
            max(c.y,d.y)<min(a.y,b.y)||
            xmult(a,c,b)*xmult(a,b,d)<0||
            xmult(c,a,d)*xmult(c,d,b)<0 )
        return 0;
        else return 1;
    }
    

    atan2

    atan2(double y,double x);

    返回坐标(x,y)的极角(-pi~pi)

    第一象限为((0,pi/2))

    第二象限为((pi/2,pi))

    第三象限为 ((-pi,-pi/2))

    第四象限为((-pi/2,0))

    极角排序

    选用一个点o为中心进行极角排序时,需要先将其他点坐标转化为相对o的坐标,即原坐标-o的坐标。

    用atan2直接排序,存在精度问题

    bool cmp1(Point a,Point b){
        if(atan2(a.y,a.x)!=atan2(b.y,b.x))
            return atan2(a.y,a.x)<atan2(b.y,b.x);
        else return a.x<b.x;
    }
    

    利用叉积排序,要求点分布在180°圆心角内

    bool cmp2(Point a,Point b){
        return sgn(a^b)>0;
    }
    

    先按照象限排序,再按照叉积排序:

    int getq(Point a) {
        if(a.x>0 && a.y>=0) return 1;
        if(a.x<=0 && a.y>0) return 2;
        if(a.x<0 && a.y<=0) return 3;
        if(a.x>=0 && a.y<0) return 4;
    }
    bool cmp(Point a,Point b){
        if(getq(a)==getq(b))
            return cmp2(a,b);
        else return getq(a)<getq(b);
    }
    

    凸包

    bool cmpA(Point a,Point b){
        return a.x<b.x||(a.x==b.x && a.y<b.y);
    }
    vector<Point> Andrew(vector<Point> p) {//输入不能有重复点,若要凸包边上没有输入点,将两个<=改为<
        sort(p.begin(),p.end(),cmpA);
        vector<Point>tb;
        for(int i=0; i<p.size(); i++) {
            while(tb.size()>=2&&sgn(xmult(tb[tb.size()-2],tb[tb.size()-1],p[i]))<=0)tb.pop_back();
            tb.push_back(p[i]);
        }
        int temp=tb.size();
        for(int i=p.size()-2; i>=0; i--) {
            while(tb.size()>temp&&sgn(xmult(tb[tb.size()-2],tb[tb.size()-1],p[i]))<=0)tb.pop_back();
            tb.push_back(p[i]);
        }
        if(p.size()>1)tb.pop_back();
        return tb;
    }
    

    半平面交

    bool onRight(Line a,Line b,Line c){
        Point jiao=getIntersectPoint(b,c);
        if(xmult(a.s,a.t,jiao)<0){
            return 1;
        }
        else{
            return 0;
        }
    }
    bool cmpHL(Line a,Line b){
        double A=a.getangle(),B=b.getangle();
        if(sgn(A-B)==0){//平行的直线将最左边的放后面,便于去重
            return xmult(a.s,a.t,b.t)>=0;
        }
        else{
            return A<B;
        }
    }
    vector<Line> getHL(vector<Line> l){
        //去除角度相同的,保留最最左的
        sort(l.begin(),l.end(),cmpHL);
        int n=l.size();
        int cnt=0;
        for(int i=0;i<=n-2;i++){
            if(sgn(l[i].getangle()-l[i+1].getangle())==0){
                continue;
            }
            l[cnt++]=l[i];
        }
        l[cnt++]=l[n-1];
    
        deque<Line> que;
        for(int i=0;i<cnt;i++){
            while(que.size()>=2&&onRight(l[i],que[que.size()-1],que[que.size()-2])) que.pop_back();
            while(que.size()>=2&&onRight(l[i],que[0],que[1])) que.pop_front();
            que.push_back(l[i]);
        }
        while(que.size()>=3&&onRight(que[0],que[que.size()-1],que[que.size()-2])) que.pop_back();
        while(que.size()>=3&&onRight(que[que.size()-1],que[0],que[1])) que.pop_front();
    
        vector<Line> hl;
        for(int i=0;i<que.size();i++){
            hl.push_back(que[i]);
        }
        return hl;
    }
    

    旋转卡壳

    //求凸包的直径
    int getMax(vector<Point> p){//要求凸包逆时针排列
        int n=p.size();
        int ans=0;
        if(n==2){
            return dist(p[1],p[0]);
        }
        int j=2;
        for(int i=0;i<n;i++){
            while(getArea(p[i],p[(i+1)%n],p[j])<getArea(p[i],p[(i+1)%n],p[(j+1)%n])){
                j=(j+1)%n;
            }
            ans=max(ans,max(dist(p[i],p[j]),dist(p[(i+1)%n],p[j])));
        }
        return ans;
    }
    

    反演变换

    Point PtP(Point a,Point p,double r){//点到点
        Point v1=a-p;
        v1.stdd();
        double len=r*r/dist(a,p);
        return p+v1*len;
    }
    Circle CtC(Circle C,Point p,double r){//圆到圆
        Circle res;
        double t = dist(C.o,p);
        double x = r*r / (t - C.r);
        double y = r*r / (t + C.r);
        res.r = (x - y) / 2.0;
        double s = (x + y) / 2.0;
        res.o = p + (C.o - p) * (s / t);
        return res;
    }
    Circle LtC(Point a,Point b,Point p,double r){//直线到过反演点的圆
        double d=distLP(a,b,p);
        d=r*r/d;
        Circle c;
        c.r=d/2;
        Point v1;
        if(xmult(a,b,p)>0)
            v1=(a-b).Rotate(PI/2);
        else
            v1=(b-a).Rotate(PI/2);
        v1.stdd();
        c.o=p+v1*c.r;
        return c;
    }
    
    

    最小圆覆盖

    Point circumcenter(Point a,Point b,Point c){
        double x1=a.x,y1=a.y,x2=b.x,y2=b.y,x3=c.x,y3=c.y;
        double a1=x2-x1,b1=y2-y1,c1=(x2*x2-x1*x1+y2*y2-y1*y1)/2;
        double a2=x3-x1,b2=y3-y1,c2=(x3*x3-x1*x1+y3*y3-y1*y1)/2;
        return {(b2*c1-b1*c2)/(a1*b2-a2*b1),(a1*c2-a2*c1)/(a1*b2-a2*b1)};
    }
    void min_cover_circle(vector<Point> p,Point &c,double &r)
    {
        random_shuffle(p.begin(),p.end());      //将n个点随机打乱
        int n=p.size();
        c=p[0]; r=0;
        for(int i=1;i<n;i++)
        {
            if(dist(p[i],c)>r+eps)   //第一个点
            {
                c=p[i]; r=0;
                for(int j=0;j<i;j++)
                    if(dist(p[j],c)>r+eps)  //第二个点
                    {
                        c.x=(p[i].x+p[j].x)/2;
                        c.y=(p[i].y+p[j].y)/2;
                        r=dist(p[j],c);
                        for(int k=0;k<j;k++)
                            if(dist(p[k],c)>r+eps)  //第三个点
                            {   //求外接圆圆心,三点必不共线
                                c=circumcenter(p[i],p[j],p[k]);
                                r=dist(p[i],c);
                            }
                    }
            }
        }
    }
    

    未整理

    #include <iostream>
    #include <cstdio>
    #include <cmath>
    #include <algorithm>
    
    using namespace std;
    const double PI = acos(-1.0);
    const double eps = 1e-10;
    
    /****************常用函数***************/
    //判断ta与tb的大小关系
    int sgn( double ta, double tb)
    {
        if(fabs(ta-tb)<eps)return 0;
        if(ta<tb)   return -1;
        return 1;
    }
    
    //点
    class Point
    {
    public:
    
        double x, y;
    
        Point(){}
        Point( double tx, double ty){ x = tx, y = ty;}
    
        bool operator < (const Point &_se) const
        {
            return x<_se.x || (x==_se.x && y<_se.y);
        }
        friend Point operator + (const Point &_st,const Point &_se)
        {
            return Point(_st.x + _se.x, _st.y + _se.y);
        }
        friend Point operator - (const Point &_st,const Point &_se)
        {
            return Point(_st.x - _se.x, _st.y - _se.y);
        }
        //点位置相同(double类型)
        bool operator == (const Point &_off)const
        {
            return  sgn(x, _off.x) == 0 && sgn(y, _off.y) == 0;
        }
    
    };
    
    /****************常用函数***************/
    //点乘
    double dot(const Point &po,const Point &ps,const Point &pe)
    {
        return (ps.x - po.x) * (pe.x - po.x) + (ps.y - po.y) * (pe.y - po.y);
    }
    //叉乘
    double xmult(const Point &po,const Point &ps,const Point &pe)
    {
        return (ps.x - po.x) * (pe.y - po.y) - (pe.x - po.x) * (ps.y - po.y);
    }
    //两点间距离的平方
    double getdis2(const Point &st,const Point &se)
    {
        return (st.x - se.x) * (st.x - se.x) + (st.y - se.y) * (st.y - se.y);
    }
    //两点间距离
    double getdis(const Point &st,const Point &se)
    {
        return sqrt((st.x - se.x) * (st.x - se.x) + (st.y - se.y) * (st.y - se.y));
    }
    
    //两点表示的向量
    class Line
    {
    public:
    
        Point s, e;//两点表示,起点[s],终点[e]
        double a, b, c;//一般式,ax+by+c=0
        double angle;//向量的角度,[-pi,pi]
    
        Line(){}
        Line( Point ts, Point te):s(ts),e(te){}//get_angle();}
        Line(double _a,double _b,double _c):a(_a),b(_b),c(_c){}
    
        //排序用
        bool operator < (const Line &ta)const
        {
            return angle<ta.angle;
        }
        //向量与向量的叉乘
        friend double operator / ( const Line &_st, const  Line &_se)
        {
            return (_st.e.x - _st.s.x) * (_se.e.y - _se.s.y) - (_st.e.y - _st.s.y) * (_se.e.x - _se.s.x);
        }
        //向量间的点乘
        friend double operator *( const Line &_st, const  Line &_se)
        {
            return (_st.e.x - _st.s.x) * (_se.e.x - _se.s.x) - (_st.e.y - _st.s.y) * (_se.e.y - _se.s.y);
        }
        //从两点表示转换为一般表示
        //a=y2-y1,b=x1-x2,c=x2*y1-x1*y2
        bool pton()
        {
            a = e.y - s.y;
            b = s.x - e.x;
            c = e.x * s.y - e.y * s.x;
            return true;
        }
        //半平面交用
        //点在向量左边(右边的小于号改成大于号即可,在对应直线上则加上=号)
        friend bool operator < (const Point &_Off, const Line &_Ori)
        {
            return (_Ori.e.y - _Ori.s.y) * (_Off.x - _Ori.s.x)
                < (_Off.y - _Ori.s.y) * (_Ori.e.x - _Ori.s.x);
        }
        //求直线或向量的角度
        double get_angle( bool isVector = true)
        {
            angle = atan2( e.y - s.y, e.x - s.x);
            if(!isVector && angle < 0)
                angle += PI;
            return angle;
        }
    
        //点在线段或直线上 1:点在直线上 2点在s,e所在矩形内
        bool has(const Point &_Off, bool isSegment = false) const
        {
            bool ff = sgn( xmult( s, e, _Off), 0) == 0;
            if( !isSegment) return ff;
            return ff
                && sgn(_Off.x - min(s.x, e.x), 0) >= 0 && sgn(_Off.x - max(s.x, e.x), 0) <= 0
                && sgn(_Off.y - min(s.y, e.y), 0) >= 0 && sgn(_Off.y - max(s.y, e.y), 0) <= 0;
        }
    
        //点到直线/线段的距离
        double dis(const Point &_Off, bool isSegment = false)
        {
            ///化为一般式
            pton();
            //到直线垂足的距离
            double td = (a * _Off.x + b * _Off.y + c) / sqrt(a * a + b * b);
            //如果是线段判断垂足
            if(isSegment)
            {
                double xp = (b * b * _Off.x - a * b * _Off.y - a * c) / ( a * a + b * b);
                double yp = (-a * b * _Off.x + a * a * _Off.y - b * c) / (a * a + b * b);
                double xb = max(s.x, e.x);
                double yb = max(s.y, e.y);
                double xs = s.x + e.x - xb;
                double ys = s.y + e.y - yb;
               if(xp > xb + eps || xp < xs - eps || yp > yb + eps || yp < ys - eps)
                    td = min( getdis(_Off,s), getdis(_Off,e));
            }
            return fabs(td);
        }
    
        //关于直线对称的点
        Point mirror(const Point &_Off)
        {
            ///注意先转为一般式
            Point ret;
            double d = a * a + b * b;
            ret.x = (b * b * _Off.x - a * a * _Off.x - 2 * a * b * _Off.y - 2 * a * c) / d;
            ret.y = (a * a * _Off.y - b * b * _Off.y - 2 * a * b * _Off.x - 2 * b * c) / d;
            return ret;
        }
        //计算两点的中垂线
        static Line ppline(const Point &_a,const Point &_b)
        {
            Line ret;
            ret.s.x = (_a.x + _b.x) / 2;
            ret.s.y = (_a.y + _b.y) / 2;
            //一般式
            ret.a = _b.x - _a.x;
            ret.b = _b.y - _a.y;
            ret.c = (_a.y - _b.y) * ret.s.y + (_a.x - _b.x) * ret.s.x;
            //两点式
            if(fabs(ret.a) > eps)
            {
                ret.e.y = 0.0;
                ret.e.x = - ret.c / ret.a;
                if(ret.e == ret. s)
                {
                    ret.e.y = 1e10;
                    ret.e.x = - (ret.c - ret.b * ret.e.y) / ret.a;
                }
            }
            else
            {
                ret.e.x = 0.0;
                ret.e.y = - ret.c / ret.b;
                if(ret.e == ret. s)
                {
                    ret.e.x = 1e10;
                    ret.e.y = - (ret.c - ret.a * ret.e.x) / ret.b;
                }
            }
            return ret;
        }
    
        //------------直线和直线(向量)-------------
        //向量向左边平移t的距离
        Line& moveLine( double t)
        {
            Point of;
            of = Point( -( e.y - s.y), e.x - s.x);
            double dis = sqrt( of.x * of.x + of.y * of.y);
            of.x= of.x * t / dis, of.y = of.y * t / dis;
            s = s + of, e = e + of;
            return *this;
        }
        //直线重合
        static bool equal(const Line &_st,const Line &_se)
        {
            return _st.has( _se.e) && _se.has( _st.s);
        }
        //直线平行
        static bool parallel(const Line &_st,const Line &_se)
        {
            return sgn( _st / _se, 0) == 0;
        }
        //两直线(线段)交点
        //返回-1代表平行,0代表重合,1代表相交
        static bool crossLPt(const Line &_st,const Line &_se, Point &ret)
        {
            if(parallel(_st,_se))
            {
                if(Line::equal(_st,_se)) return 0;
                return -1;
            }
            ret = _st.s;
            double t = ( Line(_st.s,_se.s) / _se) / ( _st / _se);
            ret.x += (_st.e.x - _st.s.x) * t;
            ret.y += (_st.e.y - _st.s.y) * t;
            return 1;
        }
        //------------线段和直线(向量)----------
        //直线和线段相交
        //参数:直线[_st],线段[_se]
        friend bool crossSL( Line &_st, Line &_se)
        {
            return sgn( xmult( _st.s, _se.s, _st.e) * xmult( _st.s, _st.e, _se.e), 0) >= 0;
        }
    
        //判断线段是否相交(注意添加eps)
        static bool isCrossSS( const Line &_st, const  Line &_se)
        {
            //1.快速排斥试验判断以两条线段为对角线的两个矩形是否相交
            //2.跨立试验(等于0时端点重合)
            return
                max(_st.s.x, _st.e.x) >= min(_se.s.x, _se.e.x) &&
                max(_se.s.x, _se.e.x) >= min(_st.s.x, _st.e.x) &&
                max(_st.s.y, _st.e.y) >= min(_se.s.y, _se.e.y) &&
                max(_se.s.y, _se.e.y) >= min(_st.s.y, _st.e.y) &&
                sgn( xmult( _se.s, _st.s, _se.e) * xmult( _se.s, _se.e, _st.s), 0) >= 0 &&
                sgn( xmult( _st.s, _se.s, _st.e) * xmult( _st.s, _st.e, _se.s), 0) >= 0;
        }
    };
    
    //寻找凸包的graham 扫描法所需的排序函数
    Point gsort;
    bool gcmp( const Point &ta, const Point &tb)/// 选取与最后一条确定边夹角最小的点,即余弦值最大者
    {
        double tmp = xmult( gsort, ta, tb);
        if( fabs( tmp) < eps)
            return getdis( gsort, ta) < getdis( gsort, tb);
        else if( tmp > 0)
            return 1;
        return 0;
    }
    
    class Polygon
    {
    public:
        const static int maxpn = 5e4+7;
        Point pt[maxpn];//点(顺时针或逆时针)
        Line dq[maxpn]; //求半平面交打开注释
        int n;//点的个数
    
    
        //求多边形面积,多边形内点必须顺时针或逆时针
        double area()
        {
            double ans = 0.0;
            for(int i = 0; i < n; i ++)
            {
                int nt = (i + 1) % n;
                ans += pt[i].x * pt[nt].y - pt[nt].x * pt[i].y;
            }
            return fabs( ans / 2.0);
        }
        //求多边形重心,多边形内点必须顺时针或逆时针
        Point gravity()
        {
            Point ans;
            ans.x = ans.y = 0.0;
            double area = 0.0;
            for(int i = 0; i < n; i ++)
            {
                int nt = (i + 1) % n;
                double tp = pt[i].x * pt[nt].y - pt[nt].x * pt[i].y;
                area += tp;
                ans.x += tp * (pt[i].x + pt[nt].x);
                ans.y += tp * (pt[i].y + pt[nt].y);
            }
            ans.x /= 3 * area;
            ans.y /= 3 * area;
            return ans;
        }
        //判断点是否在任意多边形内[射线法],O(n)
        bool ahas( Point &_Off)
        {
            int ret = 0;
            double infv = 1e20;//坐标系最大范围
            Line l = Line( _Off, Point( -infv ,_Off.y));
            for(int i = 0; i < n; i ++)
            {
                Line ln = Line( pt[i], pt[(i + 1) % n]);
                if(fabs(ln.s.y - ln.e.y) > eps)
                {
                    Point tp = (ln.s.y > ln.e.y)? ln.s: ln.e;
                    if( ( fabs( tp.y - _Off.y) < eps && tp.x < _Off.x + eps) || Line::isCrossSS( ln, l))
                        ret++;
                }
                else if( Line::isCrossSS( ln, l))
                    ret++;
            }
            return ret&1;
        }
    
        //判断任意点是否在凸包内,O(logn)
        bool bhas( Point & p)
        {
            if( n < 3)
                return false;
            if( xmult( pt[0], p, pt[1]) > eps)
                return false;
            if( xmult( pt[0], p, pt[n-1]) < -eps)
                return false;
            int l = 2,r = n-1;
            int line = -1;
            while( l <= r)
            {
                int mid = ( l + r) >> 1;
                if( xmult( pt[0], p, pt[mid]) >= 0)
                    line = mid,r = mid - 1;
                else l = mid + 1;
            }
            return xmult( pt[line-1], p, pt[line]) <= eps;
        }
    
    
    
        //凸多边形被直线分割
        Polygon split( Line &_Off)
        {
            //注意确保多边形能被分割
            Polygon ret;
            Point spt[2];
            double tp = 0.0, np;
            bool flag = true;
            int i, pn = 0, spn = 0;
            for(i = 0; i < n; i ++)
            {
                if(flag)
                    pt[pn ++] = pt[i];
                else
                    ret.pt[ret.n ++] = pt[i];
                np = xmult( _Off.s, _Off.e, pt[(i + 1) % n]);
                if(tp * np < -eps)
                {
                    flag = !flag;
                    Line::crossLPt( _Off, Line(pt[i], pt[(i + 1) % n]), spt[spn++]);
                }
                tp = (fabs(np) > eps)?np: tp;
            }
            ret.pt[ret.n ++] = spt[0];
            ret.pt[ret.n ++] = spt[1];
            n = pn;
            return ret;
        }
    
    
        /** 卷包裹法求点集凸包,_p为输入点集,_n为点的数量 **/
        void ConvexClosure( Point _p[], int _n)
        {
            sort( _p, _p + _n);
            n = 0;
            for(int i = 0; i < _n; i++)
            {
                while( n > 1 && sgn( xmult( pt[n-2], pt[n-1], _p[i]), 0) <= 0)
                    n--;
                pt[n++] = _p[i];
            }
            int _key = n;
            for(int i = _n - 2; i >= 0; i--)
            {
                while( n > _key && sgn( xmult( pt[n-2], pt[n-1], _p[i]), 0) <= 0)
                    n--;
                pt[n++] = _p[i];
            }
            if(n>1)   n--;//除去重复的点,该点已是凸包凸包起点
        }
        /****** 寻找凸包的graham 扫描法********************/
        /****** _p为输入的点集,_n为点的数量****************/
    
        void graham( Point _p[], int _n)
        {
            int cur=0;
            for(int i = 1; i < _n; i++)
                if( sgn( _p[cur].y, _p[i].y) > 0 || ( sgn( _p[cur].y, _p[i].y) == 0 && sgn( _p[cur].x, _p[i].x) > 0) )
                    cur = i;
            swap( _p[cur], _p[0]);
            n = 0, gsort = pt[n++] = _p[0];
            if( _n <= 1)   return;
            sort( _p + 1, _p+_n ,gcmp);
            pt[n++] = _p[1];
            for(int i = 2; i < _n; i++)
            {
                while(n>1 && sgn( xmult( pt[n-2], pt[n-1], _p[i]), 0) <= 0)// 当凸包退化成直线时需特别注意n
                    n--;
                pt[n++] = _p[i];
            }
        }
        //凸包旋转卡壳(注意点必须顺时针或逆时针排列)
        //返回值凸包直径的平方(最远两点距离的平方)
        pair<Point,Point> rotating_calipers()
        {
            int i = 1 % n;
            double ret = 0.0;
            pt[n] = pt[0];
            pair<Point,Point>ans=make_pair(pt[0],pt[0]);
            for(int j = 0; j < n; j ++)
            {
                while( fabs( xmult( pt[i+1], pt[j], pt[j + 1])) > fabs( xmult( pt[i], pt[j], pt[j + 1])) + eps)
                    i = (i + 1) % n;
                //pt[i]和pt[j],pt[i + 1]和pt[j + 1]可能是对踵点
                if(ret < getdis2(pt[i],pt[j]))  ret = getdis2(pt[i],pt[j]), ans = make_pair(pt[i],pt[j]);
                if(ret < getdis2(pt[i+1],pt[j+1]))  ret = getdis(pt[i+1],pt[j+1]), ans = make_pair(pt[i+1],pt[j+1]);
            }
            return ans;
        }
    
        //凸包旋转卡壳(注意点必须逆时针排列)
        //返回值两凸包的最短距离
        double rotating_calipers( Polygon &_Off)
        {
            int i = 0;
            double ret = 1e10;//inf
            pt[n] = pt[0];
            _Off.pt[_Off.n] = _Off.pt[0];
            //注意凸包必须逆时针排列且pt[0]是左下角点的位置
            while( _Off.pt[i + 1].y > _Off.pt[i].y)
                i = (i + 1) % _Off.n;
            for(int j = 0; j < n; j ++)
            {
                double tp;
                //逆时针时为 >,顺时针则相反
                while((tp = xmult(_Off.pt[i + 1],pt[j], pt[j + 1]) - xmult(_Off.pt[i], pt[j], pt[j + 1])) > eps)
                    i = (i + 1) % _Off.n;
                //(pt[i],pt[i+1])和(_Off.pt[j],_Off.pt[j + 1])可能是最近线段
                ret = min(ret, Line(pt[j], pt[j + 1]).dis(_Off.pt[i], true));
                ret = min(ret, Line(_Off.pt[i], _Off.pt[i + 1]).dis(pt[j + 1], true));
                if(tp > -eps)//如果不考虑TLE问题最好不要加这个判断
                {
                    ret = min(ret, Line(pt[j], pt[j + 1]).dis(_Off.pt[i + 1], true));
                    ret = min(ret, Line(_Off.pt[i], _Off.pt[i + 1]).dis(pt[j], true));
                }
            }
            return ret;
        }
    
        //-----------半平面交-------------
        //复杂度:O(nlog2(n))
        //获取半平面交的多边形(多边形的核)
        //参数:向量集合[l],向量数量[ln];(半平面方向在向量左边)
        //函数运行后如果n[即返回多边形的点数量]为0则不存在半平面交的多边形(不存在区域或区域面积无穷大)
        int judege( Line &_lx, Line &_ly, Line &_lz)
        {
            Point tmp;
            Line::crossLPt(_lx,_ly,tmp);
            return sgn(xmult(_lz.s,tmp,_lz.e),0);
        }
        int halfPanelCross(Line L[], int ln)
        {
            int i, tn, bot, top;
            for(int i = 0; i < ln; i++)
                L[i].get_angle();
            sort(L, L + ln);
            //平面在向量左边的筛选
            for(i = tn = 1; i < ln; i ++)
                if(fabs(L[i].angle - L[i - 1].angle) > eps)
                    L[tn ++] = L[i];
            ln = tn, n = 0, bot = 0, top = 1;
            dq[0] = L[0], dq[1] = L[1];
            for(i = 2; i < ln; i ++)
            {
                while(bot < top &&  judege(dq[top],dq[top-1],L[i]) > 0)
                    top --;
                while(bot < top &&  judege(dq[bot],dq[bot+1],L[i]) > 0)
                    bot ++;
                dq[++ top] = L[i];
            }
            while(bot < top && judege(dq[top],dq[top-1],dq[bot]) > 0)
                top --;
            while(bot < top && judege(dq[bot],dq[bot+1],dq[top]) > 0)
                bot ++;
            //若半平面交退化为点或线
            //        if(top <= bot + 1)
            //            return 0;
            dq[++top] = dq[bot];
            for(i = bot; i < top; i ++)
                Line::crossLPt(dq[i],dq[i + 1],pt[n++]);
            return n;
        }
    };
    
    
    class Circle
    {
    public:
        Point c;//圆心
        double r;//半径
        double db, de;//圆弧度数起点, 圆弧度数终点(逆时针0-360)
    
        //-------圆---------
    
        //判断圆在多边形内
        bool inside( Polygon &_Off)
        {
            if(_Off.ahas(c) == false)
                return false;
            for(int i = 0; i < _Off.n; i ++)
            {
                Line l = Line(_Off.pt[i], _Off.pt[(i + 1) % _Off.n]);
                if(l.dis(c, true) < r - eps)
                    return false;
            }
            return true;
        }
    
        //判断多边形在圆内(线段和折线类似)
        bool has( Polygon &_Off)
        {
            for(int i = 0; i < _Off.n; i ++)
                if( getdis2(_Off.pt[i],c) > r * r - eps)
                    return false;
            return true;
        }
    
        //-------圆弧-------
        //圆被其他圆截得的圆弧,参数:圆[_Off]
        Circle operator-(Circle &_Off) const
        {
            //注意圆必须相交,圆心不能重合
            double d2 = getdis2(c,_Off.c);
            double d = getdis(c,_Off.c);
            double ans = acos((d2 + r * r - _Off.r * _Off.r) / (2 * d * r));
            Point py = _Off.c - c;
            double oans = atan2(py.y, py.x);
            Circle res;
            res.c = c;
            res.r = r;
            res.db = oans + ans;
            res.de = oans - ans + 2 * PI;
            return res;
        }
        //圆被其他圆截得的圆弧,参数:圆[_Off]
        Circle operator+(Circle &_Off) const
        {
            //注意圆必须相交,圆心不能重合
            double d2 = getdis2(c,_Off.c);
            double d = getdis(c,_Off.c);
            double ans = acos((d2 + r * r - _Off.r * _Off.r) / (2 * d * r));
            Point py = _Off.c - c;
            double oans = atan2(py.y, py.x);
            Circle res;
            res.c = c;
            res.r = r;
            res.db = oans - ans;
            res.de = oans + ans;
            return res;
        }
    
        //过圆外一点的两条切线
        //参数:点[_Off](必须在圆外),返回:两条切线(切线的s点为_Off,e点为切点)
        pair<Line, Line>  tangent( Point &_Off)
        {
            double d = getdis(c,_Off);
            //计算角度偏移的方式
            double angp = acos(r / d), ango = atan2(_Off.y - c.y, _Off.x - c.x);
            Point pl = Point(c.x + r * cos(ango + angp), c.y + r * sin(ango + angp)),
                pr = Point(c.x + r * cos(ango - angp), c.y + r * sin(ango - angp));
            return make_pair(Line(_Off, pl), Line(_Off, pr));
        }
    
        //计算直线和圆的两个交点
        //参数:直线[_Off](两点式),返回两个交点,注意直线必须和圆有两个交点
        pair<Point, Point> cross(Line _Off)
        {
            _Off.pton();
            //到直线垂足的距离
            double td = fabs(_Off.a * c.x + _Off.b * c.y + _Off.c) / sqrt(_Off.a * _Off.a + _Off.b * _Off.b);
    
            //计算垂足坐标
            double xp = (_Off.b * _Off.b * c.x - _Off.a * _Off.b * c.y - _Off.a * _Off.c) / ( _Off.a * _Off.a + _Off.b * _Off.b);
            double yp = (- _Off.a * _Off.b * c.x + _Off.a * _Off.a * c.y - _Off.b * _Off.c) / (_Off.a * _Off.a + _Off.b * _Off.b);
    
            double ango = atan2(yp - c.y, xp - c.x);
            double angp = acos(td / r);
    
            return make_pair(Point(c.x + r * cos(ango + angp), c.y + r * sin(ango + angp)),
                Point(c.x + r * cos(ango - angp), c.y + r * sin(ango - angp)));
        }
    };
    
    class triangle
    {
    public:
        Point a, b, c;//顶点
        triangle(){}
        triangle(Point a, Point b, Point c): a(a), b(b), c(c){}
    
        //计算三角形面积
        double area()
        {
            return fabs( xmult(a, b, c)) / 2.0;
        }
    
        //计算三角形外心
        //返回:外接圆圆心
        Point circumcenter()
        {
            double pa = a.x * a.x + a.y * a.y;
            double pb = b.x * b.x + b.y * b.y;
            double pc = c.x * c.x + c.y * c.y;
            double ta = pa * ( b.y - c.y) - pb * ( a.y - c.y) + pc * ( a.y - b.y);
            double tb = -pa * ( b.x - c.x) + pb * ( a.x - c.x) - pc * ( a.x - b.x);
            double tc = a.x * ( b.y - c.y) - b.x * ( a.y - c.y) + c.x * ( a.y - b.y);
            return Point( ta / 2.0 / tc, tb / 2.0 / tc);
        }
    
        //计算三角形内心
        //返回:内接圆圆心
        Point incenter()
        {
            Line u, v;
            double m, n;
            u.s = a;
            m = atan2(b.y - a.y, b.x - a.x);
            n = atan2(c.y - a.y, c.x - a.x);
            u.e.x = u.s.x + cos((m + n) / 2);
            u.e.y = u.s.y + sin((m + n) / 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.x = v.s.x + cos((m + n) / 2);
            v.e.y = v.s.y + sin((m + n) / 2);
            Point ret;
            Line::crossLPt(u,v,ret);
            return ret;
        }
    
        //计算三角形垂心
        //返回:高的交点
        Point perpencenter()
        {
            Line u,v;
            u.s = c;
            u.e.x = u.s.x - a.y + b.y;
            u.e.y = u.s.y + a.x - b.x;
            v.s = b;
            v.e.x = v.s.x - a.y + c.y;
            v.e.y = v.s.y + a.x - c.x;
            Point ret;
            Line::crossLPt(u,v,ret);
            return ret;
        }
    
        //计算三角形重心
        //返回:重心
        //到三角形三顶点距离的平方和最小的点
        //三角形内到三边距离之积最大的点
        Point barycenter()
        {
            Line u,v;
            u.s.x = (a.x + b.x) / 2;
            u.s.y = (a.y + b.y) / 2;
            u.e = c;
            v.s.x = (a.x + c.x) / 2;
            v.s.y = (a.y + c.y) / 2;
            v.e = b;
            Point ret;
            Line::crossLPt(u,v,ret);
            return ret;
        }
    
        //计算三角形费马点
        //返回:到三角形三顶点距离之和最小的点
        Point fermentPoint()
        {
            Point u, v;
            double step = fabs(a.x) + fabs(a.y) + fabs(b.x) + fabs(b.y) + fabs(c.x) + fabs(c.y);
            int i, j, k;
            u.x = (a.x + b.x + c.x) / 3;
            u.y = (a.y + b.y + c.y) / 3;
            while (step > eps)
            {
                for (k = 0; k < 10; step /= 2, k ++)
                {
                    for (i = -1; i <= 1; i ++)
                    {
                        for (j =- 1; j <= 1; j ++)
                        {
                            v.x = u.x + step * i;
                            v.y = u.y + step * j;
                            if (getdis(u,a) + getdis(u,b) + getdis(u,c) > getdis(v,a) + getdis(v,b) + getdis(v,c))
                                u = v;
                        }
                    }
                }
            }
            return u;
        }
    };
    
    int main(){
        return 0;
    }
    
  • 相关阅读:
    [HNOI2002]营业额统计
    HDU 1374
    HDU 3345
    HDU 2089
    Graham扫描法
    Codeforces 1144D Deduction Queries 并查集
    Codeforces 916E Jamie and Tree 线段树
    Codeforces 1167F Scalar Queries 树状数组
    Codeforces 1167E Range Deleting
    Codeforces 749E Inversions After Shuffle 树状数组 + 数学期望
  • 原文地址:https://www.cnblogs.com/ucprer/p/12867541.html
Copyright © 2011-2022 走看看