三维凸包

复杂度 \(O(n^2)\)，可以处理多点共线、共面的情况，所有点共线返回 \(0\)，所有点共面返回一个退化成二维的凸包（两个面都有）

手动处理的共线共面，扰动点的话精度稍微要求高一些就寄了

这里有比较牛逼的数据：https://darkbzoj.cc/problem/1209

inline int sign(const double &a){return a>EPS?1:(a<-EPS?-1:0);}
struct Vector{
	double x,y,z;
	inline Vector(){x=y=z=0;}
	inline Vector(double _x,double _y,double _z){x=_x;y=_y;z=_z;}
	inline double len(){return __builtin_sqrt(x*x+y*y+z*z);}
	inline double len2(){return x*x+y*y+z*z;}
	inline Vector operator + (const Vector &a)const{return Vector(x+a.x,y+a.y,z+a.z);}
	inline Vector operator - (const Vector &a)const{return Vector(x-a.x,y-a.y,z-a.z);}
	inline Vector operator * (const double &a)const{return Vector(x*a,y*a,z*a);}
	inline Vector operator / (const double &a)const{return Vector(x/a,y/a,z/a);}
	inline double operator * (const Vector &a)const{return x*a.x+y*a.y+z*a.z;}
	inline Vector operator ^ (const Vector &a)const{return Vector(y*a.z-z*a.y,z*a.x-x*a.z,x*a.y-y*a.x);}
};
struct Plane{
	Vector a,b,c;//a->b->c counterclockwise
	int x,y,z;
	inline Plane(){}
	inline Plane(const Vector &_a,const Vector &_b,const Vector &_c,int _x=0,int _y=0,int _z=0){a=_a;b=_b;c=_c;x=_x;y=_y;z=_z;}
	inline Vector normal()const{return (c-b)^(a-b);}
	inline double area()const{return ((c-b)^(a-b)).len()/2;}
};
inline int isAbove(const Plane &p,const Vector &v){
	return sign(p.normal()*(v-p.a))==1;
}
inline int isCollinear(const Vector &a,const Vector &b,const Vector &c){
	return sign(((b-a)^(c-a)).len2())==0;
}
inline int isCoplanar(const Vector &a,const Vector &b,const Vector &c,const Vector &d){
	return isCollinear(a,b,c)||sign(Plane(a,b,c).normal()*(d-a))==0;
}
#define N 2006
inline int work(int n,Vector *a,Plane *p){
	if(n<3) return 0;
	int s[4];
	s[0]=1;s[1]=2;s[2]=3;s[3]=4;
	while(s[2]<=n&&isCollinear(a[s[0]],a[s[1]],a[s[2]])) s[2]++;
	if(s[2]>n) return 0;
	while(s[3]<=n&&isCoplanar(a[s[0]],a[s[1]],a[s[2]],a[s[3]])) s[3]++;
	if(s[3]>n) s[3]=n;
	if(!isAbove(Plane(a[s[0]],a[s[1]],a[s[2]]),a[s[3]])) std::swap(s[2],s[3]);
	int tot=0;
	p[++tot]=Plane(a[s[0]],a[s[2]],a[s[1]],s[0],s[2],s[1]);
	p[++tot]=Plane(a[s[0]],a[s[1]],a[s[3]],s[0],s[1],s[3]);
	p[++tot]=Plane(a[s[0]],a[s[3]],a[s[2]],s[0],s[3],s[2]);
	p[++tot]=Plane(a[s[1]],a[s[2]],a[s[3]],s[1],s[2],s[3]);
	static Plane u[N*2],v[N*2];
	static int tag[N][N];
	for(int i=1;i<=n;i++){
		int _u=0,_v=0;
		for(int j=1;j<=tot;j++){
			if(!isAbove(p[j],a[i])) u[++_u]=p[j];
			else{
				v[++_v]=p[j];
				tag[p[j].x][p[j].y]=tag[p[j].y][p[j].z]=tag[p[j].z][p[j].x]=i;
			}
		}
		tot=0;
		for(int j=1;j<=_v;j++){
			if(tag[v[j].y][v[j].x]!=i) p[++tot]=Plane(v[j].a,v[j].b,a[i],v[j].x,v[j].y,i);
			if(tag[v[j].z][v[j].y]!=i) p[++tot]=Plane(v[j].b,v[j].c,a[i],v[j].y,v[j].z,i);
			if(tag[v[j].x][v[j].z]!=i) p[++tot]=Plane(v[j].c,v[j].a,a[i],v[j].z,v[j].x,i);
		}
		for(int j=1;j<=_u;j++) p[++tot]=u[j];
	}
	return tot;
}