半平面交的模板

网上看到个 把半平面交转化为求凸包的方法

也即是把一条直线y=kx+b看成一个点(k,b)求凸包

详细证明及其代码见：http://trinkle.blog.uoj.ac/blog/235

#include <bits/stdc++.h>
using namespace std;
#define REP(i,a,b)  for(int i=a;i<=(int)(b);++i)
#define REPD(i,a,b) for(int i=a;i>=(int)(b);--i)
const double PI=acos(-1);
const double EPS=1e-6;
int dcmp(double x){
	if(fabs(x)<EPS) return 0;
	return x > 0 ? 1 : -1;
}
struct Point{
	double x,y;
};
typedef Point Vector;
bool operator < (const Point& p1, const Point& p2){	return p1.x<p2.x || (p1.x==p2.x && p1.y<p2.y) ;}
Vector operator / (const Point& A, double x){	return Vector{A.x/x, A.y/x};}
Vector operator * (const Vector& A, double x){	return Vector{A.x*x, A.y*x};}
Vector operator - (const Vector& A, const Vector& B){	return Vector{A.x-B.x,A.y-B.y};}
Vector operator + (const Point& A, const Vector& v){	return Point{A.x+v.x,A.y+v.y};}
Vector Rotate(const Point& p,double ang){	return Vector{p.x*cos(ang)-p.y*sin(ang),p.x*sin(ang)+p.y*cos(ang)};}
double Cross(const Vector& A, const Vector& B){	return A.x*B.y-A.y*B.x;}
double Dot(Vector A, Vector B){	return A.x*B.x+A.y*B.y;}
double Length(Vector v){	return sqrt(fabs(Dot(v,v)));}
Vector Unit(Vector v){	return v/Length(v);}

Point GetLineIntersection(Point P, Vector v, Point Q, Vector w){
	Vector u=P-Q;
	double t=Cross(w,u) / Cross(v,w);
	return P+v*t;
}
bool Onsegment(Point p, Point A, Point B){
	return dcmp(Cross(A-p, B-p)) == 0 && dcmp(Dot(A-p, B-p)) < 0 ;
}
struct Line
{
	Point p;
	Vector v;
	double ang;
	Line(){}
	Line(Point p, Vector v):p(p),v(v){ ang=atan2(v.y,v.x); }
	bool operator < (const Line& l) const {
		return ang < l.ang;
	}
};
typedef vector<Point> Polygon;
Polygon CutPolygon(Polygon poly, Point A, Point B){  //A-B切多边形poly
	Polygon newpoly;
	int n=poly.size();
	for(int i=0;i<n;i++){
		Point C=poly[i];
		Point D=poly[(i+1)%n];
		if(dcmp(Cross(B-A, C-A))>0) newpoly.push_back(C);
		if(dcmp(Cross(B-A, C-D))!=0){   //保证两线不共线
			Point ip=GetLineIntersection(A,B-A,C,D-C);
			if(Onsegment(ip,C,D)) newpoly.push_back(ip);
		}
	}
	return newpoly;
}

bool OnLeft(Point p, Line L){
	return dcmp(Cross(L.v, p-L.p)) > 0;
}

Point GetLineIntersection(Line a, Line b){
	Vector u=a.p-b.p;
	double t=Cross(b.v, u) / Cross(a.v, b.v);
	return a.p+a.v*t;
}
int HalfplaneIntersection(Line* L, int n, Point* poly){  //L组成的半平面交，前提是如果是无界的，那么首先要加个"包围框"
	sort(L, L+n);

	int first,last;
	Line* Q=new Line[n];
	Point* p=new Point[n];
	Q[first=last=0]=L[0];
	for(int i=1;i<n;i++) {
		while(first<last && !OnLeft(p[last-1], L[i])) last--;   //最后一个点不在L[i]的左侧
		while(first<last && !OnLeft(p[first], L[i])) first++;   //第一个一个点不在L[i]的左侧
		Q[++last]=L[i];
		if(dcmp(Cross(Q[last].v, Q[last-1].v)) == 0){    //最后两条线平行
			last--;
			if(OnLeft(L[i].p, Q[last])) Q[last]=L[i];   //新加入的线L[i]在最后一条线的左侧，更新最后一条线
		}
		if(first<last) p[last-1]=GetLineIntersection(Q[last-1], Q[last]);
	}

	while(first<last && !OnLeft(p[last-1], Q[first])) last--;  //确保最后一个点在第一套线的左侧

	if(last-first<=1) return 0;    //形成不了封闭的图形  
	p[last]=GetLineIntersection(Q[first], Q[last]);    //第一条线和最后一条线求交点

	int m=0;
	for(int i=first;i<=last;i++) poly[m++]=p[i];    //把答案添入poly
	return m;
}