【模板】计算几何

平面几何模板

//平面几何模板
#include <bits/stdc++.h>
using namespace std;
const double eps=1e-10;
const double PI=acos(-1.0);

struct Point		//定义点or向量 
{
	double x,y;
	Point ( double x=0,double y=0 ) : x(x),y(y) {}
};
typedef Point Vector;		//区分点和向量，因为结构体一样 

Vector operator + ( Vector A,Vector B ) { return Vector( A.x+B.x,A.y+B.y ); } 
Vector operator - ( Point A,Point B ) { return Vector( A.x-B.x,A.y-B.y ); }
Vector operator * ( Vector A,double b ) { return Vector( A.x*b,A.y*b ); }
Vector operator / ( Vector A,double b ) { return Vector( A.x/b,A.y/b ); }
bool operator < ( const Point &a,const Point &b ) { return dcmp( a.x-b.x )<0 || ( dcmp(a.x-b.x)==0 && dcmp(a.y-b.y)<0 ); }
bool operator == ( const Point &a,const Point &b ) { return dcmp( a.x-b.x )==0 && dcmp( a.y-b.y )==0; }

int dcmp( double x )
{
	if ( fabs( x )<eps ) return 0;
	else return x<0 ? -1 : 1;
}
/*
bool operator < ( const Point &a,const Point &b )
{
	return a.x<b.x || ( a.x==b.x && a.y<b.y );
}*/

double angle( Vector v )		//求极角，就是向量旋转的度数 
{
	return atan2( v.y,v.x );
}

double Dot( Vector A,Vector B )		//求点积 
{
	return A.x*B.x+A.y*B.y;
}

double Length( Vector A )		//求向量的模 
{
	return sqrt( Dot( A,A ) );
}

double Angle( Vector A,Vector B )		//求夹角 
{
	return acos( Dot(A,B)/Length(A)/Length(B) );
}

double Cross( Vector A,Vector B )		//根据向量求叉积 
{
	return A.x*B.y-A.y*B.x;
}

double Area2( Point A,Point B,Point C )		//根据点求叉积 
{
	return Cross( B-A,C-A );
}

Vector Rotate( Vector A,double rad )		//向量旋转，rad是弧度
{
	return Vector( A.x*cos(rad)-A.y*sin(rad),A.x*sin(rad)+A.y*cos(rad) );
}

Vector Normal( Vector A )		//左转九十度长度归一，注意不能是零向量
{
	double l=Length( A );
	return Vector( -A.y/l,A.x/l );
}

Point Intersection( Point P,Vector v,Point Q,Vector w )		//直线P+tv，直线Q+tw的交点
{
	Vector u=P-Q; double t=Cross( w,u )/Cross( v,w );
	return P+v*t;
}

double Toline( Point P,Point A,Point B )		//点到直线的距离
{
	Vector v1=B-A,v2=P-A;
	return fabs( Cross(v1,v2) )/Length(v1);
}

double Tosegment( Point P,Point A,Point B )		//点到线段的距离
{
	if ( A==B ) return Length( P-A );
	Vector v1=B-A,v2=P-A,v3=P-B;
	if ( dcmp(Dot(v1,v2))<0 ) return Length(v2);
	else if ( dcmp(Dot(v1,v3))>0 ) return Length( v3 );
	else return fabs( Cross( v1,v2 )/Length(v1) );
} 

Point Projection( Point P,Point A,Point B )		//点在直线上的投影
{
	Vector v=B-A;
	return A+v*( Dot(v,P-A)/Dot(v,v) );
}

bool SegmentIntersect( Point a1,Point a2,Point b1,Point b2 )		//线段相交判定
{
	double c1=Cross( a2-a1,b1-a1 ),c2=Cross( a2-a1,b2-a1 );
	double c3=Cross( b2-b1,a1-b1 ),c4=Cross( b2-b1,a2-b1 );
	return dcmp(c1)*dcmp(c2)<0 && dcmp(c3)*dcmp(c4)<0;
}

bool Onsegment ( Point p,Point a1,Point a2 )		//判断点是否在一条线段上
{
	return dcmp( Cross(a1-p,a2-p) )==0 && dcmp( Dot(a1-p,a2-p) )<0;
}

typedef vector<Point> Polygon;

