jiedai计算几何模板

计算几何模板

基本函数与常量

const double eps=1e-10;
const double PI=acos(-1);
inline int sgn(double x){return fabs(x)<eps?0:(x<0?-1:1);}

二维点类

struct point{
	double x,y;
	point(double a=0,double b=0):x(a),y(b){}
	point operator +(const point &A)const{
		return point(x+A.x,y+A.y);
	}
	point operator -(const point &A)const{
		return point(x-A.x,y-A.y);
	}
	point operator *(const double v)const{
		return point(x*v,y*v);
	}
	point operator /(const double v)const{
		return point(x/v,y/v);
	}
	bool operator ==(const point &A)const{
		return sgn(x-A.x)==0&&sgn(y-A.y)==0;
	}
	double norm(){
		return sqrt(x*x+y*y);
	}
};

点积

double dot(point a,point b){
	return a.x*b.x+a.y*b.y;
}

叉积

double det(point a,point b){
	return a.x*b.y-a.y*b.x;
}

两点之间距离

double dist(point a,point b){
	return (a-b).norm();
}

极角

注意：atan2函数比较慢

double angle(point p){
	return atan2(p.y,p.x);
}

点绕原点逆时针旋转（弧度）

点p绕原点O逆时针旋转弧度A：rotate_point(p,A)

点a绕点b逆时针旋转弧度c：rotate_point(a-b,c)+b

point rotate_point(point p,double A){
	return point(p.x*cos(A)-p.y*sin(A),p.x*sin(A)+p.y*cos(A));
}

二维线段直线类

struct line{
	point a,b;
	line(){}
	line(point x,point y):a(x),b(y){}
};

点到直线的距离

double point_to_line(point p,point s,point t){
	return fabs(det(s-p,t-p)/dist(s,t));
}

点到线段的距离

double point_to_segment(point p,point s,point t){
	if(sgn(dot(p-s,t-s))==-1)return dist(p,s);
	if(sgn(dot(p-t,s-t))==-1)return dist(p,t);
	return fabs(det(s-p,t-p)/dist(s,t));
}

点在直线上的垂足

point point_proj_line(point p,point s,point t){
	double r=dot(t-s,p-s)/dot(t-s,t-s);
	return s+(t-s)*r;
}

判断点在直线上

bool point_on_line(point p,point s,point t){
	return sgn(det(p-s,p-t))==0;
}

判断点在线段上

包含端点，不含端点则将<=改成<

bool point_on_segment(point p,point s,point t){
	return sgn(det(p-s,p-t))==0&&sgn(dot(p-s,p-t))<=0;
}

判断直线平行

bool parallel(line a,line b){
	return sgn(det(a.a-a.b,b.a-b.b))==0;
}

两条直线的交点

如果有交点则返回true，且交点保存在res中

需要注意的是，两直线重合也返回false

bool line_make_point(line a,line b,point &res){
	if(parallel(a,b))return 0;
	double s1=det(a.a-b.a,b.b-b.a);
	double s2=det(a.b-b.a,b.b-b.a);
	res=(a.b*s1-a.a*s2)/(s1-s2);
	return 1;
}

若两直线不平行，以下代码直接返回交点

point line_make_point(line a,line b){
	double s1=det(a.a-b.a,b.b-b.a);
	double s2=det(a.b-b.a,b.b-b.a);
	return (a.b*s1-a.a*s2)/(s1-s2);
}

三角形的有向面积（需*0.5）

double cross(point o,point a,point b){//*0.5
	return (a.x-o.x)*(b.y-o.y)-(a.y-o.y)*(b.x-o.x);
}

判断直线和线段是否相交

判断直线a和线段b是否相交

1表示规范相交（只有一个交点且线段端点不在直线上）

2表示不规范相交（线段的一个或两个端点在直线上）

0表示不相交

int line_and_segment(line a,line b){
	int d1=sgn(cross(a.a,a.b,b.a));
	int d2=sgn(cross(a.a,a.b,b.b));
	if((d1^d2)==-2)return 1;
	if(d1==0||d2==0)return 2;
	return 0;
}

