计算几何 模板

Posted 2020-10-26 acha

大量功能未(lan)经(de)测试等待后续测试

好吧是不知道怎么测试....

先把模板备份在这里.

大佬们心情好的话可以帮忙查查错呀

#include <cstdio>
#include <cstdlib>
#include <cctype>
#include <cmath>
#include <cstring>
#include <algorithm>
#include <cassert>
#include <vector>
#define ri rd<int>
#define rep(i, a, b) for (int i = (a), _ = (b); i <= _; ++ i)
#define per(i, a, b) for (int i = (a), _ = (b); i >= _; -- i)
#define For(i, a, b) for (int i = (a), _ = (b); i < _; ++ i)
#define forea(it, v) for (__typeof(v.begin()) it = v.begin(); it != v.end(); ++ it)
using namespace std;
typedef long double db;
const db pi = acos(-1.);
const db eps = 1e-10;

template<class T> T rd() {
    bool f = 1; char c = getchar(); for (; !isdigit(c); c = getchar()) if (c == ‘-‘) f = 0;
    T x = 0; for (; isdigit(c); c = getchar()) x = x * 10 + c - 48; return f ? x : -x;
}

int sgn(db x) { return x < -eps ? -1 : x > eps; }
bool le(db x, db y) { return sgn(y-x) != -1; }
bool lt(db x, db y) { return sgn(y-x) == 1; }
bool ge(db x, db y) { return sgn(x-y) != -1; }
bool gt(db x, db y) { return sgn(x-y) == 1; }
bool eq(db x, db y) { return sgn(x-y) == 0; }

struct Vector;
typedef Vector Point;
struct Vector {
    db x, y;
    Vector(db _x = 0, db _y = 0):x(_x), y(_y) { }

    friend Vector operator + (const Vector &a, const Vector &b) { 
        return Vector(a.x + b.x, a.y + b.y); 
    }

    friend Vector operator - (const Vector &a, const Vector &b) { 
        return Vector(a.x - b.x, a.y - b.y); 
    }

    friend Vector operator * (const Vector &a, const db k) { 
        return Vector(a.x * k, a.y * k); 
    }

    friend Vector operator * (const db k, const Vector &a) {
        return Vector(a.x * k, a.y * k); 
    }

    friend Vector operator / (const Vector &a, const db k) {
        return Vector(a.x / k, a.y / k); 
    }

    Vector operator - () { // 反向
        return Vector(-x, -y); 
    }

    friend bool operator < (const Point &a, const Point &b) { // 凸包排序用
        return eq(a.x, b.x) ? lt(a.y, b.y) : lt(a.x, b.x); 
    }

    db len() const { // 模长
        return sqrt(x * x + y * y); 
    }

    Vector unit() const { // 单位向量
        return *this / len(); 
    }

    Vector normal() const { // 法向量
        return Vector(-y, x);
    }
};

db dot(const Vector &a, const Vector &b) { 
    return a.x * b.x + a.y * b.y; 
}

db det(const Vector &a, const Vector &b) { 
    return a.x * b.y - a.y * b.x; 
}

db p_dis_p(const Point&, const Point&);

Point p_middle_p(const Point&, const Point&);

struct Line { // 直线
    Point a, b; // 两点确定一线, 方向向量定义为 a->b
    Line(Point _a, Point _b):a(_a), b(_b) { }

    bool chk_p(const Point &p) const {
        return 1;
    }
};

struct Seg:Line { // 线段
    Seg(Point _a, Point _b):Line(_a, _b) {}

    bool chk_p(const Point &p) const {
        return le(dot(p - a, p - b), 0);
    }

    Point middle() const { // 线段中点
        return (a + b) / 2;
    }
    
    Line mperp() const { // 中垂线
        Point u = middle();
        return Line(u, u + (b - a).normal());
    }

    db length() const {
        return p_dis_p(a, b);
    }
};

struct Arrow:Line { // 射线
    Arrow(Point _a, Point _b):Line(_a, _b) {}

    bool chk_p(const Point &p) const {
        return ge(dot(p - a, b - a), 0);
    }
};