double PolygonArea( Polygon poly )		//多边形有向面积
{
	double area=0; int n=poly.size();
	for ( int i=1; i<n-1; i++ )
		area+=Cross( poly[i]-poly[0],poly[(i+1)%n]-poly[0] );
	return area/2;
}

Polygon CutPolygon( Polygon poly,Point A,Point B )		//用有向直线A->B切割多边形，返回左部 
{
	Polygon newpoly; int n=poly.size();
	for ( int i=0; i<n; i++ )
	{
		Point C=poly[i],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=Intersection( A,B-A,C,D-C );
			if ( Onsegment(ip,C,D) ) newpoly.push_back( ip );
		}
	}
	return newpoly;
}

struct Circle 
{
	Point c; double r;
	Circle ( Point c,double r ) : c(c),r(r) {}
	Point point( double a ) 
	{
		return Point( c.x+cos(a)*r,c.y+sin(a)*r ); 	//根据圆心角求坐标 
	}
};

struct Line
{
	Point p; Vector v;
	Line ( Point p,Vector v ) : p(p),v(v) {}
	Point point( double t ) 
	{
		return p+v*t;
	}
	Line move( double d ) 
	{
		return Line( p+Normal(v)*d,v );
	}
};

Point LineIntersection( Line a, Line b ) 		//line版本求直线交点 
{
  return Intersection( a.p,a.v,b.p,b.v );
}

void LineGeneralEquation( const Point &p1,const Point &p2,double &a,double &b,double &c )		//两点确定直线转化为一般式ax+by+c=0 
{
	a=p2.y-p1.y; b=p1.x-p2.x; c=-a*p1.x-b*p1.y;
}

int CircleIntersection_Line( Line L,Circle C,double &t1,double &t2,vector<Point> &sol )		//直线交圆 
{
	double a=L.v.x,b=L.p.x-C.c.x,c=L.v.y,d=L.p.y-C.c.y;
	double e=a*a+c*c,f=2*(a*b+c*d),g=b*b+d*d-C.r*C.r;
	double delta=f*f-4*e*g;
	if ( dcmp(delta)<0 ) return 0;
	if ( dcmp(delta)==0 )
	{
		t1=t2=-f/(2*e); sol.push_back( L.point(t1) ); return 1;
	} 
	t1=(-f-sqrt(delta))/(2*e); sol.push_back( L.point(t1) );
	t2=(-f+sqrt(delta))/(2*e); sol.push_back( L.point(t2) );
	return 2;
}

int CircleIntersection_Circle( Circle c1,Circle c2,vector<Point> &sol )	//两圆相交
{
	double d=Length( c1.c-c2.c );
	if ( dcmp(d)==0 ) 
		if ( dcmp(c1.r-c2.r)==0 ) return -1;
		else return 0;
	if ( dcmp(c1.r+c2.r-d)<0 ) return 0;
	if ( dcmp(fabs(c1.r-c2.r)-d)>0 ) return 0;
	double a=angle( c2.c-c1.c );
	double da=acos( (c1.r*c1.r+d*d-c2.r*c2.r)/(2*c1.r*d) );
	Point p1=c1.point( a-da ),p2=c1.point(a+da);
	sol.push_back(p1); 
	if ( p1==p2 ) return 1;
	sol.push_back(p2);
	return 2;
} 

int getTangents_Point( Point p,Circle C,Vector *v )		//过定点作切线
{
	Vector u=C.c-p; double dist=Length(u);
	if ( dist<C.r ) return 0;
	else if ( dcmp(dist-C.r)==0 ) 
			{
				v[0]=Rotate( u,PI/2 ); return 1;
			}
		else
		{
			double ang=asin( C.r/dist );
			v[0]=Rotate( u,-ang ); v[1]=Rotate( u,+ang );
			return 2;
		}
} 

