二维

点和向量化为坐标 Coord 进行运算，使用 stl 中的 complex 实现。

复数相乘的几何意义为长度相乘，极角相加。

用直线上的一点 p 和方向向量 v 表示一条经过 p 的直线，直线上的所有点 q 满足 q=p+t*v，其中 t 是参数；当限制 t≥0 时，该参数方程表示射线；限制 0≤t≤1 时，该参数方程表示线段。

此外，如果已知线段端点 a1 和 a2，可以通过 Line(a1,a2-a1)来得到对应的参数形式。

Morley 定理：三角形每个内角的三等分线相交成等边三角形。

欧拉定理：平面图的点数 V、边数 E 和面数 F 满足 V+F-E=2。

typedef double lf;
typedef complex<lf> Coord;
const lf EPS = 1e-9;
#define X real()
#define Y imag()
struct Line
{
	Coord p, v;
	Line(Coord p = Coord(), Coord v = Coord()) : p(p), v(v) {}
	Coord point(lf t) const { return p + v * t; }
};
struct Circle
{
	Coord c;
	lf r;
	Circle(Coord c = Coord(), lf r = 0) : c(c), r(r) {}
	Coord point(lf t) const { return c + polar(r, t); } //t为参数，幅角
};
/*
Coord(lf x=0,lf y=0);//构造函数
lf real(Coord a);//a的实部（复平面的横坐标）,也可写作a.real()
lf imag(Coord a);//a的虚部（复平面的纵坐标）,也可写作a.imag()

lf abs(Coord a);//向量a的模长，或是点a到原点的距离
lf norm(Coord a);//abs的平方，比abs快，但是要注意浮点数精度溢出

lf arg(Coord a);//a的幅角，与atan2(a.real(),a.imag())等价
Coord polar(lf r,lf t);//极坐标生成方式，r为幅值，t为幅角

//运算符重载+、-、*、/（以及对应的赋值运算，但是赋值运算不能写在表达式中，详见参考地址）、<<、>>（输出括号形式的坐标）
*/
int sgn(lf d) { return (d > EPS) - (d < -EPS); }

bool operator!=(const Coord &A, const Coord &B) { return sgn(A.X - B.X) || sgn(A.Y - B.Y); } //不等运算符，涉及到浮点数比较要重写

bool operator==(const Coord &A, const Coord &B) { return !(A != B); }

bool cmpCoord(const Coord &A, const Coord &B) { return sgn(A.X - B.X) ? sgn(A.X - B.X) < 0 : sgn(A.Y - B.Y) < 0; } //复数没有小于运算，只能这样定义一个比较函数

bool cmpLine(const Line &A, const Line &B) { return sgn(arg(A.v) - arg(B.v)) < 0; } //按极角排序，求凸包中使用

lf Dot(const Coord &A, const Coord &B) { return A.X * B.X + A.Y * B.Y; }

lf Cross(const Coord &A, const Coord &B) { return A.X * B.Y - B.X * A.Y; }

lf Angle(const Coord &A, const Coord &B) { return acos(Dot(A, B) / abs(A) / abs(B)); }

lf Area2(const Coord &A, const Coord &B, const Coord &C) { return Cross(B - A, C - A); } //三角形ABC有向面积的两倍

Coord Rotate(const Coord &A, lf rad) { return A * polar(1.0, rad);} //向量A逆时针旋转rad弧度

Coord Normal(const Coord &A) //A的法向量，把A逆时针旋转九十度并长度化为1
{
	lf L = abs(A);
	return Coord(-A.Y / L, A.X / L);
}

bool onLeft(const Coord &P, const Line &L) { return sgn(Cross(L.v, P - L.p)) > 0; } //p是否在有向直线L左侧，不含线上

lf DistanceToLine(const Coord &P, const Line &L) { return Cross(L.v, P - L.p) / abs(L.v); } //点到直线距离（有向）

lf DistanceToLine(const Coord &P, const Coord &A, const Coord &B) { return DistanceToLine(P, Line(A, B - A)); }

