「计算几何」模板

基础模板。

#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() {
	
}
posted @ 2020-01-01 18:29  Tiw_Air_OAO  阅读(98)  评论(0编辑  收藏
把博客园图标替换成自己的图标
把博客园图标替换成自己的图标end