int getTangents_Circle( Circle A,Circle B,Point *a,Point *b )		//两圆公切线
{
	int cnt=0;
	if ( A.r<B.r ) swap( A,B ),swap( a,b );
	int d2=( A.c.x-B.c.x )*( A.c.x-B.c.x )+( A.c.y-B.c.y )*( A.c.y-B.c.y );
	int rdiff=A.r-B.r,rsum=A.r+B.r;
	if ( d2<rdiff*rdiff ) return 0;
	double base=atan2( B.c.y-A.c.y,B.c.x-A.c.x );
	if ( d2==0 && A.r==B.r ) return -1;
	if ( d2==rdiff*rdiff )
	{
		a[cnt]=A.point(base); b[cnt]=B.point(base); cnt++; return 1;
	}
	double ang=acos( (A.r-B.r)/sqrt(d2) );
	a[cnt]=A.point(base+ang); b[cnt]=B.point(base+ang); cnt++;
	a[cnt]=A.point(base-ang); b[cnt]=B.point(base-ang); cnt++;
	if ( d2==rsum*rsum )
	{
		a[cnt]=A.point(base); b[cnt]=B.point(PI+base); cnt++;
	}
	else if ( d2>rsum*rsum )
	{
		double ang=acos( (A.r+B.r)/sqrt(d2) );
		a[cnt]=A.point(base+ang); b[cnt]=B.point(PI+base+ang); cnt++;
		a[cnt]=A.point(base-ang); b[cnt]=B.point(PI+base-ang); cnt++;
	}
	return cnt;
}

//____________________________________________________________________________________________________

Circle Circumcircle( Point p1,Point p2,Point p3 )		//三角形外接圆
{
	double bx=p2.x-p1.x,by=p2.y-p1.y;
	double cx=p3.x-p1.x,cy=p3.y-p1.y;
	double d=2*(bx*cy-by*cx);
	double xx=( cy*(bx*bx+by*by)-by*(cx*cx+cy*cy) )/d+p1.x;
	double yy=( bx*(cx*cx+cy*cy)-cx*(bx*bx+by*by) )/d+p1.y;
	Point p=Point( xx,yy );
	return Circle( p,Length(p1-p) );
}

Circle Inscribedcircle( Point p1,Point p2,Point p3 )		//内切圆
{
	double a=Length( p2-p3 ),b=Length( p3-p1 ),c=Length( p1-p2 );
	Point p=( p1*a+p2*b+p3*c )/(a+b+c);
	return Circle( p,Toline( p,p1,p2 ) );
}

//—————————————————————————————————— 
vector<Point> CircleThroughPoint( Point p,Line L,double r ) 	//过定点和直线相切的，半径为r的圆
{
	vector<Point> ans; double t1,t2;
	CircleIntersection_Line( L.move(-r),Circle(p,r),t1,t2,ans ); 
	CircleIntersection_Line( L.move(r),Circle(p,r),t1,t2,ans );		//圆心到直线半径为r的两条直线和到定点距离为r的圆的交点
	return ans; 
} 
//—————————————————————————————————— 

vector<Point> CircleTangent( Line a,Line b,double r )		//和给定两直线相切的，半径为r的圆
{
	vector<Point> ans;
	Line l1=a.move(-r),l2=a.move(r),l3=b.move(-r),l4=b.move(r);		//和每条直线距离为r的4条
	ans.push_back( LineIntersection(l1,l3) );
	ans.push_back( LineIntersection(l1,l4) ); 
	ans.push_back( LineIntersection(l2,l3) ); 
	ans.push_back( LineIntersection(l2,l4) ); 
	return ans;
}

vector<Point> TwoCircleExcircle( Circle c1,Circle c2,double r )		//和两个圆外切，半径为r的圆 
{
	vector<Point> ans; Vector v=c2.c-c1.c;
	double dis=Length(v); int d=dcmp( dis-c1.r-c2.r-r*2 );
	if ( d>0 ) return ans;
	CircleIntersection_Circle( Circle(c1.c,c1.r+r),Circle(c2.c,c2.r+r),ans );		//得到两个新的圆求交点 
}

//_______________________________________________________________________________________________________________
/*
int CircleIntersection_Line( Point A,Point B,Point C,double r,double &t1,double &t2 )		//直线交圆 
{
	double a=B.x-A.x,b=A.x-C.x,c=B.y-A.y,d=A.y-C.y;
	double e=a*a+c*c,f=2*(a*b+c*d),g=b*b+d*d-r*r;
	double delta=f*f-4*e*g;
	if ( dcmp(delta)<0 ) return 0;
	if ( dcmp(delta)==0 )
	{
		t1=t2=-f/(2*e); return 1;
	} 
	t1=(-f-sqrt(delta))/(2*e); t2=(-f+sqrt(delta))/(2*e);
	return 2;
}*/