判断线段和线段是否相交

bool segment_and_segment(line a,line b){
	if(point_on_line(a.a,b.a,b.b)&&point_on_line(a.b,b.a,b.b))//四点共线
		return point_on_segment(a.a,b.a,b.b)||point_on_segment(a.b,b.a,b.b)||
			point_on_segment(b.a,a.a,a.b)||point_on_segment(b.b,a.a,a.b);
	return sgn(det(a.b-a.a,b.a-a.a)*det(a.b-a.a,b.b-a.a))<=0&&
		sgn(det(b.b-b.a,a.a-b.a)*det(b.b-b.a,a.b-b.a))<=0;
}

将直线沿法向量方向平移

沿a.a->a.b的左侧平移距离len

line move_d(line a,double len){
	point e=a.b-a.a;
	e=rotate_point(e/e.norm(),PI/2);
	return line(a.a+e*len,a.b+e*len);
}

多边形类

多边形顶点坐标按逆时针排序

struct polygon{
    static const int SIZE=1005;
	int n;
	point A[SIZE];
	polygon(){}
};

计算多边形面积

	double area(){
		double res=0;
		A[n+1]=A[1];
		for(int i=1;i<=n;i++)res+=det(A[i],A[i+1]);
		return res/2.0;
	}

计算多边形周长

	double perimeter(){
		double res=0;
		A[n+1]=A[1];
		for(int i=1;i<=n;i++)res+=(A[i]-A[i+1]).norm();
		return res;
	}

判断点是否在多边形内部

	int point_in(point p){
		int res=0;
		A[n+1]=A[1];
		for(int i=1;i<=n;i++){
			if(point_on_segment(p,A[i],A[i+1]))return 2;
			int k=sgn(det(A[i+1]-A[i],p-A[i]));
			int d1=sgn(A[i].y-p.y);
			int d2=sgn(A[i+1].y-p.y);
			if(k>0&&d1>0&&d2<=0)res++;
			if(k<0&&d2>0&&d1<=0)res--;
		}
		return res!=0;
	}

计算多边形的重心

	point mass_center(){
		point res=point(0,0);
		if(sgn(area())==0)return res;
		A[n+1]=A[1];
		for(int i=1;i<=n;i++)res=res+(A[i]+A[i+1])*det(A[i],A[i+1]);
		return res/area()/6.0;
	}

多边形边界上的格点个数

	int gcd(int a,int b){
		if(b==0)return a;
		return gcd(b,a%b);
	}
	int border_int_point(){
		int res=0;
		A[n+1]=A[1];
		for(int i=1;i<=n;i++)res+=gcd(abs(int(A[i].x-A[i+1].x)),abs(int(A[i].y-A[i+1].y)));
		return res;
	}

多边形内的格点个数

Pick定理：对于顶点均为整点的简单多边形，其内部（不含边界）的格点个数为a，其边界上的格点个数为b，其面积为S

它们满足关系式：$S=a+\frac{b}{2}-1$

变形可得：$a=S-\frac{b}{2}+1$

皮克定理与欧拉公式（$V-E+F=2$）等价

	int inside_int_point(){
		return int(area())-border_int_point()/2+1;
	}

凸多边形类

struct polygon_convex{
	vector<point>p;
	polygon_convex(int sz=0){
		p.resize(sz);
	}
};

计算凸包面积

	double area(){
		double res=0;
		for(int i=0;i<p.size();i++)res+=det(p[i],p[(i+1)%p.size()]);
		return res/2.0;
	}

计算凸包周长

	double perimeter(){
		double res=0;
		for(int i=0;i<p.size();i++)res+=(p[i]-p[(i+1)%p.size()]).norm();
		return res;
	}

求点集的凸包（逆时针顺序）