lf DistanceToSegment(const Coord &P, const Coord &A, const Coord &B) //点到线段的距离（无向）
{
	if (A == B)
		return abs(P - A);
	Coord v1 = B - A, v2 = P - A, v3 = P - B;
	if (sgn(Dot(v1, v2)) < 0)
		return abs(v2);
	if (sgn(Dot(v1, v3)) > 0)
		return abs(v3);
	return fabs(DistanceToLine(P, Line(A, B - A)));
}

Coord getLineProjection(const Coord &P, const Line &L) { return L.point(Dot(L.v, P - L.p) / norm(L.v)); } //点在直线上的投影

Coord getLineProjection(const Coord &P, const Coord &A, const Coord &B) { return getLineProjection(P, Line(A, B - A)); }

Coord getSymmetry(const Coord &P, const Coord &O) { return O + O - P; } //P关于O的对称点

Coord getSymmetry(const Coord &P, const Line &L) { return getSymmetry(P, getLineProjection(P, L)); } //P关于L的对称点

Coord getLineIntersection(const Line &L1, const Line &L2) { return L1.point(Cross(L2.v, L1.p - L2.p) / Cross(L1.v, L2.v)); } //直线交点,须确保两直线相交

Coord getLineIntersection(const Coord &A1, const Coord &A2, const Coord &B1, const Coord &B2) { return getLineIntersection(Line(A1, A2 - A1), Line(B1, B2 - B1)); }

bool SegmentProperIntersection(const Coord &A1, const Coord &A2, const Coord &B1, const Coord &B2) //线段相交判定，交点不在一条线段的端点
{
	lf C1 = Cross(A2 - A1, B1 - A1), C2 = Cross(A2 - A1, B2 - A1),
	   C3 = Cross(B2 - B1, A1 - B1), C4 = Cross(B2 - B1, A2 - B1);
	return sgn(C1) * sgn(C2) < 0 && sgn(C3) * sgn(C4) < 0;
}

bool onSegment(const Coord &P, const Coord &A1, const Coord &A2) { return sgn(Dot(A1 - P, A2 - P)) < 0 && !sgn(Cross(A1 - P, A2 - P)); } //判断点是否在线段上，不包含端点

lf PolygonArea(const vector<Coord> &p) //计算多边形的有向面积，凸多边形即为面积
{
	lf s = 0;
	for (int i = 2; i < p.size(); ++i)
		s += Area2(p[0], p[i - 1], p[i]);
	return s / 2;
}

int inPolygon(const Coord &p, const vector<Coord> &poly) //点在多边形内的判定，转角法，正值为内部，0为外部，-1在边界上
{
	int ans = 0;
	for (int i = 0, k, d1, d2, n = poly.size(); i != n; ++i)
	{
		if (onSegment(p, poly[i], poly[(i + 1) % n]))
			return -1; //在边界上
		k = sgn(Cross(poly[(i + 1) % n] - poly[i], p - poly[i]));
		d1 = sgn(poly[i].Y - p.Y);
		d2 = sgn(poly[(i + 1) % n].Y - p.Y);
		if (k > 0 && d1 <= 0 && d2 > 0)
			++ans;
		if (k < 0 && d2 <= 0 && d1 > 0)
			--ans;
	}
	return ans;
}

vector<Coord> ConvexHull(vector<Coord> p, int collineation = 1) //获得凸包，不希望凸包的边上有输入点第二个参数传0
{
	vector<Coord> ans;
	sort(p.begin(), p.end(), cmpCoord); //先比横坐标再比纵坐标
	for (int i = 0; i < p.size(); ++i)  //求出下凸包
	{
		while (ans.size() > 1 && sgn(Area2(ans[ans.size() - 2], ans[ans.size() - 1], p[i])) < collineation)
			ans.pop_back();
		ans.push_back(p[i]);
	}
	for (int i = p.size() - 2, k = ans.size(); i >= 0; --i) //求出上凸包
	{
		while (ans.size() > k && -sgn(Area2(ans[ans.size() - 1], ans[ans.size() - 2], p[i])) < collineation)
			ans.pop_back();
		ans.push_back(p[i]);
	}
	if (p.size() > 1)
		ans.pop_back();
	return ans;
}