bool CircleIntersectSegment( Point A,Point B,Point p,double R )		//圆和线段相交
{
	double t1,t2;
	int c=CircleIntersection_Line( A,B,p,R,t1,t2 );
	if ( c<=1 ) return 0;
	if ( dcmp(t1)>0 && dcmp(t1-1)<0 ) return 1;
	if ( dcmp(t2)>0 && dcmp(t2-1)<0 ) return 1;
	return 0;
} 

double Length2( Vector A )
{
	return Dot( A,A );
}

bool PointInCircle( Point p,Point center,double R )		//点在圆内，圆周不算
{
	return dcmp( Length2(p-center)-R*R )<0;
} 

double Degree( Vector v )			//格式化成角度版极角 
{
	double ang=angle( v )*180.0/PI;
	while ( dcmp(ang)<0 ) ang+=360.0;
	while ( dcmp(ang-180)>=0 ) ang-=180.0;
	return ang;
}

Line getLine( double x1,double y1,double x2,double y2 )		//点的坐标转成直线 
{
	Point p1( x1,y1 ),p2( x2,y2 );
	return Line( p1,p2-p1 );
}

double torad( double deg )				//角度转弧度 
{
	return deg/180*PI;
}

void get_coord( double R,double lat,double lng,double &x,double &y,double &z )		//经纬度（角度）转坐标 
{
	lat=torad( lat ); lng=torad( lng );
	x=R*cos(lat)*cos(lng); y=R*cos(lat)*sin(lng); z=R*sin(lat);
}

int PointInPolygon( const Point &p,const vector<Point> &poly )		//判定点在多边形内 
{
	int wn=0,n=poly.size();
	for ( int i=0; i<n; i++ )
	{
		const Point &p1=poly[i]; const Point &p2=poly[(i+1)%n];
		if ( p1==p || p2==p || Onsegment(p,p1,p2) ) return -1;		//边上 
		int d1=dcmp( p1.y-p.y ),d2=dcmp( p2.y-p.y );
		int k=dcmp( Cross( p2-p1,p-p1) );
		if ( k>0 && d1<=0 && d2>0 ) wn++;
		if ( k<0 && d2<=0 && d1>0 ) wn--; 
	}
	if ( wn!=0 ) return 1;		//在内部 
	return 0;
}

int ConvexHull( Point *p,int n,Point *ch )		//Andrew算法算凸包 
{
	sort( p,p+n ); int m=0;
	for ( int i=0; i<n; i++ )
	{
		while ( m>1 && Cross(ch[m-1]-ch[m-2],p[i]-ch[m-2])<=0 ) m--;
		ch[m++]=p[i];
	}
	int k=m;
	for ( int i=n-2; i>=0; i-- )
	{	
		while ( m>k && Cross(ch[m-1]-ch[m-2],p[i]-ch[m-2])<=0 ) m--;
		ch[m++]=p[i];
	}
	if ( n>1 ) m--;
	return m;
}

/*		
//凸包改版 
vector<Point> ConvexHull( vector<Point> p )		
{
	sort( p.begin(),p.end() ); p.erase( unique( p.begin(),p.end() ),p.end() );
	int n=p.size(); int m=0; vector<Point> ch(n+1);
	for ( int i=0; i<n; i++ )
	{
		while ( m>1 && Cross( ch[m-1]-ch[m-2],p[i]-ch[m-2] )<=0 ) m--;
		ch[m++]=p[i];
	}
	int k=m;
	for ( int i=n-2; i>=0; i-- )
	{	
		while ( m>k && Cross(ch[m-1]-ch[m-2],p[i]-ch[m-2])<=0 ) m--;
		ch[m++]=p[i];
	}
	if ( n>1 ) m--;  ch.resize( m );
	return ch;
}
*/

bool ConvexDisjoint( const vector<Point> ch1,const vector<Point> ch2 )		//判断两个凸包是否相交 
{
	int c1=ch1.size(),c2=ch2.size();
	for ( int i=0; i<c1; i++ )			//判断每个点不在另一个凸包内部 
		if ( PointInPolygon( ch1[i],ch2 )!=0 ) return 0;
	for ( int i=0; i<c2; i++ )
		if ( PointInPolygon( ch2[i],ch1 )!=0 ) return 0;
	for ( int i=0; i<c1; i++ )			//任意两条不同包的线段不相交 
	 for ( int j=0; j<c2; j++ )
	 	if ( SegmentIntersect( ch1[i],ch1[(i+1)%c1],ch2[j],ch2[(j+1)%c2] ) ) return 0;
	return 1;
}

