计算几何模版

const double eps=1e-8;

int dcmp(double x){
    if (fabs(x)<eps) {
        return 0;
    }
    else return x<0?-1:1;
}

struct Point{
    double x,y;
    Point(double x=0,double y=0):x(x),y(y){}
    int in(){
        return scanf("%lf%lf",&x,&y);
    }
};

typedef Point Vector;

struct Circle{
    Point c;
    double r;
    Circle(){}
    Circle(Point c,double r):c(c),r(r){}
    Point point(double a){
        return Point((c.x+cos(a)*r),c.y+sin(a)*r);
    }
    
    void in(){
        c.in();
        scanf("%lf",&r);
    }
    
};


inline Vector operator+ (Vector A, Vector B) { return Vector(A.x + B.x, A.y + B.y); }
inline Vector operator- (Point A, Point B) { return Vector(A.x - B.x, A.y - B.y); }
inline Vector operator* (Vector A, double p) { return Vector(A.x * p, A.y * p); }
inline Vector operator/ (Vector A, double p) { return Vector(A.x / p, A.y / p); }
inline bool operator < (Point a, Point b) { return a.x < b.x || (a.x == b.x && a.y < b.y); }
inline bool operator == (Point a, Point b) { return dcmp(a.x - b.x) == 0 && dcmp(a.y - b.y) == 0; }

inline long double Dot(Vector A, Vector B){ return A.x * B.x + A.y * B.y;}

inline long double Length(Vector A) { return sqrt(Dot(A, A)); }

inline Vector Unit(Vector x) { return x / Length(x); }  //单位向量

inline double Angle(Vector A, Vector B) { return acos(Dot(A, B) / Length(A) / Length(B)); }

inline double angle(Vector v) { return atan2(v.y, v.x); }

inline long double Cross(Vector A, Vector B) { return A.x * B.y - A.y * B.x; }

inline bool OnSegment(Point p, Point a1, Point a2)//包括端点
{
    return dcmp(Cross(a1-p, a2-p)) == 0 && dcmp(Dot(a1-p, a2-p)) <= 0;
}

inline double DistanceToLine(Point P, Point A, Point B)
{
    Vector v1 = B - A, v2 = P - A;
    return fabs(Cross(v1, v2)) / Length(v1); // 如果不取绝对值，得到的是有向距离
}

Vector Rotate(Vector A,double rad){   //旋转
    return Vector(A.x*cos(rad)-A.y*sin(rad),A.x*sin(rad)+A.y*cos(rad));
}


struct Line{
    Point p;
    Vector v;
    double ang;
    Line(){}
    Line(Point p,Vector v):p(p),v(v){
        ang=atan2(v.y, v.x);
    }
    inline Point point(double t)
    {
        return p + v * t;
    }
    bool operator <(Line a) const{
        return ang<a.ang;
    }
};

//三边构成三角形的判定
bool check_length(double a, double b, double c) {  
    return dcmp(a+b-c) > 0 && dcmp(fabs(a-b)-c) < 0;  
}  
bool isTriangle(double a, double b, double c) {  
    return check_length(a, b, c) && check_length(a, c, b) && check_length(b, c, a);  
}

//平行四边形的判定（保证四边形顶点按顺序给出）  
bool isParallelogram(Point  *p) {  
    if (dcmp(Length(p[0]-p[1]) - Length(p[2]-p[3])) || dcmp(Length(p[0]-p[3]) - Length(p[2]-p[1]))) return false;       
    Line a = Line(p[0], p[1]-p[0]);  
    Line b = Line(p[1], p[2]-p[1]);  
    Line c = Line(p[3], p[2]-p[3]);  
    Line d = Line(p[0], p[3]-p[0]);  
    return dcmp(a.ang - c.ang) == 0 && dcmp(b.ang - d.ang) == 0;  
}

//梯形的判定  顺序输入
bool isTrapezium(Point *p) {  
    Line a = Line(p[0], p[1]-p[0]);  
    Line b = Line(p[1], p[2]-p[1]);  
    Line c = Line(p[3], p[2]-p[3]);  
    Line d = Line(p[0], p[3]-p[0]);  
    return (dcmp(a.ang - c.ang) == 0 && dcmp(b.ang - d.ang)) || (dcmp(a.ang - c.ang) && dcmp(b.ang - d.ang) == 0);  
} 
//菱形的判定  顺序输入
bool isRhombus(Point *p) {  
    if (!isParallelogram(p)) return false;  
    return dcmp(Length(p[1]-p[0]) - Length(p[2]-p[1])) == 0;  
} 
//矩形的判定  
bool isRectangle(Point *p) {  
    if (!isParallelogram(p)) return false;  
    return dcmp(Length(p[2]-p[0]) - Length(p[3]-p[1])) == 0;  
}  

//正方形的判定  
bool isSquare(Point *p) {  
    return isRectangle(p) && isRhombus(p);  
}  
   
//三点共线的判定  
bool isCollinear(Point A, Point B, Point C) {  
    return dcmp(Cross(B-A, C-B)) == 0;  
} 

