计算几何 模板

大量功能未(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;
}