polygon_convex convex_hull(vector<point>vi){
	polygon_convex res(2*vi.size()+5);
	sort(vi.begin(),vi.end(),[&](point &a,point &b){
		if(sgn(a.x-b.x)!=0)return a.x<b.x;
		return a.y<b.y;
	});
	vi.erase(unique(vi.begin(),vi.end()),vi.end());
	int m=0;
	for(int i=0;i<vi.size();i++){
		while(m>1&&sgn(det(res.p[m-1]-res.p[m-2],vi[i]-res.p[m-2]))<=0)m--;
		res.p[m++]=vi[i];
	}
	int k=m;
	for(int i=int(vi.size())-2;i>=0;i--){
		while(m>k&&sgn(det(res.p[m-1]-res.p[m-2],vi[i]-res.p[m-2]))<=0)m--;
		res.p[m++]=vi[i];
	}
	res.p.resize(m);
	if(vi.size()>1)res.p.resize(m-1);
	return res;
}

判断点是否在凸包内部

需要凸包逆时针

写法一：$O(n)$，true表示点在凸包内部或边界上

bool contain_On(const polygon_convex &a,const point &b){
	int n=a.p.size(),flag=0;
	for(int i=0;i<n;i++){
		int x=sgn(det(a.p[i]-b,a.p[(i+1)%n]-b));
		if(x){
			if(flag==0)flag=x;
			else if(flag!=x)return 0;
		}
	}
	return 1;
}

写法二：$O(lo{g}_{2}n)$

半平面交

P4196

以下代码能求出半平面交所形成的闭合凸包，向量的左侧为半平面。如果闭合凸包不存在，或者半平面无交集，需要另行判断。

const double eps=1e-10;
const double PI=acos(-1);
inline int sgn(double x){return fabs(x)<eps?0:(x<0?-1:1);}
struct point{
	double x,y;
	point(double a=0,double b=0):x(a),y(b){}
	point operator -(const point &A)const{
		return point(x-A.x,y-A.y);
	}
	point operator *(const double v)const{
		return point(x*v,y*v);
	}
	point operator /(const double v)const{
		return point(x/v,y/v);
	}
};
struct line{
	point a,b;
	double ang;
	line(){}
	line(point x,point y):a(x),b(y){}
};
double det(point a,point b){
	return a.x*b.y-a.y*b.x;
}
point line_make_point(line a,line b){
	double s1=det(a.a-b.a,b.b-b.a);
	double s2=det(a.b-b.a,b.b-b.a);
	return (a.b*s1-a.a*s2)/(s1-s2);
}
struct polygon_convex{
	vector<point>p;
	polygon_convex(int sz=0){
		p.resize(sz);
	}
	double area(){
		double res=0;
		for(int i=0;i<p.size();i++)res+=det(p[i],p[(i+1)%p.size()]);
		return res/2.0;
	}
};
bool on_right(line l1,line l2,line l3){
	point o=line_make_point(l2,l3);
	return sgn(det(l1.b-l1.a,o-l1.a))<0;
}
double angle(line l){
	return atan2(l.b.y-l.a.y,l.b.x-l.a.x);
}
void halfplane_intersection(vector<line>vi,polygon_convex &res){
	res.p.clear();
	for(int i=0;i<vi.size();i++)vi[i].ang=angle(vi[i]);
	sort(vi.begin(),vi.end(),[&](line l1,line l2){
		if(sgn(l1.ang-l2.ang)==0)return det(l1.b-l1.a,l2.b-l1.a)<0;
		return l1.ang<l2.ang;
	});
	static const int SIZE=1e6+5;
	static line Q[SIZE];
	int L=0,R=0;
	for(int i=0;i<vi.size();i++){
		if(i>0&&sgn(angle(vi[i])-angle(vi[i-1]))==0)continue;
		while(L+1<R&&on_right(vi[i],Q[R-1],Q[R]))R--;
		while(L+1<R&&on_right(vi[i],Q[L+1],Q[L+2]))L++;
		Q[++R]=vi[i];
	}
	while(L+1<R&&on_right(Q[L+1],Q[R-1],Q[R]))R--;
	while(L+1<R&&on_right(Q[R],Q[L+1],Q[L+2]))L++;
	Q[R+1]=Q[L+1];
	for(int i=L+1;i<=R;i++)res.p.push_back(line_make_point(Q[i],Q[i+1]));
}

