计算几何模板（未完待续）

目前基本都是从蓝书上摘录的。

有一部分需要线性代数的知识，但是蓝书作者并没有解释，个人觉得用数学知识推出来更有助于记忆，死记硬背板子容易忘。以后有机会的话我在这里写点注解。

二维基础操作：

struct Point
{
    double x, y;
    Point(double x = 0, double y = 0):x(x), y(y){}
};

typedef Point Vector;

int dcmp(double x)//比较
{
    const double eps = 1e-10;
    if (fabs(x) < eps)    return 0;
    else    return x < 0 ? -1 : 1;
}
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); }
bool operator < (const Point& A, const Point& B) { return A.x < B.x || (A.x == B.x && A.y < B. y); }
bool operator == (const Point& A, const Point& B) { return dcmp(A.x - B.x) == 0 && dcmp(A.y - B.y) == 0; }

const double PI = acos(-1.0);
double torad(double deg) { return deg/180 * PI; }//角度转弧度

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)); }//角度

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 Cross(Vector A, Vector B) { return A.x * B.y - A.y * B.x; }//叉积

double Area2(Point A, Point B, Point C) { return Cross(B-A, C-A);}//二倍的三角形面积

Vector Normal(Vector A) { double L = Length(A); return Vector(-A.y/L, A.x/L); }//向量的单位法线，旋转的特殊情况

Vector Get_Line_Projection(Point P, Point A, Point B)//P在AB方向上的映射
{
     Vector v = B - A;
      return A + v*(Dot(v, P-A) / Dot(v, v));
}

Vector Get_Line_Intersect(Vector A, Vector v, Vector B, Vector w)//两直线交点
{
    Vector u = A - B;
    double t = Cross(w, u) / Cross(v, w);
    return A + v*t;
}

double Distance_to_Line(Point P, Point A, Point B)//点P到直线AB的距离
{
    Vector v1 = B - A, v2 = P - A;
    return fabs(Cross(v1, v2)) / Length(v1);
}

double Distance_to_Segment(Point P, Point A, Point B)//点P到线段AB的距离
{
    Vector v1 = B - A, v2 = P - A, v3 = P - B;
    if (A == B || dcmp(Dot(v1, v2)) < 0)    return Length(v2);//PA
    else if (dcmp(Dot(v1, v3)) > 0)    return Length(v3);//PB
    else    return fabs(Cross(v1, v2) / Length(v1));//PQ
}

bool On_Segment(Point P, Point A, Point B) { return dcmp(Cross(A - P, B - P)) == 0 && dcmp(Dot(A - P, B - P)) < 0; }//点在线段上

bool Segment_Proper_Intersection(Point a1, Point a2, Point b1, Point b2)//线段a1a2与b1b2相交（不包含端点）
{
    double c1 = Cross(a2 - a1, b1 - a1), c2 = Cross(a2 - a1, b2 - a1);
    double c3 = Cross(b2 - b1, a1 - b1), c2 = Cross(b2 - b1, a2 - b1);
    return dcmp(c1)*dcmp(c2) < 0 && dcmp(c3)*dcmp(c4) < 0;
}

double Polygon_Area(Point *p, int n)//多边形面积，逆时针取点0~n-1
{
    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;
}

//基于复数的几何运算，效率不是很高
#include <complex>
typedef complex<double> Point;
typedef Point Vector;

double Dot(Vector A, Vector B) { return real(conj(A) * B); }// A（x+yi）的共轭（x-yi）乘以B，取实部
double Cross(Vector A, Vector B) { return imag(conj(A) * B); }//取虚部
Vector Rotate(Vector A, double rad) { return A * exp(Point(0, rad)); }//使用了欧拉公式：e^(0 + rad*i) = cos(rad) + sin(rad)*i，求完确实是旋转

圆相关（他那个两圆公切线看着有点奇怪先不贴了）：

struct Circle
{
    Point c;
    double r;
    // Circle(Point a = (0, 0), double x = 0):c(a), r(x){}
    Point point(double seta) { return Point(c.x + cos(seta)*r, c.y + sin(seta)*r); }//已知角度求圆上某点
};

struct Line
{
    Point p;
    Vector v;
    Line(Point p, Vector v):p(p), v(v) { }

    Point point(double t) { return p + v*t; }//已知参数t求直线上一点
    Line move(double d) { return Line(p + Normal(v)*d, v); }//平移
};

int get_Tangents(Point P, Circle C, vector<Vector> &v)//过定点做圆的切线
{
    Vector u = C.c - P;
    double d = Length(u);

    if (dcmp(d - C.r) < 0)    return 0;//点在圆内部
    else if (dcmp(d - C.r) == 0)//点在圆上
    {
        v.push_back(Rotate(u, PI/2));
        return 1;
    }
    else//两条切线
    {
        double a = asin(C.r / d);
        v.push_back(Rotate(u, a));
        v.push_back(Rotate(u, -a));
        return 2;
    }
}

int get_Line_Circle_Intersection(Line L, Circle C, vector<Point> &ans)//方程求解直线与圆的交点
{
    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*e*g;//判别式
    
    if (dcmp(delta) < 0)    return 0;//相离
    else if (dcmp(delta) == 0)//相切
    {
        t1 = t2 = -f/2/e;
        ans.push_back(L.point(t1));
        return 1;
    }
    else//相交
    {
        t1 = (-f + sqrt(delta)) / 2 / e;
        t2 = (-f - sqrt(delta)) / 2 / e;
        ans.push_back(L.point(t1)), ans.push_back(L.point(t2));
        return 2;
    }
}

inline double angle(Vector A) { return atan2(A.y, A.x); }//某向量极角

int get_Circle_Circle_Intersection(Circle C1, Circle C2, vector<Point> &v)//圆与圆的交点
{
    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);
    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);

    v.push_back(p1);
    if (p1 == p2)    return 1;
    v.push_back(p2);
    return 2;
}