三维凸包 模板

只会写增量法 orz

例题：P2287
随机种子 0x383494 被卡了精度，eps=1e-10 太大了

#include <cstdio>
#include <iostream>
#include <bitset>
#include <list>
#include <random>
#include <cmath>
#include <algorithm>
#define UP(i,s,e) for(auto i=s; i!=e; ++i)
std::mt19937 rdna(0x383494);
double constexpr eps = 1e-11;  // 警钟敲烂
double reps(){ static std::uniform_real_distribution<> r(-1.0, 1.0); return r(rdna)*eps; }
constexpr int N = 2000;
namespace calc_geo{ // }{{{
struct Vec{
	double x, y, z;
	Vec(){}
	Vec(double x, double y, double z):x(x), y(y), z(z){}
	//Vec &operator=(Vec b){ x=b.x, y=b.y, z=b.z; return *this; }
	Vec operator+(Vec b){ b.x+=x, b.y+=y, b.z+=z; return b; }
	Vec operator-(Vec b){ b.x=x-b.x, b.y=y-b.y, b.z=z-b.z; return b; }
	Vec operator-(){ return Vec(-x, -y, -z); }
	double operator^(Vec b){ return x*b.x+y*b.y+z*b.z; }
	Vec operator%(Vec b){ 
		return Vec(y * b.z - z * b.y,
				z * b.x - x * b.z,
				x * b.y - y * b.x); 
	}
	double len(){ return sqrt(x*x+y*y+z*z); }
	Vec &normal(){ double l = len(); if(l){x/=l, y/=l, z/=l;} return *this; }
	void shake(){ x += reps(); y += reps(); z += reps(); }
} vecs[N];
std::bitset<N> vis[N];
struct Surface{
	Vec *vertex[3]; // outside: anticlockwise
	Vec normal(){ // point towards outside
		return (*vertex[1]-*vertex[0])%(*vertex[2]-*vertex[1]);
	}
	double area(){
		return ((*vertex[1]-*vertex[0])%(*vertex[2]-*vertex[1])).len()/2.0;
	}
};
struct Convex{
	std::list<Surface> surs;
	void init(Vec *A, Vec *B, Vec *C){
		surs.clear();
		surs.push_back({A,B,C});
		surs.push_back({C,B,A});
	}
	void addpoint(Vec *A){
		UP(BCD, surs.begin(), surs.end()){
			if((BCD->normal()^(*A-*BCD->vertex[0])) > 0){
				UP(j, 0, 3){
					vis[BCD->vertex[j]-vecs][BCD->vertex[(j+1)%3]-vecs]
						= 1;
				}
				BCD = surs.erase(BCD);
				--BCD;
				continue;
			}
		}
		auto ee = surs.end();
		UP(BCD, surs.begin(), ee){
			UP(j, 0, 3){
				int idx, idy;
				idx = (int)(BCD->vertex[j]-vecs);
				idy = (int)(BCD->vertex[(j+1)%3]-vecs);
				if(!vis[idx][idy] && vis[idy][idx]){
					surs.push_back({vecs+idy, vecs+idx, A});
					vis[idy][idx] = 0;
				}
			}
		}
	}
	double getarea(){
		double ans = 0;
		for(auto BCD:surs){
			ans += BCD.area();
		}
		return ans;
	}
} con;
} // {}}}
using std::cin;
using std::cout;
namespace m{ // }{{{
using calc_geo::vecs;
using calc_geo::con;
int in;
void work(){
	cin >> in;
	UP(i, 0, in) {
		cin >> vecs[i].x >> vecs[i].y >> vecs[i].z;
		vecs[i].shake();
	}
	con.init(vecs, vecs+1, vecs+2);
	UP(i, 3, in){
		con.addpoint(vecs+i);
	}
	char outbuf[20];
	sprintf(outbuf, "%.6lf\n", con.getarea());
	cout << outbuf;
}
} // {}}}
int main(){m::work(); return 0;}