计算几何 模板
大量功能未(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;
}