const double eps=1e-8;//精度
const int INF=0x3f3f3f3f;
const double PI=acos(-1.0);
inline int dcmp(const double& x) //判斷double等於0或。。。
{
    if(fabs(x)<eps)return 0;
    else return x<0?-1:1;
}
struct Point
{
    double x,y;
    Point() {}
    Point(double x,double y):x(x),y(y) {}
};
typedef Point Vector;
typedef vector<Point> Polygon;
inline Vector operator+(const Vector& a,const Vector& b)
{
    return Vector(a.x+b.x,a.y+b.y);    //向量+向量=向量
}
inline Vector operator-(const Point& a,const Point& b)
{
    return Vector(a.x-b.x,a.y-b.y);    //點-點=向量
}
inline Vector operator*(const Vector& a,const double& p)
{
    return Vector(a.x*p,a.y*p);    //向量*實數=向量
}
inline Vector operator/(const Vector& a,const double& p)
{
    return Vector(a.x/p,a.y/p);    //向量/實數=向量
}
inline bool operator<( const Point& A,const Point& B )
{
    return dcmp(A.x-B.x)<0||(dcmp(A.x-B.x)==0&&dcmp(A.y-B.y)<0);
}
inline bool operator==(const Point&a,const Point&b)
{
    return dcmp(a.x-b.x)==0&&dcmp(a.y-b.y)==0;
}
inline bool operator!=(const Point&a,const Point&b)
{
    return a==b?false:true;
}
struct Segment
{
    Point a,b;
    Segment() {}
    Segment(Point _a,Point _b)
    {
        a=_a,b=_b;
    }
    inline bool friend operator<(const Segment& p,const Segment& q)
    {
        return p.a<q.a||(p.a==q.a&&p.b<q.b);
    }
    inline bool friend operator==(const Segment& p,const Segment& q)
    {
        return (p.a==q.a&&p.b==q.b)||(p.a==q.b&&p.b==q.a);
    }
};
struct Circle
{
    Point c;
    double r;
    Circle() {}
    Circle(Point _c, double _r):c(_c),r(_r) {}
    Point point(double a)const
    {
        return Point(c.x+cos(a)*r,c.y+sin(a)*r);
    }
    bool friend operator<(const Circle& a,const Circle& b)
    {
        return a.r<b.r;
    }
};
struct Line
{
    Point p;
    Vector v;
    double ang;
    Line() {}
    Line(const Point &_p, const Vector &_v):p(_p),v(_v)
    {
        ang = atan2(v.y, v.x);
    }
    inline bool operator<(const Line &L)const
    {
        return  ang < L.ang;
    }
};
inline double Dot(const Vector& a,const Vector& b)
{
    return a.x*b.x+a.y*b.y;    //|a|*|b|*cosθ 點積
}
inline double Length(const Vector& a)
{
    return sqrt(Dot(a,a));    //|a| 向量長度
}
inline double Angle(const Vector& a,const Vector& b)
{
    return acos(Dot(a,b)/Length(a)/Length(b));    //向量夾角θ
}
inline double Cross(const Vector& a,const Vector& b)
{
    return a.x*b.y-a.y*b.x;    //叉積 向量圍成的平行四邊形的面積
}
inline double Area2(const Point& a,const Point& b,Point c)
{
    return Cross(b-a,c-a);    //同上 引數為三個點
}
inline double DegreeToRadius(const double& deg)
{
    return deg/180*PI;
}
inline double GetRerotateAngle(const Vector& a,const Vector& b) //向量a順時針旋轉theta度得到向量b的方向
{
    double tempa=Angle(a,Vector(1,0));
    if(a.y<0) tempa=2*PI-tempa;
    double tempb=Angle(b,Vector(1,0));
    if(b.y<0) tempb=2*PI-tempb;
    if((tempa-tempb)>0) return tempa-tempb;
    else return tempa-tempb+2*PI;
}
inline double torad(const double& deg)
{
    return deg/180*PI;    //角度化為弧度
}
inline Vector Rotate(const Vector& a,const double& rad) //向量逆時針旋轉rad弧度
{
    return Vector(a.x*cos(rad)-a.y*sin(rad),a.x*sin(rad)+a.y*cos(rad));
}
inline Vector Normal(const Vector& a) //計算單位法線
{
    double L=Length(a);
    return Vector(-a.y/L,a.x/L);
}
inline Point GetLineProjection(const Point& p,const Point& a,const Point& b) //點在直線上的投影
{
    Vector v=b-a;
    return a+v*(Dot(v,p-a)/Dot(v,v));
}
inline 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 ConvexHull(Point* p,int n,Point* sol) //計算凸包
{
    sort(p,p+n);
    int m=0;
    for(int i=0; i<n; i++)
    {
        while(m>1&&dcmp(Cross(sol[m-1]-sol[m-2],p[i]-sol[m-2]))<=0) m--;
        sol[m++]=p[i];
    }
    int k=m;
    for(int i=n-2; i>=0; i--)
    {
        while(m>k&&dcmp(Cross(sol[m-1]-sol[m-2],p[i]-sol[m-2]))<=0) m--;
        sol[m++]=p[i];
    }
    if(n>0) m--;
    return m;
}
double Heron(double a,double b,double c) //海倫公式
{
    double p=(a+b+c)/2;
    return sqrt(p*(p-a)*(p-b)*(p-c));
}
bool SegmentProperIntersection(const Point& a1,const Point& a2,const Point& b1,const Point& b2) //線段規範相交判定
{
    double c1=Cross(a2-a1,b1-a1),c2=Cross(a2-a1,b2-a1);
    double c3=Cross(b2-b1,a1-b1),c4=Cross(b2-b1,a2-b1);
    return dcmp(c1)*dcmp(c2)<0&&dcmp(c3)*dcmp(c4)<0;
}
double CutConvex(const int& n,Point* poly,const Point& a,const Point& b, vector<Point> result[3]) //有向直線a b 切割凸多邊形
{
    vector<Point> points;
    Point p;
    Point p1=a,p2=b;
    int cur,pre;
    result[0].clear();
    result[1].clear();
    result[2].clear();
    if(n==0) return 0;
    double tempcross;
    tempcross=Cross(p2-p1,poly[0]-p1);
    if(dcmp(tempcross)==0) pre=cur=2;
    else if(tempcross>0) pre=cur=0;
    else pre=cur=1;
    for(int i=0; i<n; i++)
    {
        tempcross=Cross(p2-p1,poly[(i+1)%n]-p1);
        if(dcmp(tempcross)==0) cur=2;
        else if(tempcross>0) cur=0;
        else cur=1;
        if(cur==pre)
        {
            result[cur].push_back(poly[(i+1)%n]);
        }
        else
        {
            p1=poly[i];
            p2=poly[(i+1)%n];
            p=GetLineIntersection(p1,p2-p1,a,b-a);
            points.push_back(p);
            result[pre].push_back(p);
            result[cur].push_back(p);
            result[cur].push_back(poly[(i+1)%n]);
            pre=cur;
        }
    }
    sort(points.begin(),points.end());
    if(points.size()<2)
    {
        return 0;
    }
    else
    {
        return Length(points.front()-points.back());
    }
}
double DistanceToSegment(Point p,Segment s) //點到線段的距離
{
    if(s.a==s.b) return Length(p-s.a);
    Vector v1=s.b-s.a,v2=p-s.a,v3=p-s.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);
}
inline bool isPointOnSegment(const Point& p,const Segment& s)
{
    return dcmp(Cross(s.a-p,s.b-p))==0&&dcmp(Dot(s.a-p,s.b-p))<0;
}
int isPointInPolygon(Point p, Point* poly,int n) //點與多邊形的位置關係
{
    int wn=0;
    for(int i=0; i<n; i++)
    {
        Point& p2=poly[(i+1)%n];
        if(isPointOnSegment(p,Segment(poly[i],p2))) return -1;//點在邊界上
        int k=dcmp(Cross(p2-poly[i],p-poly[i]));
        int d1=dcmp(poly[i].y-p.y);
        int d2=dcmp(p2.y-p.y);
        if(k>0&&d1<=0&&d2>0)wn++;
        if(k<0&&d2<=0&&d1>0)wn--;
    }
    if(wn) return 1;//點在內部
    else return 0;//點在外部
}
double PolygonArea(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 GetLineCircleIntersection(Line L,Circle C,Point& p1,Point& p2) //圓與直線交點 返回交點個數
{
    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)  //相切
    {
        p1=p1=L.p+L.v*(-f/(2*e));
        return 1;
    }//相交
    p1=(L.p+L.v*(-f-sqrt(delta))/(2*e));
    p2=(L.p+L.v*(-f+sqrt(delta))/(2*e));
    return 2;
}
double rotating_calipers(Point *ch,int n)//旋轉卡殼  返回凸包中距離最大的兩點的距離的平方
{
    int q=1;
    double ans=0;
    ch[n]=ch[0];
    for(int p=0; p<n; p++)
    {
        while(Cross(ch[q+1]-ch[p+1],ch[p]-ch[p+1])>Cross(ch[q]-ch[p+1],ch[p]-ch[p+1]))
            q=(q+1)%n;
        ans=max(ans,max(Length(ch[p]-ch[q]),Length(ch[p+1]-ch[q+1])));
    }
    return ans;
}
Polygon CutPolygon(Polygon poly,const Point& a,const Point& b) //用a->b切割多邊形 返回左側
{
    Polygon newpoly;
    int n=poly.size();
    for(int i=0; i<n; i++)
    {
        Point c=poly[i];
        Point d=poly[(i+1)%n];
        if(dcmp(Cross(b-a,c-a))>=0) newpoly.push_back(c);
        if(dcmp(Cross(b-a,c-d))!=0)
        {
            Point ip=GetLineIntersection(a,b-a,c,d-c);
            if(isPointOnSegment(ip,Segment(c,d))) newpoly.push_back(ip);
        }
    }
    return newpoly;
}
int GetCircleCircleIntersection(Circle c1,Circle c2,Point& p1,Point& p2) //求兩圓相交
{
    double d=Length(c1.c-c2.c);
    if(dcmp(d)==0)
    {
        if(dcmp(c1.r-c2.r)==0) return -1;//兩圓重合
        return 0;
    }
    if(dcmp(c1.r+c2.r-d)<0) return 0;
    if(dcmp(fabs(c1.r-c2.r)-d)>0) return 0;
    double a=Angle(c2.c-c1.c,Vector(1,0));
    double 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;
}
inline bool isPointOnleft(Point p,Line L)
{
    return dcmp(Cross(L.v,p-L.p))>0;    //點在直線左邊 線上不算，算的話要>=
}
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&&!isPointOnleft(p[last-1],L[i])) last--;
        while(first<last&&!isPointOnleft(p[first],L[i])) first++;
        q[++last]=L[i];
        if(dcmp(Cross(q[last].v,q[last-1].v))==0)
        {
            last--;
            if(isPointOnleft(L[i].p,q[last])) q[last]=L[i];
        }
        if(first<last) p[last-1]=GetLineIntersection(q[last-1].p,q[last-1].v,q[last].p,q[last].v);
    }
    while(first<last&&!isPointOnleft(p[last-1],q[first])) last--;
    if(last-first<=1) return 0;
    p[last]=GetLineIntersection(q[last].p,q[last].v,q[first].p,q[first].v);
    int m=0;
    for(int i=first; i<=last; i++) poly[m++]=p[i];
    return m;
}
//兩點式化為一般式A = b.y-a.y, B = a.x-b.x, C = -a.y*(B)-a.x*(A);,點與圓的切線、圓與圓的共切線、三角形的外接圓和內接圓在白皮書p266