vector<Coord> cutPolygon(const vector<Coord> &poly, const Coord &A, const Coord &B) //用有向直线A->B切割多边形poly， 返回“左侧”。 如果退化，可能会返回一个单点或者线段，复杂度O(n^2)
{
	vector<Coord> newpoly;
	for (int i = 0, n = poly.size(); i != n; ++i)
	{
		Coord C = poly[i], D = poly[(i + 1) % n];
		if (sgn(Cross(B - A, C - A)) >= 0)
			newpoly.push_back(C);
		if (!sgn(Cross(B - A, C - D)))
		{
			Coord ip = getLineIntersection(Line(A, B - A), Line(C, D - C));
			if (onSegment(ip, C, D))
				newpoly.push_back(ip);
		}
	}
	return newpoly;
}

vector<Coord> getHalfPlaneIntersection(vector<Line> L) //半平面交
{
	sort(L.begin(), L.end(), cmpLine);  //按极角排序
	vector<Coord> p(L.size(), Coord()); //p[i]为q[i]和q[i+1]的交点
	int first = 0, last = 0;			//双端队列的第一个元素和最后一个元素
	vector<Line> q(L.size(), Line());   //双端队列
	q[0] = L[0];						//队列初始化为只有一个半平面L[0]

	for (int i = 0, n = L.size(); i != n; ++i)
	{
		while (first < last && !onLeft(p[last - 1], L[i]))
			--last;
		while (first < last && !onLeft(p[first], L[i]))
			++first;

		q[++last] = L[i];
		if (!sgn(Cross(q[last].v, q[last - 1].v)))
		{
			--last;
			if (onLeft(L[i].p, q[last]))
				q[last] = L[i];
		}
		if (first < last)
			p[last - 1] = getLineIntersection(q[last - 1], q[last]);
	}
	while (first < last && !onLeft(p[last - 1], q[first]))
		--last; //删除无用平面

	if (last - first <= 1)
		return vector<Coord>(); //空集
	p[last] = getLineIntersection(q[last], q[first]);
	return vector<Coord>(p.begin() + first, p.begin() + last + 1); //从deque复制到输出中
}

int getLineCircleIntersection(const Line &L, const Circle &C, vector<Coord> &sol)
{
	lf a = L.v.X,
	   b = L.p.X - C.c.X,
	   c = L.v.Y,
	   d = L.p.Y - C.c.Y,
	   e = a * a + c * c,
	   f = 2 * (a * b + c * d),
	   g = b * b + d * d - C.r * C.r,
	   delta = f * f - 4 * e * g;

	if (sgn(delta) < 0)
		return 0;
	if (!sgn(delta))
		return sol.push_back(L.point(-f / (2 * e))), 1;
	sol.push_back(L.point((-f - sqrt(delta)) / (2 * e)));
	sol.push_back(L.point((-f + sqrt(delta)) / (2 * e)));
	return 2;
}

int getCircleIntersection(const Circle &C1, const Circle &C2, vector<Coord> &sol)
{
	lf d = abs(C1.c - C2.c);

	if (!sgn(d))
		return sgn(C1.r - C2.r) ? 0 : -1; //重合返回-1

	if (sgn(C1.r + C2.r - d) < 0 || sgn(fabs(C1.r - C2.r) - d) > 0) //外离或内含
		return 0;

	lf a = arg(C2.c - C1.c),
	   da = acos((C1.r * C1.r + d * d - C2.r * C2.r) / (2 * C1.r * d));

	Coord p1 = C1.point(a - da), p2 = C1.point(a + da);

	sol.push_back(p1);

	if (p1 == p2)
		return 1; //相切
	return sol.push_back(p2), 2;
}