//点在多边形内的判定(多边形顶点需按逆时针排列)  
bool isInPolygon(Point p,Point *poly,int n) {  
    for(int i = 0; i < n; i++) {  
        //若允许点在多边形边上，可关闭下行注释  
        // if (isOnSegment(p, poly[(i+1)%n], poly[i])) return true;  
        if (Cross(poly[(i+1)%n]-poly[i], p-poly[i]) < 0) return false;  
    }  
    return true;  
}

//过定点作圆的切线  
int getTangents(Point P, Circle C, Line *L,int &n) {  
    Vector u = C.c - P;  
    double dis = Length(u);  
    if (dcmp(dis - C.r) < 0) return 0;  
    if (dcmp(dis - C.r) == 0) {  
        L[n++]={P, Rotate(u, pi / 2.0)};  
        return 1;  
    }  
    double ang = asin(C.r / dis);  
    L[n++]=Line(P, Rotate(u, ang));  
    L[n++]=Line(P, Rotate(u, -ang));  
    return 2;  
} 

//两圆位置关系判定  
int GetCircleLocationRelation(Circle A, Circle B){ 
    double d = Length(A.c-B.c);  
    double sum = A.r + B.r;  
    double sub = fabs(A.r - B.r);  
    if (dcmp(d) == 0) return dcmp(sub) != 0;  
    if (dcmp(d - sum) > 0) return 2;//XiangLi;  
    if (dcmp(d - sum) == 0) return 1;//WaiQie;  
    if (dcmp(d - sub) > 0 && dcmp(d - sum) < 0) return 0;//INTERSECTING;  
    if (dcmp(d - sub) == 0) return -1;//NeiQie;  
    if (dcmp(d - sub) < 0) return -2;//NeiHan;  
}

//两圆相交的面积
double GetCircleIntersectionArea(Circle A, Circle B) {
    int rel = GetCircleLocationRelation(A, B);
    if (rel < 0/*INTERSECTING*/) return min(A.r*A.r*pi, B.r*B.r*pi);
    if (rel > 0/*INTERSECTING*/) return 0;
    double dis = Length(A.c - B.c);
    double ang1 = acos((A.r*A.r + dis*dis - B.r*B.r) / (2.0*A.r*dis));
    double ang2 = acos((B.r*B.r + dis*dis - A.r*A.r) / (2.0*B.r*dis));
    return ang1*A.r*A.r + ang2*B.r*B.r - A.r*dis*sin(ang1); 
}

inline Point GetLineProjection(Point P, Point A, Point B)//点在直线上的投影
{
    Vector v = B - A;
    return A + v * (Dot(v, P - A) / Dot(v, v));
}

inline Point getmiorr(Point P,Line a){
    Point pp=GetLineProjection(P,a.p,a.p+a.v);
    return pp*2-P;
}


inline bool Seg_inter_Line(Line l1,Point a,Point b){//判断直线l1 是否相交线段ab
    return dcmp(Cross(l1.p-a,l1.p-a+l1.v))*dcmp(Cross(l1.p-b, l1.p+l1.v-b))<=0;
}

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

int getLineCircleIntersection(Line L, Circle C, double& t1, double& t2, vector<Point>& sol)
{
    double a = L.v.x, b = L.p.x - C.c.x, c = L.v.y, d = L.p.y - C.c.y;
    double e = a*a + c*c, f = 2*(a*b + c*d), g = b*b + d*d - C.r*C.r;
    double delta = f*f - 4*e*g; // 判别式
    if(dcmp(delta) < 0) return 0; // 相离
    if(dcmp(delta) == 0)
    {
        // 相切
        t1 = t2 = -f / (2 * e);
        sol.push_back(L.point(t1));
        return 1;
    }
    // 相交
    t1 = (-f - sqrt(delta)) / (2 * e);
    sol.push_back(L.point(t1));
    t2 = (-f + sqrt(delta)) / (2 * e);
    sol.push_back(L.point(t2));
    return 2;
}

long 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);
    else if(dcmp(Dot(v1, v3))>0) return Length(v3);
    else return fabs(Cross(v1, v2))/Length(v1);
}

bool SegmentProperIntersection(Point a1,Point a2,Point b1,Point b2){ //线段相交
    double c1=Cross(a2-a1,b1-a1),c2=Cross(a2-a1, b2-a1),c3=Cross(b2-b1, a1-b1),c4=Cross(b2-b1, a2-b1);
    return dcmp(c1)*dcmp(c2)<=0&&dcmp(c3)*dcmp(c4)<=0;
}

//二维凸包
int n,res[maxn],top;

bool cmp(Point a,Point b){
    if(a.x==b.x)return a.y<b.y;
    return a.x<b.x;
}

int Cross(Point a,Point b){
    return a.x*b.y-a.y*b.x;
}

