Geometry/include.hpp

• View this file on GitHub
• Last update: 2023-04-05 23:10:22+09:00
• Include: #include "Geometry/include.hpp"

Depends on

• Geometry/Angle.hpp
• Geometry/CircumscribedCircle.hpp
• Geometry/ClosestPair.hpp
• Geometry/Contain.hpp
• Geometry/Convex.hpp
• Geometry/CounterClockWise.hpp
• Geometry/Cross.hpp
• Geometry/CrossPoint.hpp
• Geometry/Distance.hpp
• Geometry/Dot.hpp
• Geometry/InscribedCircle.hpp
• Geometry/Intersect.hpp
• Geometry/MinimumBoundingCircle.hpp
• Geometry/Projection.hpp
• Geometry/Rotate.hpp
• Geometry/Tangent.hpp
• Geometry/template.hpp

Code

#include "./template.hpp"
#include "./Rotate.hpp"
#include "./Dot.hpp"
#include "./Cross.hpp"
#include "./CounterClockWise.hpp"
#include "./Projection.hpp"
#include "./Intersect.hpp"
#include "./Distance.hpp"
#include "./CrossPoint.hpp"
#include "./Angle.hpp"
#include "./InscribedCircle.hpp"
#include "./CircumscribedCircle.hpp"
#include "./Tangent.hpp"
#include "./Contain.hpp"
#include "./MinimumBoundingCircle.hpp"
#include "./ClosestPair.hpp"
#include "./Convex.hpp"

#line 1 "Geometry/template.hpp"
// Real
using Real=double;
const Real EPS=1e-6;
const Real pi=acosl(-1);

// Point
using Point=complex<Real>;
istream &operator>>(istream &is,Point &p){
    Real a,b;
    is>>a>>b;
    p=Point(a,b);
    return is;
}
ostream &operator<<(ostream &os,Point &p){
    return os<<fixed<<setprecision(12)<<p.real()<<' '<<p.imag();
}
inline bool eq(Real a,Real b){
    return fabs(a-b)<EPS;
}
Point operator*(const Point &p,const Real &d){
    return Point(real(p)*d,imag(p)*d);
}

// Line
struct Line{
    Point p1,p2;
    Line()=default;
    Line(Point p1,Point p2):p1(p1),p2(p2){}
    //Ax + By = C
    Line(Real A,Real B,Real C){
        if(eq(A,0))     p1=Point(0,C/B),p2=Point(1,C/B);
        else if(eq(B,0))p1=Point(C/A,0),p2=Point(C/A,1);
        else            p1=Point(0,C/B),p2=Point(C/A,0);
    }
};

// Segment
struct Segment:Line{
    Segment()=default;
    Segment(Point p1,Point p2):Line(p1,p2){}
};
struct Circle{
    Point center;
    Real r;
    Circle()=default;
    Circle(Point center,Real r):center(center),r(r){}
};