Line getTangent(const Coord &C, const Coord &P) { return Line(P, Normal(C - P)); } //圆心C，圆上一点P处切线

int getTangents(const Coord &p, const Circle &C, vector<Coord> &sol) //点到圆的切点，返回个数
{
	Coord u = p - C.c;
	lf d = abs(u);
	if (d < C.r)
		return 0;	  //点在圆内
	if (!sgn(d - C.r)) //点在圆上
		return sol.push_back(p), 1;
	lf base = arg(u), ang = acos(C.r / d);
	sol.push_back(C.point(base + ang));
	sol.push_back(C.point(base - ang));
	return 2;
}

int getTangents(Circle A, Circle &B, vector<Coord> &a, vector<Coord> &b) //公共切线的切点
{
	int cnt = 0;

	if (A.r < B.r)
		swap(A, B), swap(a, b); //有时需标记交换

	lf d = abs(A.c - B.c),
	   rdiff = A.r - B.r,
	   rsum = A.r + B.r;

	if (sgn(d - rdiff) < 0)
		return 0; //内含

	lf base = arg(B.c - A.c);

	if (!sgn(d) && !sgn(rdiff))
		return -1; //重合，无穷多条切线

	if (!sgn(d - rdiff)) //内切，外公切线
	{
		a.push_back(A.point(base));
		b.push_back(B.point(base));
		return 1;
	}

	//有外公切线的情形
	lf ang = acos(rdiff / d);
	a.push_back(A.point(base + ang));
	b.push_back(B.point(base + ang));
	a.push_back(A.point(base - ang));
	b.push_back(B.point(base - ang));
	cnt += 2;

	if (!sgn(d - rsum))
	{
		a.push_back(A.point(base));
		b.push_back(B.point(base + M_PI));
		++cnt;
	}
	else if (sgn(d - rsum) > 0)
	{
		lf ang_in = acos(rsum / d);
		a.push_back(A.point(base + ang_in));
		b.push_back(B.point(base + ang_in + M_PI));
		a.push_back(A.point(base - ang_in));
		b.push_back(B.point(base - ang_in + M_PI));
		cnt += 2;
	}
	return cnt;
}

lf AreaCircleWithTriangle(const Circle &C, Coord A, Coord B) //C和三角形OAB的相交面积，如果三角形顶点不在O上则把圆和三角形同时平移，直到有一个顶点在O上
{
	int sg = sgn(Cross(A, B));
	if (!sg || A == C.c || B == C.c)
		return 0;

	lf OA = abs(A - C.c), OB = abs(B - C.c), angle = Angle(A, B),
	   d = DistanceToLine(Coord(), A, B);

	if (sgn(OA - C.r) <= 0 && sgn(OB - C.r) <= 0)
		return Cross(A, B) / 2;

	if (sgn(OA - C.r) >= 0 && sgn(OB - C.r) >= 0 && sgn(d - C.r) >= 0)
		return sg * C.r * C.r * angle / 2;
	if (sgn(OA - C.r) >= 0 && sgn(OB - C.r) >= 0 && sgn(d - C.r) < 0)
	{
		Coord prj = getLineProjection(Coord(), A, B);
		if (!onSegment(prj, A, B))
			return sg * C.r * C.r * angle / 2;

		vector<Coord> p;
		Line L = Line(A, B - A);
		getLineCircleIntersection(L, C, p);

		lf s1 = C.r * C.r * angle / 2,
		   s2 = C.r * C.r * Angle(p[0], p[1]) / 2;
		s2 -= fabs(Cross(p[0], p[1]) / 2);
		s1 = s1 - s2;

		return sg * s1;
	}
	if (sgn(OB - C.r) < 0)
		swap(A, B);

	Line L = Line(A, B - A);
	vector<Coord> inter;
	getLineCircleIntersection(L, C, inter);
	Coord inter_point = inter[!onSegment(inter[0], A, B)];

	lf s = fabs(Cross(inter_point, A) / 2);
	s += C.r * C.r * Angle(inter_point, B) / 2;

	return s * sg;
}

