计算几何模版

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<vector>
#include<algorithm>
#include<cassert>
#define PI acos(-1.0)
#define inf 0x3f3f3f3f
using namespace std;
const double eps = 1e-6;
struct Point{
    double x,y;
    Point(double tx=0,double ty=0):x(tx),y(ty){}
};
typedef Point Vctor;

//向量的加减乘除
Vctor operator + (Vctor A,Vctor B){return Vctor(A.x+B.x,A.y+B.y);}
Vctor operator - (Point A,Point B){return Vctor(A.x-B.x,A.y-B.y);}
Vctor operator * (Vctor A,double p){return Vctor(A.x*p,A.y*p);}
Vctor operator / (Vctor A,double p){return Vctor(A.x/p,A.y/p);}
bool operator < (Point A,Point B){return A.x < B.x || (A.x == B.x && A.y < B.y);}

struct Line{
    Point p;
    Vctor v;
    Line(Point p=Point(0,0),Vctor v=Vctor(0,0)):p(p),v(v){}
    Point point(double t){return p + v*t;} //获得直线上的距离p点t个单位长度的点
};
struct Circle{
    Point c;
    double r;
    Circle(Point tc=Point(0,0),double tr=0):c(tc),r(tr){}
    Point point(double a){return Point( c.x + cos(a)*r , c.y + sin(a)*r);}
};

int dcmp(double x)
{
    if(fabs(x)<eps) return 0;
    else return (x<0)?(-1):(1);
}
bool operator == (Point A,Point B){return dcmp(A.x-B.x)==0 && dcmp(A.y-B.y)==0;}

//向量的点积，长度，夹角
double Dot(Vctor A,Vctor B){return A.x*B.x+A.y*B.y;}
double Length(Vctor A){return sqrt(Dot(A,A));}
double Angle(Vctor A,Vctor B){return acos(Dot(A,B)/Length(A)/Length(B));}

//叉积，三角形面积
double Cross(Vctor A,Vctor B){return A.x*B.y-A.y*B.x;}
double TriangleArea(Point A,Point B,Point C){return Cross(B-A,C-A);}

//向量的旋转,求向量的单位法线（即左转90度，然后长度归一）
Vctor Rotate(Vctor A,double rad){return Vctor( A.x*cos(rad) - A.y*sin(rad) , A.x*sin(rad) + A.y*cos(rad) );}
Vctor Normal(Vctor A)
{
    double L = Length(A);
    return Vctor(-A.y/L, A.x/L);
}

//直线的交点
Point getLineIntersection(Line L1,Line L2)
{
    Vctor u = L1.p-L2.p;
    double t = Cross(L2.v,u)/Cross(L1.v,L2.v);
    return L1.p + L1.v*t;
}

//点到直线的距离
double DistanceToLine(Point P,Line L)
{
    return fabs(Cross(P-L.p,L.v))/Length(L.v);
}
double DistanceToLine(Point P,Point A,Point B){//点P到线段A，B的距离
    Vctor v1=B-A;
    Vctor v2=P-A;
    return fabs(Cross(v1,v2))/Length(v1);
}


//点到线段的距离
double DistanceToSegment(Point P,Point A,Point B)
{
    if(A==B) return Length(P-A);
    Vctor 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);
}
double DistanceTosegment(Point P,Point A,Point B){//点P到射线A,B的距离
    if(A==B)return Length(P-A);
    Vctor v1=B-A;
    Vctor 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);
}


//点到直线的映射
Point getLineProjection(Point P,Line L)
{
    return L.v + L.v*Dot(L.v,P-L.p)/Dot(L.v,L.v);
}
Point GetLineProjection(Point P,Point A,Point B){//点P到线A，B的投影的点
    Vctor v=B-A;
    return A+v*(Dot(v,P-A)/Dot(v,v));
}


//判断线段是否规范相交
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;
}

//判断点是否在一条线段上
bool OnSegment(Point P,Point a1,Point a2)
{
    return dcmp(Cross(a1 - P,a2 - P))==0 && dcmp(Dot(a1 - P,a2 - P))<0;
}

