「计算几何」模板
基础模板。
#include <bits/stdc++.h>
#define fr friend
const double PI = acos(-1);
const double EPS = 1E-9;
const double INF = 1E9;
const int N = 100000;
int dcmp(double x) {return (fabs(x) <= EPS ? 0 : (x > 0 ? 1 : -1));}
// 这里还是建议写成 fabs(x) <= EPS,因为当 EPS = 0 时良定义
struct point{
double x, y;
fr point operator + (const point &a, const point &b) {return (point){a.x + b.x, a.y + b.y};}
fr point operator - (const point &a, const point &b) {return (point){a.x - b.x, a.y - b.y};}
fr point operator * (const point &a, double k) {return (point){a.x * k, a.y * k};}
fr point operator / (const point &a, double k) {return (point){a.x / k, a.y / k};}
fr void operator += (point &a, const point &b) {a = a + b;}
fr void operator -= (point &a, const point &b) {a = a - b;}
fr void operator *= (point &a, double k) {a = a * k;}
fr void operator /= (point &a, double k) {a = a / k;}
fr double operator * (const point &a, const point &b) {return a.x * b.x + a.y * b.y;} // dot
fr double operator ^ (const point &a, const point &b) {return a.x * b.y - a.y * b.x;} // cross
fr bool operator == (const point &a, const point &b) {return dcmp(a.x - b.x) == 0 && dcmp(a.y - b.y) == 0;}
fr bool operator < (const point &a, const point &b) {
return dcmp(a.x - b.x) == 0 ? dcmp(a.y - b.y) < 0 : dcmp(a.x - b.x) < 0;
}
fr double angle(const point &a) {return atan2(a.y, a.x);}
fr double length(const point &a) {return sqrt(a * a);}
fr double dist(const point &a, const point &b) {return length(a - b);}
fr double area(const point &a, const point &b, const point &c) {return fabs((b - a) ^ (c - a)) / 2;}
fr double area(point *a, int n) {
double ret = a[0] ^ a[n - 1];
for(int i=1;i<n;i++) ret += a[i] ^ a[i - 1];
return fabs(ret) / 2;
}
fr point normal(const point &a) {return (point){a.y, -a.x};}
fr point unit(const point &a) {return a / length(a);}
fr point bisector(const point &a, const point &b) {return a * length(b) + b * length(a);}
fr point rotate(const point &a, double t) {
return (point){cos(t)*a.x - sin(t)*a.y, cos(t)*a.y + sin(t)*a.x};
}
fr void read(point &p) {scanf("%lf%lf", &p.x, &p.y);}
};
struct line{
point a, ab; double ang; // a + ab*t
line() {} line(point _a, point _b) : a(_a), ab(_b - _a), ang(angle(ab)) {}
point b() const &{return a + ab;}
fr bool operator < (const line &a, const line &b) {return a.ang < b.ang;}
fr bool on_left(const line &l, const point &p) {return dcmp(l.ab ^ (p - l.a)) > 0;}
fr bool on_line(const line &l, const point &p) {return dcmp((p - l.a) ^ (p - l.b())) == 0;}
fr bool on_segment(const line &l, const point &p) {
return on_line(l, p) && dcmp((p - l.a) * l.ab) >= 0 && dcmp((p - l.b()) * l.ab) <= 0;
}
fr bool same_side(const line &l, const point &p1, const point &p2) {
return dcmp(((p1 - l.a) ^ (p1 - l.b())) * ((p2 - l.a) ^ (p2 - l.b()))) > 0;
} // proper same_side
fr bool is_intersect_segment(const line &l1, const line &l2) {
return (!same_side(l1, l2.a, l2.b())) && (!same_side(l2, l1.a, l1.b()));
} // proper intersect
fr double dist_p_line(const line &l, const point &p) {
return fabs(l.ab ^ (p - l.a)) / length(l.ab);
}
fr double dist_p_segment(const line &l, const point &p) {
if( dcmp(l.ab * (p - l.b())) >= 0 ) return dist(p, l.b());
else if( dcmp(l.ab * (p - l.a)) <= 0 ) return dist(p, l.a);
else return dist_p_line(l, p);
}
fr double dist_s_segment(const line &l1, const line &l2) {
if( is_intersect_segment(l1, l2) ) return 0;
return std::min(std::min(dist_p_segment(l2, l1.a), dist_p_segment(l2, l1.b())),
std::min(dist_p_segment(l1, l2.a), dist_p_segment(l1, l2.b())));
}
fr point intersect_line(const line &l1, const line &l2) {
return l1.a + l1.ab * ((l2.a - l1.a) ^ l2.ab) / (l1.ab ^ l2.ab);
}
fr point projection(const line &l, const point &p) {
return l.a + l.ab * ((p - l.a) * l.ab) / (l.ab * l.ab);
}
fr point reflection(const line & l, const point &p) {
return projection(l, p) * 2 - p;
}
};
void convex(point *a, int n, point *b, int &m) {
static point t[N + 5]; int k = n;
for(int i=0;i<k;i++) t[i] = a[i]; std::sort(t, t + k);
int o = 1; m = 0;
for(int i=0;i<k;i++) {
while( m > o && ((b[m - 1] - b[m - 2]) ^ (t[i] - b[m - 1])) >= 0 )
m--;
b[m++] = t[i];
}
o = m;
for(int i=k-2;i>=0;i--) {
while( m > o && ((b[m - 1] - b[m - 2]) ^ (t[i] - b[m - 1])) >= 0 )
m--;
b[m++] = t[i];
}
b[--m] = (point){0, 0};
}
/*
求凸包注意两个常常需要讨论的点:
(1)N 个点共线。
(2)凸包上的边是否允许多点共线(代码中给的是不允许共线的版本)。
*/
namespace half_plane{
line l[N + 5], ql[N + 5]; point qp[N + 5]; int s, t;
void insert(line x) {
while( s < t && !on_left(x, qp[t-1]) ) t--;
while( s < t && !on_left(x, qp[s]) ) s++;
//注意这里必须先弹队尾再弹队首
if( dcmp(x.ab ^ ql[t].ab) == 0 ) {
if( on_left(ql[t], x.a) ) ql[t] = x;
}
else ql[++t] = x;
if( s < t ) qp[t-1] = intersect_line(ql[t-1], ql[t]);
}
bool half_plane(line *A, int n) {
for(int i=0;i<n;i++) l[i] = A[i];
l[n++] = line((point){-INF, -INF}, (point){1, 0});
l[n++] = line((point){-INF, INF}, (point){0, -1});
l[n++] = line((point){INF, INF}, (point){-1, 0});
l[n++] = line((point){INF, -INF}, (point){0, 1});
std::sort(l, l + n); ql[s = t = 1] = l[0];
for(int i=1;i<n;i++) insert(l[i]);
while( s < t && !on_left(ql[s], qp[t-1]) ) t--; //模拟插入队首直线
return t - s + 1 >= 3;
}//只判断是否存在半平面交,其他用途同理
/*
用半平面交时注意一个细节:直线上的点是否在半平面上。
主要是最后半平面交可能退化成单点,如果按上面代码来写就是无解的
*/
}
int main() {
}