lf AreaCircleWithPolygon(const Circle &C, const vector<Coord> &p)
{
	lf ans = 0;
	for (int i = 0; i < p.size(); ++i)
		ans += AreaCircleWithTriangle(C, p[i], p[(i + 1) % p.size()]);
	return fabs(ans);
}

Coord getGravityCenter(const vector<Coord> &p) //多边形重心
{
	Coord a(0, 0);
	lf am = 0, mj;
	for (int i = 0; i < p.size(); ++i)
	{
		mj = Cross(p[i], p[(i + 1) % p.size()]);
		a += mj * (p[i] + p[(i + 1) % p.size()]);
		am += mj;
	}
	return a / am / 3.0;
}

三维

typedef double lf;
const lf EPS = 1e-9, INF = 1e9;
int sgn(lf d) { return (d > EPS) - (d < -EPS); }
struct Coord3
{
	lf X, Y, Z;
	friend bool operator!=(const Coord3 &a, const Coord3 &b) { return sgn(a.X - b.X) || sgn(a.Y - b.Y) || sgn(a.Z - b.Z); }
	friend bool operator==(const Coord3 &a, const Coord3 &b) { return !(a != b); }
	Coord3 &operator+=(const Coord3 &b) { return X += b.X, Y += b.Y, Z += b.Z, *this; }
	friend Coord3 operator+(Coord3 a, const Coord3 &b) { return a += b; }
	Coord3 &operator-=(const Coord3 &b) { return X -= b.X, Y -= b.Y, Z -= b.Z, *this; }
	friend Coord3 operator-(Coord3 a, const Coord3 &b) { return a -= b; }
	Coord3 &operator*=(lf d) { return X *= d, Y *= d, Z *= d, *this; }
	friend Coord3 operator*(Coord3 a, lf d) { return a *= d; }
	friend Coord3 operator*(lf d, Coord3 a) { return a *= d; }
	Coord3 &operator/=(lf d) { return X /= d, Y /= d, Z /= d, *this; }
	friend Coord3 operator/(Coord3 a, lf d) { return a /= d; }
	friend lf Dot(const Coord3 &A, const Coord3 &B) { return A.X * B.X + A.Y * B.Y + A.Z * B.Z; }
	friend Coord3 Cross(const Coord3 &A, const Coord3 &B) { return {A.Y * B.Z - A.Z * B.Y, A.Z * B.X - A.X * B.Z, A.X * B.Y - A.Y * B.X}; }
	friend lf norm(const Coord3 &A) { return Dot(A, A); }
	friend lf abs(const Coord3 &A) { return sqrt(norm(A)); }
	friend lf Angle(const Coord3 &A, const Coord3 &B) { return acos(Dot(A, B) / abs(A) / abs(B)); }
	friend lf Area2(Coord3 A, Coord3 B, Coord3 C) { return abs(Cross(B - A, C - A)); }
	friend lf Volume6(Coord3 A, Coord3 B, Coord3 C, Coord3 D) { return Dot(D - A, Cross(B - A, C - A)); }	 //四面体体积
	friend Coord3 Centroid(Coord3 A, Coord3 B, Coord3 C, Coord3 D) { return (A + B + C + D) / 4.0; }		  //四面体的重心
	friend lf DistanceToPlane(Coord3 p, Coord3 p0, const Coord3 &n) { return Dot(p - p0, n) / abs(n); }		  //点p到平面p0-n的有向距离
	friend Coord3 getPlaneProjection(Coord3 p, Coord3 p0, const Coord3 &n) { return p - n * Dot(p - p0, n); } //点p在平面p0-n上的投影。n必须为单位向量
	friend Coord3 LinePlaneIntersection(Coord3 p1, Coord3 p2, Coord3 p0, Coord3 n)							  //直线p1-p2 与平面p0-n的交点，假设交点唯一存在
	{
		Coord3 v = p2 - p1;
		lf t = Dot(n, p0 - p1) / Dot(n, p2 - p1); //分母为0，直线与平面平行或在平面上
		return p1 + v * t;						  //如果是线段 判断t是否在0~1之间
	}
	friend lf DistanceToLine(Coord3 P, Coord3 A, Coord3 B) //点P到直线AB的距离
	{
		Coord3 v1 = B - A, v2 = P - A;
		return abs(Cross(v1, v2)) / abs(v1);
	}
	friend lf DistanceToSeg(Coord3 P, Coord3 A, Coord3 B) //点到线段的距离
	{
		if (A == B)
			return abs(P - A);
		Coord3 v1 = B - A, v2 = P - A, v3 = P - B;
		if (sgn(Dot(v1, v2)) < 0)
			return abs(v2);
		if (sgn(Dot(v1, v3)) > 0)
			return abs(v3);
		return fabs(DistanceToLine(P, A, B));
	}
	friend bool LineDistance3D(Coord3 p1, Coord3 u, Coord3 p2, Coord3 v, lf &s) //求异面直线 p1+s*u与p2+t*v的公垂线对应的s，如果平行|重合，返回0
	{
		lf b = Dot(u, u) * Dot(v, v) - Dot(u, v) * Dot(u, v);
		if (!sgn(b))
			return 0;
		lf a = Dot(u, v) * Dot(v, p1 - p2) - Dot(v, v) * Dot(u, p1 - p2);
		return s = a / b, 1;
	}
	friend bool SameSide(Coord3 p1, Coord3 p2, Coord3 a, Coord3 b) { return sgn(Dot(Cross(b - a, p1 - a), Cross(b - a, p2 - a))) >= 0; } //p1和p2是否在线段a-b的同侧
	friend bool PointInTri(Coord3 PP, Coord3 P[3])																						 //点P在三角形P0,P1,p中
	{
		return SameSide(PP, P[0], P[1], P[2]) &&
			   SameSide(PP, P[1], P[0], P[2]) &&
			   SameSide(PP, P[2], P[0], P[1]);
	}
	friend bool TriSegIntersection(Coord3 P[3], Coord3 A, Coord3 B, Coord3 &PP) //三角形P0P1p是否和线段AB相交，如有则为PP
	{
		Coord3 n = Cross(P[1] - P[0], P[2] - P[0]);
		if (sgn(Dot(n, B - A)) == 0)
			return false;						 //线段A-B和平面P0P1p平行或共面
		lf t = Dot(n, P[0] - A) / Dot(n, B - A); //平面A和直线P1-p有惟一交点
		if (sgn(t) < 0 || sgn(t - 1) > 0)
			return false; //不在线段AB上
		return PointInTri(PP = A + (B - A) * t, P);
	}
	friend bool TriTriIntersection(Coord3 T1[3], Coord3 T2[3]) //空间两三角形是否相交
	{
		Coord3 P;
		for (int i = 0; i < 3; ++i)
			if (TriSegIntersection(T1, T2[i], T2[(i + 1) % 3], P) ||
				TriSegIntersection(T2, T1[i], T1[(i + 1) % 3], P))
				return 1;
		return 0;
	}
	friend lf SegSegDistance(Coord3 a1, Coord3 b1, Coord3 a2, Coord3 b2, Coord3 &ans1, Coord3 &ans2) //空间两直线上最近点对 返回最近距离 点对保存在ans1 ans2中
	{
		Coord3 v1 = (a1 - b1), v2 = (a2 - b2);
		Coord3 N = Cross(v1, v2);
		Coord3 ab = (a1 - a2);
		lf ans = Dot(N, ab) / abs(N);
		Coord3 d1 = b1 - a1, d2 = b2 - a2, cd = Cross(d1, d2);
		lf nd = norm(cd), t1 = Dot(Cross(a2 - a1, d2), cd) / nd, t2 = Dot(Cross(a2 - a1, d1), cd) / nd;
		return ans1 = a1 + (b1 - a1) * t1, ans2 = a2 + (b2 - a2) * t2, fabs(ans);
	}
	friend bool InsideWithMinDistance(Coord3 PP, Coord3 *P, lf dist) //判断PP是否在三角形P中，并且到三条边的距离都至少为dist。保证P,A,B,C共面
	{
		return PointInTri(PP, P) &&
				   DistanceToLine(PP, P[0], P[1]) >= dist ||
			   DistanceToLine(PP, P[1], P[2]) >= dist ||
			   DistanceToLine(PP, P[2], P[0]) >= dist;
	}
	friend lf minBall(const vector<Coord3> &data, const lf eps = EPS * 1e-3) //模拟退火求最小球覆盖，EPS玄学调整；返回半径
	{
		lf step = 1, ans = INF;
		for (Coord3 z{0, 0, 0}; step > eps; step *= 0.99)
		{
			int s = 0;
			for (int i = 0; i < data.size(); ++i)
				if (abs(z - data[s]) < abs(z - data[i]))
					s = i;
			ans = min(ans, abs(z - data[s]));
			z -= (z - data[s]) * step;
		}
		return ans;
	}
};
struct Sphere
{
	Coord3 o;
	lf r;
	Coord3 point(lf lat, lf lng) const //经纬度确定球面上一点
	{
		lat *= M_PI / 180;
		lng *= M_PI / 180;
		return o + r * Coord3{cos(lat) * cos(lng), cos(lat) * sin(lng), sin(lat)};
	}
};
struct ConvexPolyhedron //空间多边形和凸包问题
{
	struct Face
	{
		int v[3];
		Face(int a, int b, int c) { v[0] = a, v[1] = b, v[2] = c; }
		Coord3 Normal(const vector<Coord3> &P) const { return Cross(P[v[1]] - P[v[0]], P[v[2]] - P[v[0]]); }
		bool CanSee(const vector<Coord3> &P, int i) const { return Dot(P[i] - P[v[0]], Normal(P)) > 0; } //f是否能看见P[i]
	};
	vector<Face> faces;
	vector<Coord3> p;
	ConvexPolyhedron(vector<Coord3> P) : p(P)
	{
		for (int i = 0; i < p.size(); ++i)
			P[i] += Coord3{randEPS(), randEPS(), randEPS()};
		vector<vector<int>> vis(P.size(), vector<int>(P.size()));
		faces.push_back(Face(0, 1, 2)); //由于已经进行扰动，前三个点不共线
		faces.push_back(Face(2, 1, 0));
		for (int i = 3; i < P.size(); ++i)
		{
			vector<Face> next;
			for (int j = 0; j < faces.size(); ++j) //计算每条边的“左面”的可见性
			{
				Face &f = faces[j];
				int res = f.CanSee(P, i);
				if (!res)
					next.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 < faces.size(); ++j)
				for (int k = 0; k < 3; ++k)
				{
					int a = faces[j].v[k], b = faces[j].v[(k + 1) % 3];
					if (vis[a][b] != vis[b][a] && vis[a][b]) //(a,b)是分界线，左边对P[i]可见
						next.push_back(Face(a, b, i));
				}
			swap(faces, next);
		}
	}
	lf randEPS() { return (rand() / lf(RAND_MAX) - 0.5) * EPS; }
	Coord3 centroid() //三维凸包重心
	{
		Coord3 C = p[0], tot{0, 0, 0};
		lf totv = 0;
		for (int i = 0; i < faces.size(); ++i)
		{
			Coord3 p1 = p[faces[i].v[0]], p2 = p[faces[i].v[1]], p3 = p[faces[i].v[2]];
			lf v = -Volume6(p1, p2, p3, C);
			totv += v;
			tot += v * Centroid(p1, p2, p3, C);
		}
		return tot / totv;
	}
	lf dist(Coord3 C) //凸包内一点到表面最近距离
	{
		lf ans = INF;
		for (int i = 0; i < faces.size(); ++i)
		{
			Coord3 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;
	}
};