struct Line2		//有向直线
{
	Point P; Vector v; double ang;
	Line2() {}
	Line2( Point P,Vector v ) :P(P),v(v) {ang=atan2( v.y,v.x );}
	bool operator <( const Line2 &L ) const
	{
		return ang<L.ang;
	}
};

bool Onleft( Line2 L,Point p )		//点p在有向直线左边 
{
	return Cross( L.v,p-L.P ) >0;
}

Point GetLineIntersection( Line2 a,Line2 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( Line2 *L,int n,Point *poly )		//半平面交 
{
	sort( L,L+n );
	int first,last; Point *p=new Point[n]; Line2 *q=new Line2[n];
	q[first=last=0]=L[0];
	for ( int i=1; i<n; i++ )
	{
		while ( first<last && !Onleft( L[i],p[last-1] ) ) last--;
		while ( first<last && !Onleft( L[i],p[first] ) ) first++;
		q[++last]=L[i];
		if ( fabs( Cross( q[last].v,q[last-1].v ) )<eps )
		{
			last--;
			if ( Onleft( q[last],L[i].P ) ) q[last]=L[i];
		}
		if ( first<last ) p[last-1]=GetLineIntersection( q[last-1],q[last] );
	}
	while (first<last && !Onleft( q[first],p[last-1] ) ) last--;
	if ( last-first<=1 ) return 0;
	p[last]=GetLineIntersection( q[last],q[first] );
	int m=0;
	for ( int i=first; i<=last; i++ )
		poly[m++]=p[i];
	return m;
}

/*
vector<Point> HalfPlaneIntersection( vector<Line2> L )		//半平面交改版 
{
	int n=L.size(); sort( L.begin(),L.end() );
	int first,last; vector<Point> p(n); vector<Line2> q(n); vector<Point> ans;
	q[first=last=0]=L[0];
	for ( int i=1; i<n; i++ )
	{
		while ( first<last && !Onleft( L[i],p[last-1] ) ) last--;
		while ( first<last && !Onleft( L[i],p[first] ) ) first++;
		q[++last]=L[i];
		if ( fabs( Cross( q[last].v,q[last-1].v ) )<eps )
		{
			last--;
			if ( Onleft( q[last],L[i].P ) ) q[last]=L[i];
		}
		if ( first<last ) p[last-1]=GetLineIntersection( q[last-1],q[last] );
	}
	while (first<last && !Onleft( q[first],p[last-1] ) ) last--;
	if ( last-first<=1 ) return ans;
	p[last]=GetLineIntersection( q[last],q[first] );
	for ( int i=first; i<=last; i++ )
		ans.push_back( p[i] );
	return ans;
}
*/

struct edge
{
	int from,to; double ang;
};
const int N=10010;

struct PSLG		//平面直线图PSLG 
{
	int n,m,face_cnt;
	double x[N],y[N],area[N];
	vector<edge> edges; vector<int> g[N]; vector<Polygon> faces;
	int vis[N*2],left[N*2],prev[N*2];		//每条边是否访问过，左边面编号，相同起点的上一条边编号 
	
	void init( int n )
	{
		this->n=n; 
		for ( int i=0; i<n; i++ )
			g[i].clear();
		edges.clear(); faces.clear();
	}
	
	double getAngle( int from,int to )	//极角
	{
		return atan2( y[to]-y[from],x[to]-x[from] );
	} 
	
	void addedge( int from,int to )
	{
		edges.push_back( (edge){from,to,getAngle(from,to)} );
		edges.push_back( (edge){to,from,getAngle(to,from)} );
		m=edges.size();
		g[from].push_back( m-2 ); g[to].push_back( m-1 );
	}
	
	void build()		//找faces，计算面积 
	{
		for ( int u=0; u<n; u++ )
		{
			int d=g[u].size();
			for ( int i=0; i<d; i++ )		//排序复杂度可以改进 
			 for ( int j=i+1; j<d; j++ )
			 	if ( edges[g[u][i]].ang>edges[g[u][j]].ang ) swap( g[u][i],g[u][j] );
			for ( int i=0; i<d; i++ )
				prev[g[u][(i+1)%d]]=g[u][i];
		}
		memset( vis,0,sizeof(vis) ); face_cnt=0;
		for ( int u=0; u<n; u++ )
		 for ( int i=0; i<g[u].size(); i++ )
		 {
		 	int e=g[u][i];
		 	if ( !vis[e] )
		 	{
		 		face_cnt++; Polygon poly;
		 		for ( ; ; )
		 		{
		 			vis[e]=1; left[e]=face_cnt;
		 			int from=edges[e].from;
		 			poly.push_back( Point(x[from],y[from]) );
		 			e=prev[e^1];
		 			if ( e==g[u][i] ) break;
		 			assert( vis[e]==0 );
				}
				faces.push_back( poly );
			}
		 }
		for ( int i=0; i<faces.size(); i++ )
			area[i]=PolygonArea( faces[i] );
	}
}g;

const int maxp=105;
int n,c;
Point P[maxp],V[maxp*(maxp-1)/2+maxp];

int getid( Point p )
{
	return lower_bound( V,V+c,p )-V;
}

Polygon simplify( const Polygon &poly )			//删除三点共线的情况
{
	Polygon ans; int n=poly.size();
	for ( int i=0; i<n; i++ )
	{
		Point a=poly[i],b=poly[(i+1)%n],c=poly[(i+2)%n];
		if ( dcmp(Cross(a-b,c-b))!=0 ) ans.push_back(b);
	}
	return ans;
}

void BuildGraph()		//输出平面直线图外轮廓 
{
	c=n; vector<double> dis[maxp];
	for ( int i=0; i<n; i++ )
		V[i]=P[i];
	for ( int i=0; i<n; i++ )
	 for ( int j=i+1; j<n; j++ )
	 	if ( SegmentIntersect( P[i],P[(i+1)%n],P[j],P[(j+1)%n] ) )
	 	{
	 		Point p=Intersection( P[i],P[(i+1)%n]-P[i],P[j],P[(j+1)%n]-P[j] );
	 		V[c++]=p; 
	 		dis[i].push_back( Length(p-P[i]) ); dis[j].push_back( Length(p-P[j]) );
		}
	sort( V,V+c ); c=unique( V,V+c )-V; 		//保证很接近的点看做同一个
	g.init( c ); 
	for ( int i=0; i<c; i++ )
		g.x[i]=V[i].x,g.y[i]=V[i].y;
	for ( int i=0; i<n; i++ )
	{
		Vector v=P[(i+1)%n]-P[i]; double len=Length(v);
		dis[i].push_back( 0 ); dis[i].push_back( len );
		sort( dis[i].begin(),dis[i].end() );
		int sz=dis[i].size();
		for ( int j=1; j<sz; j++ )
		{
			Point a=P[i]+v*(dis[i][j-1]/len); 
			Point b=P[i]+v*(dis[i][j]/len);
			if ( a==b ) continue;
			g.addedge( getid(a),getid(b) );
		}
	}
	g.build(); Polygon poly;
	for ( int i=0; i<g.faces.size(); i++ )
		if ( g.area[i]<0 )		//面积小于0的无限面 
		{
			poly=g.faces[i]; 
			reverse( poly.begin(),poly.end() ); 		//无限面对于内部来说，顶点顺时针
			poly=simplify( poly );
			break; 
		}
	int m=poly.size();
	printf( "%d\n",m );
	int start=0;
	for ( int i=0; i<m; i++ )
		if ( poly[i]<poly[start] ) start=i;
	for ( int i=start; i<m; i++ )
		printf( "%.4lf %.4lf\n",poly[i].x,poly[i].y );
	for ( int i=0; i<start; i++ )
		printf( "%.4lf %.4lf\n",poly[i].x,poly[i].y );
}

int main()
{
	
}

三维几何模板

/*
三维计算几何模板 by lyc  = = 
*/
#include <bits/stdc++.h>
using namespace std;
const double eps=1e-18;
struct Point3
{
	double x,y,z;
	Point3 (double x=0,double y=0,double z=0) : x(x),y(y),z(z) {}
};
typedef Point3 Vector3;

Vector3 operator + (Vector3 A,Vector3 B) { return Vector3(A.x+B.x,A.y+B.y,A.z+B.z); }
Vector3 operator - (Point3 A,Point3 B) { return Vector3(A.x-B.x,A.y-B.y,A.z-B.z); }
Vector3 operator * (Vector3 A,double p) { return Vector3(A.x*p,A.y*p,A.z*p); }
Vector3 operator / (Vector3 A,double p) { return Vector3(A.x/p,A.y/p,A.z/p); }
int dcmp( double x) { if ( fabs(x)<eps ) return 0; else return x<0 ? -1 : 1; }
bool operator == (Point3 p1,Point3 p2) 
{ return (dcmp(p1.x-p2.x)==0) && (dcmp(p1.y-p2.y)==0) && (dcmp(p1.z-p2.z)==0); }
double Dot(Vector3 A,Vector3 B) { return A.x*B.x+A.y*B.y+A.z*B.z; }		//三维点积 
double Length( Vector3 A ) { return sqrt( Dot(A,A) ); }
double Angle( Vector3 A,Vector3 B ) { return acos( Dot(A,B)/Length(A)/Length(B) ); }
//转化成弧度
double torad( double deg ) { return deg/180.0*PI; }
//点p到平面p0-n的距离（n为平面法向量，垂直于平面） 
double DistanceToPlane (Point3 p,Point3 p0,Vector3 n) 
{ return fabs(Dot(p-p0,n))/Length(n); }
//p在平面p0-n的投影 感觉书上是错的
Point3 GetPlaneProjection (Point3 p,Point3 p0,Vector3 n) { return p-n*Dot(p-p0,n); }
//直线p1-p2到平面p0-n的交点（假定唯一） 
Point3 LinePlaneIntersection (Point3 p1,Point3 p2,Point3 p0,Vector3 n)
{ Vector3 v=p2-p1; double t=( Dot(n,p0-p1)/Dot(n,p2-p1) ); return p1+v*t; }
//三维叉积
Vector3 Cross (Vector3 A,Vector3 B)
{ return Vector3( A.y*B.z-A.z*B.y,A.z*B.x-A.x*B.z,A.x*B.y-A.y*B.x ); } 
double Area2( Point3 A,Point3 B,Point3 C )
{ return Length( Cross(B-A,C-A) ); }
//点p在△p0p1p2
bool PointInTri ( Point3 P,Point3 p0,Point3 p1,Point3 p2 )
{ 
	double area1=Area2(P,p0,p1),area2=Area2(P,p1,p2),area3=Area2(P,p2,p0);
  	return dcmp(area1+area2+area3-Area2(p0,p1,p2))==0; 
} 
//△p0p1p2是否和线段AB相交
bool TriSegIntersection( Point3 p0,Point3 p1,Point3 p2,Point3 A,Point3 B,Point3 &p )
{
	Vector3 n=Cross( p1-p0,p2-p0 );
	if ( dcmp(Dot(n,B-A))==0 ) return 0; 		//平行或者共面
	else
	{
		double t=Dot(n,p0-A)/Dot(n,B-A);
		if ( dcmp(t)<0 || dcmp(t-1)>0 ) return 0; 		//交点不在线段上
		p=A+(B-A)*t; return PointInTri(p,p0,p1,p2);		//判断交点是否在三角形内 
	} 
} 
//点P到直线AB的距离
double DistanceToLine ( Point3 P,Point3 A,Point3 B ) 
{
	Vector3 v1=B-A,v2=P-A; return Length( Cross(v1,v2)/Length(v1) );
}
//点p到线段AB的距离
double DistanceToSegment( Point3 P,Point3 A,Point3 B )
{
	if ( A==B ) return Length(P-A);
	Vector3 v1=B-A,v2=P-A,v3=P-B;
	if ( dcmp(Dot(v1,v2))<0 ) return Length(v2);
	else if ( dcmp(Dot(v1,v3))>0 ) return Length(v3);
	else return Length(Cross(v1,v2))/Length(v1);
}
//线段到线段最短距离
double MinSegmentToSegmentDis( Point3 p1,Point3 p2,Point3 p3,Point3 p4) 
{
    Vector3 v=Cross(p2-p1,p4-p3);
    double A=Dot(Cross(p2-p1,v),p3-p1);
    double B=Dot(Cross(p2-p1,v),p4-p1);
    double C=Dot(Cross(p3-p4,v),p1-p4);
    double D=Dot(Cross(p3-p4,v),p2-p4);
    if ( dcmp(A*B)<=0 && dcmp(C*D)<=0 ) return fabs(Dot(p1-p3,v))/Length(v);
    return min(min(DistanceToSegment(p1,p3,p4),DistanceToSegment(p2,p3,p4)),min(DistanceToSegment(p3,p1,p2),DistanceToSegment(p4,p1,p2)));   
}
//AB，AC，AD的混合积，就是四面体ABCD有向体积6倍
double Volume6( Point3 A,Point3 B,Point3 C,Point3 D )
{
	return Dot( D-A,Cross(B-A,C-A) );
}

struct Face
{
	int v[3];
	Face( int a,int b,int c ) { v[0]=a; v[1]=b; v[2]=c; }
	Vector3 normal( Point3 *P ) const { return Cross( P[v[1]]-P[v[0]],P[v[2]]-P[v[0]] ); }
	int cansee( Point3 *P,int i ) const { return Dot( P[i]-P[v[0]],normal(P) ) >0 ? 1 : 0; }
};
//增量法求三维凸包，不考虑特殊情况（如四点共面），实践中要对输入的点进行微小扰动
vector<Face> CH3D( Point3 *P,int n )
{
	vector<Face> cur; cur.push_back( (Face){{0,1,2}} ); cur.push_back( (Face){{2,1,0}} );
	for ( int i=3; i<n; i++ )
	{
		vector<Face> nxt;
		for ( int j=0; j<cur.size(); j++ )
		{
			Face &f=cur[j]; int res=f.cansee(P,i);
			if ( !res ) nxt.push_back(f);
			for ( int k=0; k<3; k++ )
				vis[f.v[k]][f.v[(k+1)%3]]=res;
		}
		for ( int j=0; j<cur.size(); j++ )
		 for ( int k=0; k<3; k++ )
		 {
		 	int a=cur[j].v[k],b=cur[j].v[(k+1)%3];
		 	if ( vis[a][b]!=vis[b][a] && vis[a][b] ) nxt.push_back( (Face){{a,b,i}} );
		 	//(a,b)是分界线，左边对p[i]可见 
		 }
		cur=nxt;
	}
	return cur;
} 
//扰动
double rand01() { return rand()/(double)RAND_MAX; }		//0-1之间的 
double randeps() { return (rand01()-0.5)*eps; } 		//-eps/2~eps/2
Point3 add_noise( Point3 p ) { return Point3(p.x+randeps(),p.y+randeps(),p.z+randeps() ); } 
//求异面直线p1+su和p2+tv的公垂线对应的s，如果平行/重合为false
bool LineDistance3D( Point3 p1,Vector3 u,Point3 p2,Vector3 v,double &s )
{
	double b=Dot(u,u)*Dot(v,v)-Dot(u,v)*Dot(u,v);
	if ( dcmp(b)==0 ) return 0;
	double a=Dot(u,v)*Dot(v,p1-p2)-Dot(v,v)*Dot(u,p1-p2);
	s=a/b; return 1;
} 
//三棱锥重心：各个坐标平均值 
Point3 Centroid( const Point3& A,const Point3& B,const Point3& C,const Point3& D ) 
{ 
	return (A+B+C+D)/4.0; 
}
//多面体 
struct ConvexPolyhedron 
{
    int n;
    vector<Point3> P, P2;
    vector<Face> faces;
    //读入 
    bool read() 
	{
        if ( scanf("%d", &n)!=1 ) return 0;
        P.resize(n); P2.resize(n);
        for ( int i=0; i<n; i++ ) 
		{ 
			P[i]=read_point3(); P2[i]=add_noise(P[i]); 
		}
        faces=CH3D(P2); return 1;
    }
    //重心：按面把多面体分成若干个三棱锥，求重心，重心坐标加权（质量）为多面体重心。 
    Point3 centroid() 
	{
        Point3 C=P[0],tot(0,0,0); double totv=0;
        for ( int i=0; i<faces.size(); i++ ) 
		{
            Point3 p1=P[faces[i].v[0]],p2=P[faces[i].v[1]],p3=P[faces[i].v[2]];
            double v=-Volume6(p1, p2, p3, C);
            totv+=v; tot=tot+Centroid(p1, p2, p3, C)*v;
        }
        return tot/totv;
    }
    //各个面到点C的最短距离 
    double mindist( Point3 C ) 
	{
        double ans=1e30;
        for( int i=0; i<faces.size(); i++ ) 
		{
            Point3 p1=P[faces[i].v[0]],p2=P[faces[i].v[1]],p3=P[faces[i].v[2]];
            ans=min(ans,fabs(-Volume6(p1,p2,p3,C)/Area2(p1, p2, p3)));
        }
        return ans;
    }
};

int main()
{
}