// Polygon
using Polygon=vector<Point>;
#line 1 "Geometry/Rotate.hpp"
Point rotate(Real theta,Point p){
    return Point(cos(theta)*real(p)-sin(theta)*imag(p),sin(theta)*real(p)+cos(theta)*imag(p));
}
#line 1 "Geometry/Dot.hpp"
// Dot
Real dot(Point a,Point b) {
    return real(a)*real(b)+imag(a)*imag(b);
}
#line 1 "Geometry/Cross.hpp"
// Cross
Real cross(Point a,Point b){
    return real(a)*imag(b)-imag(a)*real(b);
}
#line 1 "Geometry/CounterClockWise.hpp"
// ccw (counter clockwise) (Requires: cross, dot)
//https://onlinejudge.u-aizu.ac.jp/courses/library/4/CGL/all/CGL_1_C
int ccw(Point a,Point b,Point c){
   b-=a;c-=a;
   if(cross(b,c)>EPS)       return  1;//COUNTER CLOCKWISE
   else if(cross(b,c)<-EPS) return -1;//CLOCKWISE
   else if(dot(b,c)<0)      return  2;//c--a--b ONLINE BACK
   else if(norm(b)<norm(c)) return -2;//a--b--c ONLINE FRONT
   else                     return  0;//a--c--b ON SEGMENT
}
#line 1 "Geometry/Projection.hpp"
// Projection (Requires: dot)
Point projection(Line l,Point p){
    // ベクトルl乗に点pからおろした垂線の足
    Real k=dot(l.p1-l.p2,p-l.p1)/norm(l.p1-l.p2);
    return l.p1+(l.p1-l.p2)*k;
}
Point projection(Segment l,Point p){
    Real k=dot(l.p1-l.p2,p-l.p1)/norm(l.p1-l.p2);
    return l.p1+(l.p1-l.p2)*k;
}
#line 1 "Geometry/Intersect.hpp"
// Intersect (Requires : ccw, Dots, Cross, Projection)
bool intersect(Line l,Point p){
    return abs(ccw(l.p1,l.p2,p))!=1;
}
//直線の交差判定，外積
bool intersect(Line l1,Line l2){
    return abs(cross(l1.p2-l1.p1,l2.p2-l2.p1))>EPS or abs(cross(l1.p2-l1.p1,l2.p2-l1.p1))<EPS;
}
//線分に点が乗るかの判定，ccw
bool intersect(Segment s,Point p){
    return ccw(s.p1,s.p2,p)==0;
}
//直線と線分の交差判定
bool intersect(Line l,Segment s){
    return cross(l.p2-l.p1,s.p1-l.p1)*cross(l.p2-l.p1,s.p2-l.p1)<EPS;
}
//円と直線の交差判定
bool intersect(Circle c,Line l){
    return abs(c.center-projection(l,c.center))<=c.r+EPS;
}
//円上かどうか，内部かどうかではない
bool intersect(Circle c,Point p){
    return abs(abs(p-c.center)-c.r)<EPS;
}
//線分と線分の交差判定
bool intersect(Segment s,Segment t){
    return ccw(s.p1,s.p2,t.p1)*ccw(s.p1,s.p2,t.p2) <=0 and ccw(t.p1,t.p2,s.p1)*ccw(t.p1,t.p2,s.p2)<=0;
}
//線分と円の交差判定，交点の個数を返す
int intersect(Circle c,Segment l){
    Point h=projection(l,c.center);
    //直線まるっと円の外側
    if(norm(h-c.center)-c.r*c.r>EPS)    return 0;
    Real d1=abs(c.center-l.p1),d2=abs(c.center-l.p2);
    //線分が円内
    if(d1<c.r+EPS and d2<c.r+EPS) return 0;
    if((d1<c.r-EPS and d2>c.r+EPS) or (d2<c.r-EPS and d1>c.r+EPS)) return 1;
    //円の外部にまるまるはみ出ていないか
    if(dot(l.p1-h,l.p2-h)<0) return 2;
    return 0;
}
//円と円の位置関係，共通接線の個数を返す
int intersect(Circle c1,Circle c2){
    if(c1.r<c2.r) swap(c1,c2);
    Real d=abs(c1.center-c2.center);
    //2円が離れている
    if(c1.r+c2.r<d)     return 4;
    //2円が外接する
    if(eq(c1.r+c2.r,d)) return 3;
    //2円が交わる
    if(c1.r-c2.r<d)     return 2;
    //円が内接する
    if(eq(c1.r-c2.r,d)) return 1;
    //内包
                        return 0;
}
#line 1 "Geometry/Distance.hpp"
// Distance (Requires: Projection, Intersect)
Real dis(Point a,Point b){
    return abs(a-b);
}
Real dis(Line l,Point p){
    return abs(p-projection(l,p));
}
Real dis(Segment s,Point p){
    Point r=projection(s,p);
    if(intersect(s,r)) return abs(r-p);
    return min(abs(s.p1-p),abs(s.p2-p));
}
Real dis(Segment a,Segment b){
    if(intersect(a,b)) return 0;
    return min({dis(a,b.p1),dis(a,b.p2),dis(b,a.p1),dis(b,a.p2)});
}
Real dis(Polygon a,Polygon b){
    Real ret=-10;
    int n=(int)a.size(),m=(int)b.size();
    for(int i=0;i<n;i++)for(int j=0;j<m;j++){
        Real d=dis(Segment(a[i],a[(i+1)%n]),Segment(b[j],b[(j+1)%m]));
        if(ret<0) ret=d;
        else      ret=min(ret,d);
    }
    return ret;
}
Real dis(Polygon poly,Point p){
    Real ret=-10;
    int n=(int)poly.size();
    for(int i=0;i<n;i++){
        Real d=dis(Segment(poly[i],poly[(i+1)%n]),p);
        if(ret<0) ret=d;
        else      ret=min(ret,d);
    }
    return ret;
}
#line 1 "Geometry/CrossPoint.hpp"
//intersectをチェックすること
//v
Point crosspoint(Line l,Line m){
   Real A=cross(m.p2-m.p1,m.p1-l.p1);
   Real B=cross(m.p2-m.p1,l.p2-l.p1);
   if(eq(A,0) and eq(B,0)) return l.p1;
   if(eq(B,0))             throw "NAI";
   return l.p1+A/B*(l.p2-l.p1);   
}
Point crosspoint(Segment l,Segment m){
   return crosspoint(Line(l),Line(m));
}
vector<Point> crosspoint(Circle c,Line l){
   vector<Point> ret;
   Point h=projection(l,c.center);
   Real d=sqrt(c.r*c.r-norm(h-c.center));
   Point e=(l.p2-l.p1)*(1/abs(l.p2-l.p1));
   if(c.r*c.r+EPS<norm(h-c.center)) return ret;
   if(eq(dis(l,c.center),c.r)){
       ret.push_back(h);
       return ret;
   }
   ret.push_back(h+e*d);ret.push_back(h-e*d);
   return ret;
}
//要verify，
vector<Point> crosspoint(Circle c,Segment s){
   Line l=Line(s.p1,s.p2);
   int ko=intersect(c,s);
   if(ko==2) return crosspoint(c,l);
   vector<Point> ret;
   if(ko==0) return ret;
   ret=crosspoint(c,l);
   if(ret.size()==1) return ret;
   vector<Point> rret;
   //交点で挟める方を返す
   if(dot(s.p1-ret[0],s.p2-ret[0])<0)  rret.push_back(ret[0]);
   else                                rret.push_back(ret[1]);
   return rret;
}
//v
vector<Point> crosspoint(Circle c1,Circle c2){
   vector<Point> ret;
   int isec=intersect(c1,c2);
   if(isec==0 or isec==4) return ret;
   Real d=abs(c1.center-c2.center);
   Real a=acos((c1.r*c1.r+d*d-c2.r*c2.r)/(2*c1.r*d));
   Real t=atan2(c2.center.imag()-c1.center.imag(),c2.center.real()-c1.center.real());
   ret.push_back(c1.center+Point(cos(t+a)*c1.r,sin(t+a)*c1.r));
   ret.push_back(c1.center+Point(cos(t-a)*c1.r,sin(t-a)*c1.r));
   return ret;
}
#line 1 "Geometry/Angle.hpp"
// angle of a-b-c
Real get_smaller_angle(Point a,Point b,Point c){
    Point v=a-b,w=c-b;
    auto A=atan2(imag(v),real(v));
    auto B=atan2(imag(w),real(w));
    if(A>B) swap(A,B);
    Real res=B-A;
    return min(res,pi*2.0-res);
}
#line 1 "Geometry/InscribedCircle.hpp"
// 内接円
Circle inscribed_circle(Point a,Point b,Point c){
    Real A,B;
    {
        Point t=c-a;
        t*=conj(b-a);
        t/=norm(b-a);
        A=atan2(imag(t),real(t));
    }
    {
        Point t=a-b;
        t*=conj(c-b);
        t/=norm(c-b);
        B=atan2(imag(t),real(t));
    }
    Line Amid=Line(a,a+rotate(A*0.5,b-a)),Bmid=Line(b,b+rotate(B*0.5,c-b));
    auto center=crosspoint(Amid,Bmid);
    auto h=projection(Line(a,b),center);
    return Circle(center,dis(h,center));
}
#line 1 "Geometry/CircumscribedCircle.hpp"
// 外接円
Circle circumscribed_circle(Point a,Point b,Point c){
    Line orth_ab((a+b)*0.5,(a+b)*0.5+Point(-imag(b-a),real(b-a)));
    Line orth_bc((b+c)*0.5,(b+c)*0.5+Point(-imag(c-b),real(c-b)));
    Point center=crosspoint(orth_ab,orth_bc);
    Real r=dis(a,center);
    return Circle(center,r);
}
#line 1 "Geometry/Tangent.hpp"
//v
//点pから引いた円cの接線の接点を返す
vector<Point> tangent(Circle c,Point p){
    return crosspoint(c,Circle(p,sqrt(norm(c.center-p)-c.r*c.r)));
}
//v
//二円の共通接線，Lineの2点は接点を表す
vector<Line> tangent(Circle c1,Circle c2){
    vector<Line> ret;
    if(c1.r<c2.r) swap(c1,c2);
    Real g=norm(c1.center-c2.center);
    if(eq(g,0)) return ret;
    Point u=(c2.center-c1.center)/sqrt(g);
    Point v=rotate(pi*0.5,u);
    for(int s:{-1,1}){
        Real h=(c1.r+s*c2.r)/sqrt(g);
        if(eq(1-h*h,0)){
            ret.push_back(Line(c1.center+u*c1.r,c1.center+(u+v)*c1.r));
        }
        else if(1-h*h>0){
            Point uu=u*h,vv=v*sqrt(1-h*h);
            ret.push_back(Line(c1.center+(uu+vv)*c1.r,c2.center-(uu+vv)*c2.r*s));
            ret.push_back(Line(c1.center+(uu-vv)*c1.r,c2.center-(uu-vv)*c2.r*s));
        }
    }
    return ret;
}
#line 1 "Geometry/Contain.hpp"
// out 0, on 1, in 2
int contains(Polygon poly,Point p){
    int res=0;
    int n=(int)poly.size();
    for(int i=0;i<n;i++){
        Point a=poly[i]-p,b=poly[(i+1)%n]-p;
        if(imag(a)>imag(b)) swap(a,b);
        if(imag(a)<=0 and 0<imag(b) and cross(a,b)<0) res^=1;
        if(eq(cross(a,b),0) and (dot(a,b)<0 or eq(dot(a,b),0))) return 1;
    }
    if(res) res=2;
    return res;
}
#line 1 "Geometry/MinimumBoundingCircle.hpp"
//最小包含円を返す　計算量は期待値O(n)
Circle MinimumBoundingCircle(vector<Point> v){
    int n=v.size();
    //ランダムシャッフル．いぢわるされたくないもんだ
    mt19937 mt(time(0));
    shuffle(v.begin(),v.end(),mt);
    Circle ret(0,0);

    auto make_circle2=[&](Point a,Point b){
        return Circle((a+b)*0.5,dis(a,b)/2);
    };

    auto make_circle3=[&](Point A,Point B,Point C){
        Point cent=circumscribed_circle(A,B,C).center;
        return Circle(cent,dis(cent,A));
    };

    auto isIn=[&](Point a){
        return dis(ret.center,a)<ret.r+EPS;
    };

    ret=make_circle2(v[0],v[1]);
    for(int i=2;i<n;i++){
        //v[i]が円に入っていないなら
        if(!isIn(v[i])){
            //円内にないなら点v[i]は必ず円周上に来る
            ret=make_circle2(v[0],v[i]);
            for(int j=1;j<i;j++){
                if(!isIn(v[j])){
                    //この時iとjが円周上を考える
                    ret=make_circle2(v[i],v[j]);
                    //最後の1点の決定
                    for(int k=0;k<j;k++)if(!isIn(v[k])) ret=make_circle3(v[i],v[j],v[k]);
                }
            }
        }
    }
    return ret;
}
#line 1 "Geometry/ClosestPair.hpp"
// 最近点対
// O(NlogN)
Real closest_pair(vector<Point> ps){
    sort(ALL(ps),[&](Point a,Point b){
        return real(a)<real(b);
    });
    function<Real(int,int)> rec=[&](int l,int r){
        if(r-l<=1) return (Real)1e18;
        int m=(l+r)/2;
        Real x=real(ps[m]);
        Real ret=min(rec(l,m),rec(m,r));
        inplace_merge(begin(ps)+l,begin(ps)+m,begin(ps)+r,[&](Point a,Point b){
            return imag(a)<imag(b);
        });
        // 分割を跨いで最小距離があるか調べる
        vector<Point> b;
        for(int i=l;i<r;i++){
            if(abs(real(ps[i])-x)>=ret) continue;
            for(int j=(int)b.size()-1;j>=0;j--){
                if(abs(imag(ps[i]-b[j]))>=ret) break;
                ret=min(ret,abs(ps[i]-b[j]));
            }
            b.push_back(ps[i]);
        }
        return ret;
    };
    return rec(0,(int)ps.size());
}
#line 1 "Geometry/Convex.hpp"
// 凸多角形系統
// 凸多角形の頂点は反時計周りに訪れる順序
// v
// 頂点は反時計周りに訪れる順序，時計回りとなるような3点があるとfalse
bool is_convex(const vector<Point> &ps){
    int n=(int)ps.size();
    for(int i=0;i<n;i++)if(ccw(ps[(i+n-1)%n],ps[i],ps[(i+1)%n])==-1)return false;
    return true;
}

// 凸包，あんまりよくわかってない．直線状に頂点をのせない場合(↑)，のせる場合(↓)
vector<Point> convex_hull(vector<Point> p){
    int n=(int)p.size(),k=0;
    if(n<=2)return p;
    sort(begin(p),end(p),[](Point a,Point b){
        return real(a)!=real(b)?real(a)<real(b):imag(a)<imag(b);
    });
    vector<Point>ch(2*n);
    for(int i=0;i<n;ch[k++]=p[i++]){
        // while(k>=2 and cross(ch[k-1]-ch[k-2],p[i]-ch[k-1])<EPS)k--;
        while(k>=2 and cross(ch[k-1]-ch[k-2],p[i]-ch[k-1])<0)k--;
    }
    for(int i=n-2,t=k+1;i>=0;ch[k++]=p[i--]){
        // while(k>=t and cross(ch[k-1]-ch[k-2],p[i]-ch[k-1])<EPS)k--;
        while(k>=t and cross(ch[k-1]-ch[k-2],p[i]-ch[k-1])<0)k--;
    }
    ch.resize(k-1);
    return ch;
}
#line 18 "Geometry/include.hpp"

Back to top page