平面最远点对（直线距离）

double convex_diameter(polygon_convex &a,int &first,int &second){
	vector<point>&vi=a.p;
	int n=vi.size();
	if(n<=1){
		first=second=0;
		return 0.0;
	}
	double res=0;
	for(int i=0,j=0;i<n;i++){
		while((j+1)%n!=i&&sgn(det(vi[(i+1)%n]-vi[i],vi[j]-vi[i])-det(vi[(i+1)%n]-vi[i],vi[(j+1)%n]-vi[i]))<=0)j=(j+1)%n;
		double d=dist(vi[i],vi[j]);
		if(d>res)res=d,first=i,second=j;
		d=dist(vi[(i+1)%n],vi[j]);
		if(d>res)res=d,first=(i+1)%n,second=j;
	}
	return res;
}

三维点类

struct point3{
	double x,y,z;
	point3(double a=0,double b=0,double c=0):x(a),y(b),z(c){}
	point3 operator +(const point3 &A)const{
		return point3(x+A.x,y+A.y,z+A.z);
	}
	point3 operator -(const point3 &A)const{
		return point3(x-A.x,y-A.y,z-A.z);
	}
	point3 operator *(const double v)const{
		return point3(x*v,y*v,z*v);
	}
	point3 operator /(const double v)const{
		return point3(x/v,y/v,z/v);
	}
	bool operator ==(const point3 &A)const{
		return sgn(x-A.x)==0&&sgn(y-A.y)==0&&sgn(z-A.z)==0;
	}
};

向量长度

	double length(){
		return sqrt(x*x+y*y+z*z);
	}

单位向量

	point3 unit(){
		return *this/length();
	}

点积

double dot(point3 a,point3 b){
	return a.x*b.x+a.y*b.y+a.z*b.z;
}

叉积

point3 det(point3 a,point3 b){
	return point3(a.y*b.z-a.z*b.y,a.z*b.x-a.x*b.z,a.x*b.y-a.y*b.x);
}

混合积

除以6就是a,b,c三个向量所构成的四面体的体积

double mix(point3 a,point3 b,point3 c){
	return dot(a,det(b,c));
}

两点之间距离

double dist(point3 a,point3 b){
	return (a-b).length();
}

三维直线与平面类

struct line3{
	point3 a,b;
	line3(){}
	line3(point3 x,point3 y):a(x),b(y){}
	double length(){
		return (a-b).length();
	}
};

struct plane3{
	point3 a,b,c;
	plane3(){}
	plane3(point3 x,point3 y,point3 z):a(x),b(y),c(z){}
};

平面法向量

point3 normal(point3 a,point3 b,point3 c){
	return det(a-b,b-c);
}

point3 normal(plane3 s){
	return det(s.a-s.b,s.b-s.c);
}

判断三点一线

bool points_one_line(point3 a,point3 b,point3 c){
	return sgn((det(a-b,b-c)).length())==0;
}

判断四点共面

bool points_one_plane(point3 a,point3 b,point3 c,point3 d){
	return sgn(dot(normal(a,b,c),d-a))==0;
}

判断直线平行

bool parallel(line3 l1,line3 l2){
	return sgn(det(l1.a-l1.b,l2.a-l2.b).length())==0;
}

判断直线垂直

bool perpendicular(line3 l1,line3 l2){
	return sgn(dot(l1.a-l1.b,l2.a-l2.b))==0;
}

判断平面内两点在直线同侧

bool same_side(point3 a,point3 b,line3 l){
	return dot(det(l.a-l.b,a-l.b),det(l.a-l.b,b-l.b))>eps;
}

判断平面内两点在直线异侧

bool opposite_side(point3 a,point3 b,line3 l){
	return dot(det(l.a-l.b,a-l.b),det(l.a-l.b,b-l.b))<-eps;
}

判断点在线段上（包含端点）