struct Circle {
    Point o; db r;
    Circle() {}
    Circle(Point _o, db _r):o(_o), r(_r) { } // 点径型构造函数 
    Circle(Point a, Point b) { // 直径型构造函数
        o = p_middle_p(a, b);
        r = p_dis_p(a, o);
    } 

    bool chk_p(const Point &p) const {
        return eq(p_dis_p(p, o), r);
    }
};

struct Arc:Circle {
    Point a, b; // a 逆时针转至 b
    Arc(Point _o, db _r, Point _a, Point _b):Circle(_o, _r), a(_a), b(_b) { } 
    Arc(Circle _c, Point _a, Point _b):a(_a), b(_b) { o = _c.o, r = _c.r; } 

    bool chk_p(const Point &p) const {
        return eq(p_dis_p(p, o), r) && le(dot(p - a, p - b), 0);
    }

    Arc bu() const { // 补弧
        return Arc(o, r, b, a);
    }

    db bad() const { // 劣弧 | 半圆 ?
        return ge(det(a - o, b - o), 0);
    }

    db length() const { // 弧长
        if (bad()) {
            Point u = p_middle_p(a, b);
            db x = p_dis_p(o, u);
            db y = p_dis_p(a, u);
            return atan2(y, x) * r;
        }
        else return 2 * pi * r - bu().length();
    }

    db shan() const { // 扇形
        return length() * r / 2;
    }

    db gong() const { // 弓形
        Point u = p_middle_p(a, b);
        db x = p_dis_p(o, u);
        db y = p_dis_p(a, u);
        return atan2(y, x) * r * r / 2 - x * y;
    }
};

template<class T1, class T2> bool l_parallel_l(const T1 a, const T2 b) {
    return eq(det(a.b - a.a, b.b - b.a), 0);
}

Point p_middle_p(const Point &a, const Point &b) {
    return (a + b) / 2;
}

template<class T> Point p_projection_l(const Point &p, const T &l) {
    Vector v = (l.b - l.a).unit();
    return l.a + dot(p - l.a, v) * v;
}

template<class T> Seg seg_projection_l(const Seg &s, const T &l) {
    return Seg(p_projection_l(s.a, l), p_projection_l(s.b, l));
}

Point p_over_p(const Point &a, const Point &b) {
    return 2 * b - a; 
}

template<class T> Point p_over_l(const Point &p, const T &l) {
    return p_over_p(p_projection_l(p, l));
}

db p_dis_p(const Point &a, const Point &b) { 
    return (b - a).len(); 
}

template<class T> db p_dis_l(const Point &p, const T &l) {
    Point q = p_projection_l(p, l);
    if (l.chk_p(q)) return p_dis_p(p, q);
    else return min(p_dis_p(p, l.a), p_dis_p(p, l.b));
}

int c_relative_c(Circle a, Circle b) { // 代码中统一不考虑两圆重合的情况
    if (lt(a.r, b.r)) swap(a, b);
    db d = p_dis_p(a.o, b.o);
    if (gt(d, a.r + b.r)) return 2; // 相离
    if (eq(d, a.r + b.r)) return 1; // 外切
    if (eq(d, a.r - b.r)) return -1; // 内切
    if (lt(d, a.r - b.r)) return -2; // 内含
    return 0; // 相交
}

#define PUT(x) if (1) { Point _ = x; if (a.chk_p(_) && b.chk_p(_)) S.push_back(_); }
template<class T1, class T2> vector<Point> l_intersection_l(const T1 &a, const T2 &b) { // 代码中统一不考虑两线在同一直线上的情况
    vector<Point> S;
    if (l_parallel_l(a, b)) return S;
    db s1 = det(b.a - a.a, b.b - a.a);
    db s2 = det(b.b - a.b, b.a - a.b);
    PUT(a.a + (a.b - a.a) * s1 / (s1 + s2)); return S;
}