//多边形面积
double PolgonArea(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,vector<Point> &sol)
{
    double t1,t2;
    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.0*e*g;
    if(dcmp(delta)<0) return 0;
    else if(dcmp(delta)==0)
    {
        t1 = t2 = -f/(2.0*e);
        sol.push_back(L.point(t1));
        return 1;
    }
    else
    {
        t1 = (-f-sqrt(delta))/(2.0*e); sol.push_back(L.point(t1));
        t2 = (-f+sqrt(delta))/(2.0*e); sol.push_back(L.point(t2));
        return 2;
    }
}

//判断并求出两圆的交点
double angle(Vctor v){return atan2(v.y,v.x);}
int getCircleIntersection(Circle C1,Circle C2,vector<Point> &sol)
{
    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.point(a - da), p2 = C1.point(a + da);
    sol.push_back(p1);
    if(p1==p2) return 1;
    sol.push_back(p2);
    return 2;
}

//求点到圆的切线
int getTangents(Point p,Circle C,vector<Line> &sol)
{
    Vctor u=C.c-p;
    double dis=Length(u);
    if(dis<C.r) return 0;
    else if(dcmp(dis-C.r) == 0)
    {
        sol.push_back(Line(p,Rotate(u,M_PI/2)));
        return 1;
    }
    else
    {
        double ang=asin(C.r/dis);
        sol.push_back(Line(p,Rotate(u,-ang)));
        sol.push_back(Line(p,Rotate(u,ang)));
        return 2;
    }
}

//求两圆的切线
int getCircleTangents(Circle A,Circle B,Point *a,Point *b)
{
    int cnt = 0;
    if(A.r<B.r){swap(A,B);swap(a,b);}
    //圆心距的平方
    double d2 = (A.c.x - B.c.x)*(A.c.x - B.c.x) + (A.c.y - B.c.y)*(A.c.y - B.c.y);
    double rdiff = A.r - B.r;
    double rsum = A.r + B.r;
    double base = angle(B.c - A.c);
    //重合有无限多条
    if(d2 == 0 && dcmp(A.r - B.r) == 0) return -1;
    //内切
    if(dcmp(d2-rdiff*rdiff) == 0)
    {
        a[cnt] = A.point(base);
        b[cnt] = B.point(base);
        cnt++;
        return 1;
    }
    //有外公切线
    double ang = acos((A.r - B.r) / sqrt(d2));
    a[cnt] = A.point(base + ang); b[cnt] = B.point(base + ang); cnt++;
    a[cnt] = A.point(base - ang); b[cnt] = B.point(base - ang); cnt++;
    
    //一条内切线
    if(dcmp(d2-rsum*rsum) == 0)
    {
        a[cnt] = A.point(base);
        b[cnt] = B.point(M_PI + base);
        cnt++;
    }//两条内切线
    else if(dcmp(d2-rsum*rsum) > 0)
    {
        double ang = acos((A.r + B.r) / sqrt(d2));
        a[cnt] = A.point(base + ang); b[cnt] = B.point(base + ang); cnt++;
        a[cnt] = A.point(base - ang); b[cnt] = B.point(base - ang); cnt++;
    }
    return cnt;
}
//求外切圆,三点不能共线
Circle CircumscribedCircle(Point p1,Point p2,Point p3){
    double Bx=p2.x-p1.x,By=p2.y-p1.y;
    double Cx=p3.x-p1.x,Cy=p3.y-p1.y;
    double D=2*(Bx*Cy-By*Cx);
    double cx=(Cy*(Bx*Bx+By*By)-By*(Cx*Cx+Cy*Cy))/D+p1.x;
    double cy=(Bx*(Cx*Cx+Cy*Cy)-Cx*(Bx*Bx+By*By))/D+p1.y;
    Point p=Point(cx,cy);
    return Circle(p,Length(p1-p));
}//求内切圆,三点不能共线
Circle InscribedCircle(Point p1,Point p2,Point p3){
    double a=Length(p2-p3);
    double b=Length(p3-p1);
    double c=Length(p1-p2);
    Point p=(p1*a+p2*b+p3*c)/(a+b+c);
    return Circle(p,DistanceToLine(p, p1, p2));
}