bool point_on_segment_in(point3 p,line3 l){
	return sgn(det(p-l.a,p-l.b).length())==0&&
		(l.a.x-p.x)*(l.b.x-p.x)<eps&&
		(l.a.y-p.y)*(l.b.y-p.y)<eps&&
		(l.a.z-p.z)*(l.b.z-p.z)<eps;
}

判断点在线段上（不含端点）

bool point_on_segment_ex(point3 p,line3 l){
	return point_on_segment_in(p,l)&&
		(sgn(p.x-l.a.x)||sgn(p.y-l.a.y)||sgn(p.z-l.a.z))&&
		(sgn(p.x-l.b.x)||sgn(p.y-l.b.y)||sgn(p.z-l.b.z));
}

判断线段是否相交（包含端点）

bool segment_and_segment_in(line3 l1,line3 l2){
	if(points_one_plane(l1.a,l1.b,l2.a,l2.b)==0)return 0;
	if(points_one_line(l1.a,l2.a,l2.b)&&points_one_line(l1.b,l2.a,l2.b))
		return point_on_segment_in(l1.a,l2)||point_on_segment_in(l1.b,l2)||
			point_on_segment_in(l2.a,l1)||point_on_segment_in(l2.b,l1);
	return same_side(l1.a,l1.b,l2)==0&&same_side(l2.a,l2.b,l1)==0;
}

判断线段是否相交（不含端点）

注意：当线段共线时输出0

bool segment_and_segment_ex(line3 l1,line3 l2){
	return points_one_plane(l1.a,l1.b,l2.a,l2.b)&&opposite_side(l1.a,l1.b,l2)&&opposite_side(l2.a,l2.b,l1);
}

两条直线的交点（需保证相交）

point3 line_make_point(line3 l1,line3 l2){
	double t=((l1.a.x-l2.a.x)*(l2.a.y-l2.b.y)-(l1.a.y-l2.a.y)*(l2.a.x-l2.b.x))
		/((l1.a.x-l1.b.x)*(l2.a.y-l2.b.y)-(l1.a.y-l1.b.y)*(l2.a.x-l2.b.x));
	return l1.a+(l1.b-l1.a)*t;
}

点到直线的距离

double point_to_line(point3 p,line3 l){
	return det(p-l.a,l.b-l.a).length()/l.length();
}

直线到直线的距离（不平行）

double line_to_line(line3 l1,line3 l2){
	point3 n=det(l1.a-l1.b,l2.a-l2.b);
	return fabs(dot(l1.a-l2.a,n))/n.length();
}

两条直线的夹角的cos值

double angle_cos(line3 l1,line3 l2){
	return dot(l1.a-l1.b,l2.a-l2.b)/(l1.a-l1.b).length()/(l2.a-l2.b).length();
}

点绕向量逆时针旋转（弧度）

点p绕向量Ol逆时针旋转弧度A：rotate_point(p,l,A)

点p绕直线l逆时针旋转弧度A：rotate_point(p-l.a,l.b-l.a,A)+l.a

point3 rotate_point(point3 p,point3 l,double A){
	l=l.unit();
	return p*cos(A)+det(l,p)*sin(A)+l*dot(l,p)*(1-cos(A));
}

判断两点在平面同侧

bool same_side(point3 a,point3 b,plane3 s){
	return dot(normal(s),a-s.a)*dot(normal(s),b-s.a)>eps;
}

判断两点在平面异侧

bool opposite_side(point3 a,point3 b,plane3 s){
	return dot(normal(s),a-s.a)*dot(normal(s),b-s.a)<-eps;
}

判断直线与平面平行

bool parallel(line3 l,plane3 s){
	return sgn(dot(l.a-l.b,normal(s)))==0;
}

判断两平面平行

bool parallel(plane3 s1,plane3 s2){
	return sgn(det(normal(s1),normal(s2)).length())==0;
}

判断直线与平面垂直

bool perpendicular(line3 l,plane3 s){
	return sgn(det(l.a-l.b,normal(s)).length())==0;
}