void ConvexHull(){
    int len;
    top=0;
    sort(p,p+n,cmp);
    for(int i=0;i<n;i++){
        while(top>1&&Cross(p[res[top-1]]-p[res[top-2]],p[i]-p[res[top-2]])<=0)
            top--;
        res[top++]=i;
    }
    len=top;
    for(int i=n-2;i>=0;i--){
        while(top>len&&Cross(p[res[top-1]]-p[res[top-2]],p[i]-p[res[top-2]])<=0)top--;
        res[top++]=i;
    }
    if(n>1)--top;
    //  cout<<top<<endl;
}

// 半平面交
bool OnLeft(Line L,Point p){
    return Cross(L.v, p-L.p)>0;
}

Point GetIntersection(Line a,Line b){
    Vector u=a.p-b.p;
    double t=Cross(b.v, u)/Cross(a.v, b.v);
    return a.p+a.v*t;
}

double ConvexPolygonArea(Point *p,int n){
    double area=0;
    for (int i=1; i<n-1; ++i)
        area+=Cross(p[i]-p[0], p[i+1]-p[0]);
    return area/2;
}
//左上角 逆时针平面
int HalfplaneIntersection(Line *L,int n,Point *poly){
    sort(L, L+n);
    int first,last;
    Point *p=new Point[n];
    Line *q=new Line[n];
    q[first=last=0]=L[0];
    for (int i=1; i<n; ++i) {
        while (first<last&&!OnLeft(L[i], p[last-1])) {
            last--;
        }
        while (first<last&&!OnLeft(L[i], p[first])) {
            first++;
        }
        q[++last]=L[i];
        if(fabs(Cross(q[last].v,q[last-1].v))<eps) {
            --last;
            if(OnLeft(q[last], L[i].p))
                q[last]=L[i];
        }
        if (first<last)
            p[last-1]=GetIntersection(q[last-1], q[last]);
    }
    while (first<last&&!OnLeft(q[first], p[last-1])) {
        --last;
    }
    if (last-first<=1) {
        return 0;
    }
    p[last]=GetIntersection(q[last], q[first]);
    int m=0;
    for (int i=first; i<=last; ++i) {
        poly[m++]=p[i];
    }
    delete []p;
    delete []q;
    return m;
}

double Rotation_Calipers(int len){      //旋转卡壳求凸包最大三角形

    double ans = 0;
    for(int i = 0 ; i < len ; i++){
        int j = (i + 1) % len;
        int k = (j + 1) % len;
        while(j != i && k != i){
            while(dcmp(fabs(Cross(p[j]-p[i],p[(k+1)%len]-p[i])) - fabs(Cross(p[j]-p[i],p[k]-p[i])))>0)
                k = (k + 1) % len;
            ans = max(ans,fabs(Cross(p[j]-p[i],p[k]-p[i])));
            // 最大 三个点 i,j,k
            j = (j + 1) % len;
        }
    }
    return ans;
}
// 求凸包 最远点对
double dis(Point a, Point b)
{
    return (a.x-b.x)*(a.x-b.x) + (a.y-b.y)*(a.y-b.y);
}

double rotating_calipers(Point *p, int n)///旋转卡壳求凸包的直径，平面最远的点对
{
    int k=1;
    double ans = 0;
    p[n]=p[0];
    for(int i=0; i<n; i++)
    {
        while(fabs(Cross(p[i+1]-p[i],p[k]-p[i]))<fabs(Cross(p[i+1]-p[i],p[k+1]-p[i])))
            k=(k+1)%n;
        ans = max(ans,max(dis(p[i],p[k]),dis(p[i+1],p[k])));///注意是用的函数dis
    }
    return ans;
}

//定圆覆盖最大点数
struct alpha
{
    double v;
    bool flag;
    bool operator < (alpha b) const //排序专用偏序关系
    {
        return v < b.v;
    }
} alp[Mn * Mn + 5]; //C(N,2)*2的可能交点（可能极角）


double dis(Point a,Point b)
{
    return sqrt((a.x - b.x)*(a.x - b.x) + (a.y - b.y)*(a.y - b.y));
}

double R;//定长

int solve(int n)
{
    int MAX = 0;
    for(int i = 0; i < n; i++)
    {
        int t = 0;
        for(int j = 0; j < n; j++)
        {
            if(i == j)
                continue;
            double theta,phi,D;
            D = dis(p[i],p[j]);
            if(D > 2.0 * R)
                continue;
            theta = atan2(p[j].y - p[i].y,p[j].x - p[i].x);
            if(theta < 0)
                theta += 2 * pi;
            phi = acos(D / (2.0 * R));
            alp[t].v = theta - phi + 2 * pi;
            alp[t].flag = true;
            alp[t + 1].v = theta + phi + 2 * pi;
            alp[t + 1].flag = false;
            t += 2;
        }
        sort(alp,alp + t);
        int sum = 0;
        for(int j = 0; j < t; j++)
        {
            if(alp[j].flag)
                sum ++;
            else
                sum --;
            if(sum > MAX)
                MAX = sum;
        }
    }
    return MAX+1;
}