hdu 4273 Rescue（三维凸包 三维中心）

http://acm.hdu.edu.cn/showproblem.php?pid=4273

倍增法求三维凸包，求三维凸包中心。。。各种模板

#include<cstdio>
#include<algorithm>
#include<vector>
#include<cmath>
#include<queue>
#pragma comment(linker, "/STACK:1024000000,1024000000")
using namespace std;
// http://acm.hdu.edu.cn/showproblem.php?pid=4273
const int MAXN = 105;
const double eps = 1e-8;
inline int cmp(double a)
{
	return (a> eps) - (a<-eps);
}
struct _Point
{
	double x, y, z;
	_Point(double a, double b, double c)
	{x = a, y = b, z = c;}
	_Point() {}
	_Point operator - (const _Point & a)
	{
		return _Point(x-a.x, y-a.y, z-a.z);
	}
	_Point operator + (const _Point & a)
	{
		return _Point(x+a.x, y+a.y, z+a.z);
	}
	double operator ^ (const _Point & a)
	{
		return x*a.x + y*a.y + z*a.z;
	}
	_Point operator * (const _Point & a)
	{
		return _Point(y*a.z-z*a.y, z*a.x-x*a.z, x*a.y-y*a.x);
	}
	_Point operator * (const double a)
	{
		return _Point(x*a, y*a, z*a);
	}
	_Point operator / (const double a)
	{
		return _Point(x/a, y/a, z/a);
	}
	double dis()
	{
		return sqrt(x*x + y*y + z*z);
	}
};
struct _face
{
	int a, b, c;
	int ok;
	_face(){}
	_face(int x, int y, int z, int s): a(x),b(y),c(z),ok(s) {}
};
int np, nf;
_Point p[MAXN];
_face f[MAXN*8];
int g[MAXN][MAXN];
double fpcmp(_face & ff, _Point & pp)
{
	_Point b = p[ff.b] - p[ff.a];
	_Point c = p[ff.c] - p[ff.a];
	_Point a = pp - p[ff.a];
	return (b*c)^a;
}
void dfs(int pp, int cur);
void deal(int pp, int a, int b)
{
	int df = g[a][b];
	_face add;
	if (f[df].ok)
	{
		if (fpcmp(f[df], p[pp]) > eps)
			dfs(pp, df);
		else
		{
			add.a = b;
			add.b = a;
			add.c = pp;
			add.ok = 1;
			g[b][a] = g[a][pp] = g[pp][b] = nf;
			f[nf++] = add;
		}
	}
}
void dfs(int pp, int cur)
{
	f[cur].ok = 0;
	deal(pp, f[cur].b, f[cur].a);
	deal(pp, f[cur].c, f[cur].b);
	deal(pp, f[cur].a, f[cur].c);
}
double getVolume(_Point o, _Point a, _Point b, _Point c)
{
	return (a-o)*(b-o)^(c-o);
}
_Point getCenter()
{
	_Point res(0,0,0), o;
	o = res;
	double all = 0;
	for (int i = 0; i< nf; ++i)
	{
		double val = getVolume(o, p[f[i].a],p[f[i].b],p[f[i].c]);
		res = res + (o+p[f[i].a]+p[f[i].b]+p[f[i].c])/4.0*val;
		all += val;
	}
	res = res/all;
	return res;
}
void build()
{
	for (int i = 2; i< np; ++i)
	{
		if (((p[0]-p[1])*(p[1]-p[i])).dis() > eps)
		{
			swap(p[2], p[i]);
			break;
		}
	}
	for (int i = 3; i< np; ++i)
	{
		if (cmp( ((p[0]-p[1])*(p[0]-p[2])) ^ (p[0]-p[i]) ))
		{
			swap(p[3], p[i]);
			break;
		}
	}
	for (int i = 0; i< 4; ++i)
	{
		_face add((i+1)%4,(i+2)%4,(i+3)%4, 1);
		if (fpcmp(add, p[i]) > 0) swap(add.b, add.c);
		g[add.a][add.b] = g[add.b][add.c] = g[add.c][add.a] = nf;
		f[nf++] = add;
	}
	for (int i = 4; i< np; ++i)
	{
		for (int j = 0; j< nf; ++j)
		{
			if (f[j].ok && fpcmp(f[j], p[i]) > eps)
			{
				dfs(i, j);
				break;
			}
		}
	}
	int tp = nf;
	nf = 0;
	for (int i=0; i< tp; ++i)
	{
		if (f[i].ok) f[nf++] = f[i];
	}
}
double getArea(_Point & a, _Point & b, _Point & c)
{
	return ((b-a)*(c-a)).dis();
}
double dis_pf(_Point pp, int ff)
{
	return getVolume(p[f[ff].a], p[f[ff].b], p[f[ff].c], pp)
			/getArea(p[f[ff].a], p[f[ff].b], p[f[ff].c]);
}
int main()
{
#ifndef ONLINE_JUDGE
    freopen("in.txt", "r", stdin);
#endif
	while (scanf("%d", &np) !=EOF)
	{
		for (int i = 0; i< np; ++i)
			scanf("%lf%lf%lf", &p[i].x, &p[i].y, &p[i].z);
		// 前四点不共面;
		nf = 0;
		build();
		_Point o = getCenter();
		double mn = 1e100;
		for (int i = 0; i< nf; ++i)
		{
			mn = min(mn, fabs(dis_pf(o, i)));
		}
		printf("%.3lf\n", mn);
	}
	return 0;
}