判断两平面垂直

bool perpendicular(plane3 s1,plane3 s2){
	return sgn(dot(normal(s1),normal(s2)))==0;
}

判断点在三角形内（包含边界）

bool point_in_plane_in(point3 p,plane3 s){
	return sgn(det(s.a-s.b,s.a-s.c).length()
		-det(p-s.a,p-s.b).length()
		-det(p-s.b,p-s.c).length()
		-det(p-s.c,p-s.a).length())==0;
}

判断点在三角形内（不含边界）

bool point_in_plane_ex(point3 p,plane3 s){
	return point_in_plane_in(p,s)&&
		det(p-s.a,p-s.b).length()>eps&&
		det(p-s.b,p-s.c).length()>eps&&
		det(p-s.c,p-s.a).length()>eps;
}

判断线段与三角形相交（包含边界）

bool segment_plane_in(line3 l,plane3 s){
	return same_side(l.a,l.b,s)==0&&
		same_side(s.a,s.b,plane3(l.a,l.b,s.c))==0&&
		same_side(s.b,s.c,plane3(l.a,l.b,s.a))==0&&
		same_side(s.c,s.a,plane3(l.a,l.b,s.b))==0;
}

判断线段与三角形相交（不含边界）

bool segment_plane_ex(line3 l,plane3 s){
	return opposite_side(l.a,l.b,s)&&
		opposite_side(s.a,s.b,plane3(l.a,l.b,s.c))&&
		opposite_side(s.b,s.c,plane3(l.a,l.b,s.a))&&
		opposite_side(s.c,s.a,plane3(l.a,l.b,s.b));
}

直线与平面的交点

point3 line_plane_make_point(line3 l,plane3 s){
	point3 n=normal(s);
	double t=(n.x*(s.a.x-l.a.x)+n.y*(s.a.y-l.a.y)+n.z*(s.a.z-l.a.z))
		/(n.x*(l.b.x-l.a.x)+n.y*(l.b.y-l.a.y)+n.z*(l.b.z-l.a.z));
	return l.a+(l.b-l.a)*t;
}

两平面的交线

bool plane_make_line(plane3 s1,plane3 s2,line3 &res){
	if(parallel(s1,s2))return 0;
	if(parallel(line3(s2.a,s2.b),s1))res.a=line_plane_make_point(line3(s2.b,s2.c),s1);
	else res.a=line_plane_make_point(line3(s2.a,s2.b),s1);
	res.b=res.a+det(normal(s1),normal(s2));
	return 1;
}

line3 plane_make_line(plane3 s1,plane3 s2){
	line3 res;
	if(parallel(line3(s2.a,s2.b),s1))res.a=line_plane_make_point(line3(s2.b,s2.c),s1);
	else res.a=line_plane_make_point(line3(s2.a,s2.b),s1);
	res.b=res.a+det(normal(s1),normal(s2));
	return res;
}

点到面的距离

double point_to_plane(point3 p,plane3 s){
	return fabs(dot(normal(s),p-s.a))/normal(s).length();
}

两平面的夹角的cos值

需要注意，区分正负

double angle_cos(plane3 s1,plane3 s2){
	return dot(normal(s1),normal(s2))/normal(s1).length()/normal(s2).length();
}

直线与平面的夹角的sin值

需要注意，区分正负

double angle_sin(line3 l,plane3 s){
	return dot(l.a-l.b,normal(s))/(l.a-l.b).length()/normal(s).length();
}

圆类

struct circle{
	point p;
	double r;
	circle(){}
	circle(point p,double r):p(p),r(r){}
};

inline double mysqrt(double x){return sqrt(max(0.0,x));}

圆与直线的交点

