# @总结 - 14@ 计算几何相关模板

#include <cmath>
#include <cstdio>
#include <algorithm>
using namespace std;

const double PI = acos(-1.0);
const double INF = 1E9;
const double EPS = 1E-9;
const int MAXN = 100000;

int dcmp(double x) {return (fabs(x) <= EPS ? 0 : (x > 0 ? 1 : -1));}

struct point{
double x, y;
point() : x(0), y(0) {}
point(double _x, double _y) : x(_x), y(_y) {}

friend point operator - (point A, point B) {return point(A.x - B.x, A.y - B.y);}
friend point operator + (point A, point B) {return point(A.x + B.x, A.y + B.y);}
friend point operator * (point A, double k) {return point(A.x * k, A.y * k);}
friend point operator / (point A, double k) {return point(A.x / k, A.y / k);}
friend double operator * (point A, point B) {return A.x * B.x + A.y * B.y;} //点积
friend double operator ^ (point A, point B) {return A.x * B.y - A.y * B.x;} //叉积
friend bool operator == (point A, point B) {return dcmp(A.x - B.x) == 0 && dcmp(A.y - B.y) == 0;}
friend bool operator < (point A, point B) {return (A.x == B.x ? A.y < B.y : A.x < B.x);}

friend double length(point A) {return sqrt(A * A);}
friend double angle(point P) {
double a = atan2(P.y, P.x);
return (a < 0 ? a + 2*PI : a);
}
friend double angle(point P, point Q) {
double a = atan2(P ^ Q, P * Q);
a = (a < 0 ? a + 2*PI : a);
return (a > PI ? 2*PI - a : a);
}
friend double area(point pnt[], int N) {
double sum = 0;
for(int i=0;i<N;i++)
sum += (pnt[i] ^ pnt[(i + 1) % N]);
return fabs(sum) / 2;
}
friend double dist(point A, point B) {
return length(A - B);
}

friend point normal(point A) {return point(A.y, -A.x);} //法向量
friend point unit(point A) {return A / length(A);}
friend point rotate(point A, double t) {
return point(cos(t)*A.x - sin(t)*A.y, cos(t)*A.y + sin(t)*A.x);
}
friend point bisector(point A, point B) { //角平分线
return A * length(B) + B * length(A);
}

friend double slope(point a, point b) {
if( a.x == b.x )
return a.y < b.y ? INF : -INF;
else return (a.y - b.y) / (a.x - b.x);
}
friend void read(point &p) {scanf("%lf%lf", &p.x, &p.y);}
};

struct line{
point a, ab; //直线 a + ab*t
double ang;
line() : a(), ab(), ang() {}
line(point _a, point _ab) : a(_a), ab(_ab), ang(atan2(_ab.y, _ab.x)) {}
point b() {return a + ab;}
line(double _a, double _b, double _c) {
ab = point(-_b, _a);
a = (_a == 0 ? point(0, -_c/_b) : point(-_c/_a, 0));
}//直线标准式 Ax + By + C = 0 的转化

friend bool on_left(line L, point P) {
return dcmp(L.ab ^ (P - L.a)) > 0;
}
friend bool operator < (line a, line b) {return a.ang < b.ang;}
friend bool is_on_line(line L, point P) {
return dcmp(((P - L.a)^(P - L.b()))) == 0;
}
friend bool is_on_ray(line L, point P) {
return is_on_line(L, P) && (dcmp((P - L.a)*L.ab) >= 0 );
}
friend bool is_on_segment(line L, point P) {
return is_on_ray(L, P) && (dcmp((P - L.b())*L.ab) <= 0);
}
friend bool is_same_side(line L, point P1, point P2) {
return (dcmp(((P1 - L.a)^(P1 - L.b()))*((P2 - L.a)^(P2 - L.b()))) > 0);
}
friend bool is_intersect_ray(line L1, line L2) {
return dcmp(((L2.a - L1.a)^L2.ab)/(L1.ab^L2.ab)) > 0 &&
dcmp(((L1.a - L2.a)^L1.ab)/(L2.ab^L1.ab)) > 0;
}
friend bool is_intersect_segment(line L1, line L2) {
return (!is_same_side(L1, L2.a, L2.b())) && (!is_same_side(L2, L1.a, L1.b()));
}
friend double dist_line(line L, point P) {
return fabs(L.ab ^ (P - L.a)) / length(L.ab);
}
friend double dist_ray(line L, point P) {
if( dcmp(L.ab * (P - L.a)) <= 0 ) return dist(P, L.a);
else return dist_line(L, P);
}
friend double dist_segment(line L, point P) {
if( dcmp(L.ab * (P - L.b())) >= 0 ) return dist(P, L.b());
else return dist_ray(L, P);
}
friend double dist_ray(line L1, line L2) {
if( is_intersect_ray(L1, L2) ) return 0;
return min(dist_ray(L2, L1.a), dist_ray(L1, L2.a));
}
friend double dist_segment(line L1, line L2) {
if( is_intersect_segment(L1, L2) ) return 0;
return min(min(dist_segment(L2, L1.a), dist_segment(L2, L1.b())),
min(dist_segment(L1, L2.a), dist_segment(L1, L2.b())));
}
friend point intersect_line(line L1, line L2) {
return L1.a + L1.ab*((L2.a - L1.a)^L2.ab)/(L1.ab^L2.ab);
}
friend point projection(line L, point P) {
return L.a + L.ab*((P - L.a) * L.ab)/(L.ab * L.ab);
}//投影
friend point reflection(line L, point P) {
return projection(L, P)*2 - P;
}//对称点
friend void transform(line L, double &A, double &B, double &C) {
A = L.ab.y, B = -L.ab.x, C = L.a.y*L.ab.x - L.a.x*L.ab.y;
}//直线标准式 Ax + By + C = 0 的转化
l = line(P1, P2 - P1);
}
};

void convex(point *A, int n, point *B, int &m) {
static point t[MAXN + 5], s[MAXN + 5];
for(int i=0;i<n;i++) t[i] = A[i]; sort(t, t + n);

int cnt = 0, tp = 0;
for(int i=0;i<n;i++) {
if( i && A[i] == A[i-1] ) continue; // 注意判重点
while( tp >= 2 && slope(s[tp - 1], s[tp]) > slope(s[tp], t[i]) )
tp--;
s[++tp] = t[i];
}
for(int i=1;i<=tp;i++) B[cnt++] = s[i];
tp = 0;
for(int i=0;i<n;i++) {
if( i && A[i] == A[i-1] ) continue; // 注意判重点
while( tp >= 2 && slope(s[tp - 1], s[tp]) < slope(s[tp], t[i]) )
tp--;
s[++tp] = t[i];
}
for(int i=tp-1;i>=2;i--) B[cnt++] = s[i];
m = cnt;
}
/*

（1）N 个点共线。
（2）凸包上的边是否允许多点共线（代码中给的是允许共线的版本）。

update in 2020/06/03：写凸包可以不用求斜率，直接判叉积。
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]; 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};
}
*/

line l[MAXN + 5], ql[MAXN + 5]; point qp[MAXN + 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));

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;
}//只判断是否存在半平面交，其他用途同理
/*

*/

posted @ 2020-01-01 18:29  Tiw_Air_OAO  阅读(66)  评论(0编辑  收藏