template<class T1, class T2> vector<Point> l_intersection_c(const T1 &a, const T2 &b) {
    vector<Point> S;
    Point u = p_projection_l(b.o, a);
    db x = p_dis_p(b.o, u);
    if (lt(b.r, x)) return S;
    if (eq(b.r, x)) { PUT(u); return S; }
    Vector v = (a.b - a.a).unit();
    db y = sqrt(b.r * b.r - x * x);
    PUT(u + v * y); PUT(u - v * y); return S;
}

template<class T1, class T2> vector<Point> c_intersection_c(T1 a, T2 b) {
    if (lt(a.r, b.r)) swap(a, b);
    vector<Point> S;
    db d = p_dis_p(a.o, b.o);
    int kd = c_relative_c(a, b);
    if (abs(kd) == 2) return S;
    else if (abs(kd) == 1) { PUT(a.o + (b.o - a.o).unit() * a.r); return S; }
    db x = fabs((d + (a.r * a.r - b.r * b.r) / d) / 2); // 分类讨论两圆相交的各种情况, 式子相同, 符号不一
    db z = sqrt(a.r * a.r - x * x);
    Vector u = (b.o - a.o).unit(), v = u.normal();
    Point c = a.o + u * x;
    PUT(c + v * z); PUT(c - v * z); return S;
}
#undef PUT

vector<Point> p_tangentp_c(const Point &p, const Circle &c) {
    Circle b(p, c.o);
    return c_intersection_c(b, c);
}

vector<Line> p_tangentl_c(const Point &p, const Circle &c) {
    vector<Point> S = p_tangentp_c(p, c);
    vector<Line> T; forea (it, S) T.push_back(Line(*it, p)); return T;
}

// 求两圆的共切点的话, 四个点关系也不知道, 意义不大, 直接求线
vector<Line> c_itangentl_c(Circle a, Circle b) { // 内共切线
    vector<Line> T;
    if (lt(a.r, b.r)) swap(a, b);
    int kd = c_relative_c(a, b);
    if (kd <= 0) return T;
    if (kd == 1) {
        Point u = a.o + (b.o - a.o).unit() * a.r;
        T.push_back(Line(u, u + (b.o - a.o).normal()));
        return T;
    }
    Circle c(a.o, a.r + b.r);
    Circle d(a.o, b.o);
    vector<Point> S = c_intersection_c(c, d); 
    forea (it, S) {
        Point u = a.o + (*it - a.o).unit() * a.r;
        T.push_back(Line(u, u + (b.o - *it)));
    }
    return T;
}

vector<Line> c_extangentl_c(Circle a, Circle b) { // 外公切线
    vector<Line> T;
    if (eq(a.r, b.r)) {
        Point u = a.o + (b.o - a.o).unit().normal() * a.r;
        T.push_back(Line(u, u + (b.o - a.o)));
        u = a.o - (b.o - a.o).unit().normal() * a.r;
        T.push_back(Line(u, u + (b.o - a.o)));
        return T;
    }
    if (lt(a.r, b.r)) swap(a, b);
    int kd = c_relative_c(a, b);
    if (kd == -2) return T;
    if (kd == -1) {
        Point u = a.o + (b.o - a.o).unit() * a.r;
        T.push_back(Line(u, u + (b.o - a.o).normal()));
        return T;
    }
    Circle c(a.o, a.r - b.r);
    Circle d(a.o, b.o);
    vector<Point> S = c_intersection_c(c, d); 
    forea (it, S) {
        Point u = a.o + (*it - a.o).unit() * a.r;
        T.push_back(Line(u, u + (b.o - *it)));
    }
    return T;
}

int main() {
// test data
    Circle a(Point(10, 10), 10);
    Circle b(Point(16, 14), 11.313708498984761);

    int kd = c_relative_c(a, b);
    kd = c_relative_c(b, a);

    vector<Point> S = c_intersection_c(a, b);
    S = c_intersection_c(b, a);

    vector<Line> T1 = c_extangentl_c(a, b);
    T1 = c_extangentl_c(b, a);

    vector<Line> T2 = c_itangentl_c(a, b);
    T2 = c_itangentl_c(b, a);

    return 0;
}