inline double mysqrt(double x){return sqrt(max(0.0,x));}
vector<point>circle_line_make_point(point a,point b,point o,double r){
	double dx=b.x-a.x,dy=b.y-a.y;
	double A=dx*dx+dy*dy;
	double B=2*dx*(a.x-o.x)+2*dy*(a.y-o.y);
	double C=(a.x-o.x)*(a.x-o.x)+(a.y-o.y)*(a.y-o.y)-r*r;
	double delta=B*B-4*A*C;
	vector<point>vi;
	if(sgn(delta)>=0){
		double t1=(-B+mysqrt(delta))/(2*A);
		double t2=(-B-mysqrt(delta))/(2*A);
		vi.push_back(point(a.x+t1*dx,a.y+t1*dy));
		if(sgn(delta)>0)vi.push_back(point(a.x+t2*dx,a.y+t2*dy));
	}
	return vi;
}

最小圆覆盖

平面上n个点，求一个半径最小的圆，能覆盖所有的点

时间复杂度期望为 $O(n)$

struct point{
	double x,y;
	point(double a=0,double b=0):x(a),y(b){}
	point operator +(const point &A)const{
		return point(x+A.x,y+A.y);
	}
	point operator -(const point &A)const{
		return point(x-A.x,y-A.y);
	}
	point operator /(const double v)const{
		return point(x/v,y/v);
	}
	double norm2(){
		return x*x+y*y;
	}
};
point circle_center(point a,point b,point c){
	double a1=b.x-a.x,b1=b.y-a.y,c1=(a1*a1+b1*b1)/2;
	double a2=c.x-a.x,b2=c.y-a.y,c2=(a2*a2+b2*b2)/2;
	double d=a1*b2-a2*b1;
	return point(a.x+(c1*b2-c2*b1)/d,a.y+(a1*c2-a2*c1)/d);
}
void min_circle_cover(int n,point A[],point &o,double &r){
    srand(time(NULL));
	random_shuffle(A+1,A+n+1);
	o=point(0,0),r=0.0;
	for(int i=1;i<=n;i++)if((A[i]-o).norm2()>r){
		o=A[i],r=0;
		for(int j=1;j<i;j++)if((A[j]-o).norm2()>r){
			o=(A[i]+A[j])/2,r=(A[j]-o).norm2();
			for(int k=1;k<j;k++)if((A[k]-o).norm2()>r)
				o=circle_center(A[i],A[j],A[k]),r=(A[k]-o).norm2();
		}
	}
	r=sqrt(r);
}

其他知识或技巧

三角形重心、外心、垂心、内心

point triangle_mass_center(point a,point b,point c){//重心
	return (a+b+c)/3;
}
point triangle_circle_center(point a,point b,point c){//外心
	double a1=b.x-a.x,b1=b.y-a.y,c1=(a1*a1+b1*b1)/2;
	double a2=c.x-a.x,b2=c.y-a.y,c2=(a2*a2+b2*b2)/2;
	double d=a1*b2-a2*b1;
	return point(a.x+(c1*b2-c2*b1)/d,a.y+(a1*c2-a2*c1)/d);
}
point triangle_orthocenter(point a,point b,point c){//垂心
	return triangle_mass_center(a,b,c)*3-triangle_circle_center(a,b,c)*2;
}
point triangle_innercenter(point a,point b,point c){//内心
	double la=(b-c).norm();
	double lb=(c-a).norm();
	double lc=(a-b).norm();
	return point((la*a.x+lb*b.x+lc*c.x)/(la+lb+lc),(la*a.y+lb*b.y+lc*c.y)/(la+lb+lc));
}

四面体体积

给定四面体六条棱的长度，求四面体的体积

已知四面体棱 $ab,ac,ad,bc,bd,cd$ 的长度分别为 $l,n,a,m,b,c$

double tetrahedral_volume(double l,double n,double a,double m,double b,double c){
	double x,y;
	x=4*a*a*b*b*c*c-a*a*(b*b+c*c-m*m)*(b*b+c*c-m*m)-b*b*(c*c+a*a-n*n)*(c*c+a*a-n*n);
	y=c*c*(a*a+b*b-l*l)*(a*a+b*b-l*l)-(a*a+b*b-l*l)*(b*b+c*c-m*m)*(c*c+a*a-n*n);
	return sqrt(x-y)/12;
}