几何模板

double eps=1e-9;
const double pi=acos(-1);

struct point{
    double x,y;
    point(){}
    point(double x,double y):x(x),y(y){}
};
typedef point Vector;
typedef point Point;
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);
}
int dcmp(double x){//精度三态函数(>0,<0,=0)
    if (fabs(x) < eps)return 0;
    else if (x > 0)return 1;
    return -1;
}
bool operator == (const Point &a, const Point &b){//向量相等
    return dcmp(a.x - b.x) == 0 && dcmp(a.y - b.y) == 0;
}
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 Cross(Vector a, Vector b){//外积
    return a.x*b.y - a.y*b.x;
}
Vector Rotate(Vector a, double rad){//逆时针旋转 易损失精度
    return Vector(a.x*cos(rad) - a.y*sin(rad), a.x*sin(rad) + a.y*cos(rad));
}
double Distance(Point a, Point b){//两点间距离
    return sqrt((a.x - b.x)*(a.x - b.x) + (a.y - b.y)*(a.y - b.y));
}
double Area(Point a, Point b, Point c){//三角形面积
    return fabs(Cross(b - a, c - a) / 2);
}
bool Intersect(Point a, Point b, Point c, Point d){//线段相交（不包括端点）
    double t1 = Cross(c - a, d - a)*Cross(c - b, d - b);
    double t2 = Cross(a - c, b - c)*Cross(a - d, b - d);
    return dcmp(t1) < 0 && dcmp(t2) < 0;
}
bool StrictIntersect(Point a, Point b, Point c, Point d){ //线段相交（包括端点）
    return
        dcmp(max(a.x, b.x) - min(c.x, d.x)) >= 0
        && dcmp(max(c.x, d.x) - min(a.x, b.x)) >= 0
        && dcmp(max(a.y, b.y) - min(c.y, d.y)) >= 0
        && dcmp(max(c.y, d.y) - min(a.y, b.y)) >= 0
        && dcmp(Cross(c - a, d - a)*Cross(c - b, d - b)) <= 0
        && dcmp(Cross(a - c, b - c)*Cross(a - d, b - d)) <= 0;
}
double DistanceToLine(Point A, Point M, Point N){//点A到直线MN的距离,Error:MN=0
    return fabs(Cross(A - M, A - N) / Distance(M, N));
}
point intersection(point u1,point u2,point v1,point v2)//两直线的交点
{
    point ret=u1;
    double t=((u1.x-v1.x)*(v1.y-v2.y)-(u1.y-v1.y)*(v1.x-v2.x))
             /((u1.x-u2.x)*(v1.y-v2.y)-(u1.y-u2.y)*(v1.x-v2.x));
    ret.x+=(u2.x-u1.x)*t;
    ret.y+=(u2.y-u1.y)*t;
    return ret;
}
point symmetric_point(point p1, point l1, point l2)//点关于直线的对称点
{
    point u=l1-l2;
    u=point(u.y,-u.x);
    point N=p1+u/Length(u)*DistanceToLine(p1,l1,l2)*2.0;
    if(Cross(l1-N,l2-N)*Cross(l1-p1,l2-p1)<eps)return N;
    return p1+u/Length(u)*DistanceToLine(p1,l1,l2)*(-2.0);
}

//求多边形面积（叉积和计算法）
/*模板说明：P[]为多边形的所有顶点，下标为0~n-1，n为多边形边数*/
/*
Point P[1005];
int n;
double PolygonArea(){
    double sum = 0;
    Point O = Point(0, 0);
    for (int i = 0; i < n; i++)
        sum += Cross(P[i] - O, P[(i + 1) % n] - O);
    if (sum < 0)sum = -sum;
    return sum / 2;
}
*/

//求凸包
/*
point p[200009],ans[200009];
int n,top;
point Tmp;//选好的起点
bool cmp(point a,point b){
    F ans=Cross(a-Tmp,b-Tmp);
    if(dcmp(ans)==0)return dcmp(Distance(a,Tmp)-Distance(b,Tmp))<0;
    return ans>0;//表示a到b是逆时针转
}
void Graham(point p[],int n,point ans[],int &top){
    for(int i=2;i<=n;i++){
        if(p[i].y<p[1].y||(dcmp(p[i].y-p[1].y)==0&&p[i].x<p[1].x))
            swap(p[1],p[i]);
    }
    Tmp=p[1];
    sort(p+2,p+1+n,cmp);
    ans[1]=p[1],ans[2]=p[2],top=2;
    for(int i=3;i<=n;i++){
        while(Cross(ans[top]-ans[top-1],p[i]-ans[top])<0)top--;
        ans[++top]=p[i];
    }
    //最后几个点如果是和p[1]在同一直线上的，那么会漏掉除最远的点以外的点
    for(int i=n-1;;i--){
        if(Cross(ans[top]-ans[1],p[i]-ans[1])!=0)break;
        ans[++top]=p[i];
    }
}
*/


struct Circle{
    Point c;
    double r;
    Point gpoint(double a){//基于圆心角求圆上一点坐标
        return Point(c.x + cos(a)*r, c.y + sin(a)*r);
    }
};
double Angle(Vector v1){
    if (v1.y >= 0)return Angle(v1, Vector(1.0, 0.0));
    else return 2 * pi - Angle(v1, Vector(1.0, 0.0));
}
/*
point a,b;
int GetCC(Circle C1, Circle C2){//求两圆交点
    double d = Length(C1.c - C2.c);
    if (dcmp(d) == 0)
    {
        if (dcmp(C1.r - C2.r) == 0)return -1;//重合
        else 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);
    double da = acos((C1.r*C1.r + d * d - C2.r*C2.r) / (2 * C1.r*d));
    Point p1 = C1.gpoint(a - da), p2 = C1.gpoint(a + da);
    if (p1 == p2){a=b=p1;return 1;}
    else {a=p1,b=p2;return 2;}
}
*/