uva 11836(两个四面体最短距离）

题意：给出两个四面体的四个顶点坐标(保证四两个面体不相交，且体积大于0)，输出两个四面体的最近距离。
题解：要求最近距离分3种情况：
1.顶点到另一个四面体的某面距离
2.顶点到另一个四面体的某条边距离
3.两个四面体的边与边的最短距离(异面线段最短距离)
异面线段最短距离又分两种情况，两条线段所在直线的公垂线的最小距离点对是否都在线段上，如果在那么公垂线的距离就是异面线段最短距离，否则就是两对端点计算四个点距的最短距离了。

#include <cstdio>#include <cmath>#include <algorithm>using namespace std;const double eps = 1e-9;const double PI = acos(-1);struct Point3 {    double x, y, z;    Point3(double x = 0, double y = 0, double z = 0):x(x), y(y), z(z) {}};typedef Point3 Vector3;int dcmp(double x) {    if (fabs(x) < eps)        return 0;    return x < 0 ? -1 : 1;}Vector3 operator + (Vector3 A, Vector3 B) { return Vector3(A.x + B.x, A.y + B.y, A.z + B.z); }Vector3 operator - (Vector3 A, Vector3 B) { return Vector3(A.x - B.x, A.y - B.y, A.z - B.z); }Vector3 operator * (Vector3 A, double b) { return Vector3(A.x * b, A.y * b, A.z * b); }Vector3 operator / (Vector3 A, double b) { return Vector3(A.x / b, A.y / b, A.z / b); }bool operator == (Vector3 A, Vector3 B) { return !dcmp(A.x - B.x) && !dcmp(A.y - B.y); }double Dot (Vector3 A, Vector3 B) { return A.x * B.x + A.y * B.y + A.z * B.z; }double Length(Vector3 A) { return sqrt(Dot(A, A)); }double Angle(Vector3 A, Vector3 B) { return acos(Dot(A, B) / Length(A) / Length(B)); }double torad(double deg) { return deg / 180.0 * PI; }Vector3 Cross(Vector3 A, Vector3 B) { return Vector3(A.y * B.z - A.z * B.y, A.z * B.x - A.x * B.z, A.x * B.y - A.y * B.x); }double Area2(Point3 A, Point3 B, Point3 C) { return Length(Cross(B - A, C - A)); }double DistanceToPlane(const Point3& p, const Point3& p0, const Vector3& n) {    return fabs(Dot(p - p0, n));}Point3 GetPlaneProjection(const Point3& p, const Point3& p0, const Vector3& n) {    return p - n * Dot(p - p0, n);}double DistanceToSegment(Point3 P, Point3 A, Point3 B) {    if (A == B)        return Length(P - A);    Vector3 v1 = B - A, v2 = P - A, v3 = P - B;    if (dcmp(Dot(v1, v2)) < 0)        return Length(v2);    if (dcmp(Dot(v1, v3)) > 0)        return Length(v3);    return Length(Cross(v1, v2)) / Length(v1);}bool PointInTri(Point3 P, Point3 P0, Point3 P1, Point3 P2) {    double area1 = Area2(P, P0, P1);    double area2 = Area2(P, P1, P2);    double area3 = Area2(P, P2, P0);    return dcmp(area1 + area2 + area3 - Area2(P0, P1, P2)) == 0;}double MinSegmentToSegmentDis(Point3 p1, Point3 p2, Point3 p3, Point3 p4) {    Vector3 v = Cross(p2 - p1, p4 - p3); //公垂直线    double A = Dot(Cross(p2 - p1, v), p3 - p1);    double B = Dot(Cross(p2 - p1, v), p4 - p1);    double C = Dot(Cross(p3 - p4, v), p1 - p4);    double D = Dot(Cross(p3 - p4, v), p2 - p4);    if (dcmp(A * B) <= 0 && dcmp(C * D) <= 0)        return fabs(Dot(p1 - p3, v)) / Length(v);    return min(min(DistanceToSegment(p1, p3, p4), DistanceToSegment(p2, p3, p4)), min(DistanceToSegment(p3, p1, p2), DistanceToSegment(p4, p1, p2)));   }Point3 P1[4], P2[4];int main() {    int t;    scanf("%d", &t);    while (t--) {        for (int i = 0; i < 4; i++)            scanf("%lf%lf%lf", &P1[i].x, &P1[i].y, &P1[i].z);        for (int i = 0; i < 4; i++)            scanf("%lf%lf%lf", &P2[i].x, &P2[i].y, &P2[i].z);        double res = 1e18;        for (int i = 0; i < 4; i++) {            for (int j = 0; j < 4; j++) {                Vector3 v1 = P2[(j + 1) % 4] - P2[j], v2 = P2[(j + 2) % 4] - P2[j];                double z = (v2.y * v1.x - v1.y * v2.x) / (v1.y * v2.z - v2.y * v1.z);                double y;                if (dcmp(v1.y))                    y = (-z * v1.z - v1.x) / v1.y;                else                    y = (-z * v2.z - v2.x) / v2.y;                Vector3 n = Vector3(1.0, y, z);                n = n / Length(n);                Point3 p = GetPlaneProjection(P1[i], P2[j], n);                if (PointInTri(p, P2[j], P2[(j + 1) % 4], P2[(j + 2) % 4]))                    res = min(res, DistanceToPlane(P1[i], P2[j], n));            }            for (int j = i + 1; j < 4; j++)                res = min(res, DistanceToSegment(P1[i], P2[i], P2[j]));        }        for (int i = 0; i < 4; i++) {            for (int j = 0; j < 4; j++) {                Vector3 v1 = P1[(j + 1) % 4] - P1[j], v2 = P1[(j + 2) % 4] - P1[j];                double z = (v2.y * v1.x - v1.y * v2.x) / (v1.y * v2.z - v2.y * v1.z);                double y;                if (dcmp(v1.y))                    y = (-z * v1.z - v1.x) / v1.y;                else                    y = (-z * v2.z - v2.x) / v2.y;                Vector3 n = Vector3(1.0, y, z);                n = n / Length(n);                Point3 p = GetPlaneProjection(P2[i], P1[j], n);                if (PointInTri(p, P1[j], P1[(j + 1) % 4], P1[(j + 2) % 4]))                    res = min(res, DistanceToPlane(P2[i], P1[j], n));            }            for (int j = i + 1; j < 4; j++)                res = min(res, DistanceToSegment(P2[i], P1[i], P1[j]));        }        for (int i = 0; i < 4; i++)            for (int j = i + 1; j < 4; j++)                for (int k = 0; k < 4; k++)                    for (int l = k + 1; l < 4; l++)                        res = min(res, MinSegmentToSegmentDis(P1[i], P1[j], P2[k], P2[l]));        printf("%.2lf\n", res);    }    return 0;}