UVALive4973 CERC2010A Ardenia
分类讨论的情况不难想，难点在于判断各种垂线垂足是否在线段上。设bl1、bl2为两个线段上公垂线垂足位置的比例值，x为p0的公垂线垂足X坐标，则：
x = (p1.x - p0.x) * bl1 + p0.x
同理可得其他坐标。公垂线向量与两线段向量点积为0可得两个方程，求得bl1和bl2皆在0~1范围内则公垂线垂足都在线段上。
1 #include<stdio.h> 2 #include<stdlib.h> 3 #include<string.h> 4 typedef long long LL; 5 const double eps = 1e-8; 6 int dcmp(double d) 7 { 8 if(d < -eps) return -1; 9 return d > eps ? 1 : 0; 10 } 11 struct Point3 12 { 13 Point3() 14 { 15 x = y = z = 0; 16 } 17 Point3(LL a, LL b, LL c) 18 { 19 x = a, y = b, z = c; 20 } 21 Point3 operator*(const Point3 &p)const 22 { 23 return Point3(y * p.z - z * p.y, 24 z * p.x - x * p.z, 25 x * p.y - y * p.x); 26 } 27 Point3 operator-(const Point3 &p)const 28 { 29 return Point3(x - p.x, y - p.y, z - p.z); 30 } 31 Point3 operator+(const Point3 &p)const 32 { 33 return Point3(x + p.x, y + p.y, z + p.z); 34 } 35 Point3 operator-()const 36 { 37 return Point3(-x, -y, -z); 38 } 39 LL dot(const Point3 &p)const 40 { 41 return x * p.x + y * p.y + z * p.z; 42 } 43 LL x, y, z; 44 } p[4]; 45 LL gcd(LL a, LL b) 46 { 47 return b ? gcd(b, a % b) : a; 48 } 49 struct Dis 50 { 51 LL fz, fm; 52 void yf() 53 { 54 if(fz == 0) 55 { 56 fm = 1; 57 return; 58 } 59 LL t = gcd(fz, fm); 60 fz /= t; 61 fm /= t; 62 } 63 bool operator<(const Dis &p)const 64 { 65 return fz * p.fm < p.fz * fm; 66 } 67 }; 68 inline LL Sqr(LL a) 69 { 70 return a * a; 71 } 72 bool Paral(Point3 a, Point3 b, Point3 c, Point3 d) 73 { 74 Point3 tmp = (b - a) * (d - c); 75 return !tmp.dot(tmp); 76 } 77 bool JudgeCZ(Point3 a, Point3 b, Point3 c, Point3 d) 78 { 79 LL A1, B1, C1, A2, B2, C2; 80 A1 = (b - a).dot(b - a); 81 B1 = -(d - c).dot(b - a); 82 C1 = -(a - c).dot(b - a); 83 A2 = (b - a).dot(d - c); 84 B2 = -(d - c).dot(d - c); 85 C2 = -(a - c).dot(d - c); 86 double bl1 = dcmp(A2 * B1 - A1 * B2) ? ((double)C2 * B1 - C1 * B2) / (A2 * B1 - A1 * B2) : (A1 ? (double)C1 / A1 : (A2 ? (double)C2 / A2 : 0)); 87 double bl2 = dcmp(B2 * A1 - B1 * A2) ? ((double)C2 * A1 - C1 * A2) / (B2 * A1 - B1 * A2) : (B1 ? (double)C1 / B1 : (B2 ? (double)C2 / B2 : 0)); 88 return bl1 > -eps && bl1 < 1 + eps && bl2 > -eps && bl2 < 1 + eps; 89 } 90 Dis CalPtoL(Point3 p, Point3 a, Point3 b) 91 { 92 Point3 t = (a - p) * (b - a); 93 Dis tmp; 94 tmp.fz = t.dot(t); 95 tmp.fm = (b - a).dot(b - a); 96 tmp.yf(); 97 return tmp; 98 } 99 Dis CalPtoP(Point3 a, Point3 b) 100 { 101 Dis tmp; 102 tmp.fz = (b - a).dot(b - a); 103 tmp.fm = 1; 104 tmp.yf(); 105 return tmp; 106 } 107 Dis min(Dis a, Dis b) 108 { 109 return a < b ? a : b; 110 } 111 int main() 112 { 113 int t, i; 114 Dis ans; 115 for(scanf("%d", &t); t -- ;) 116 { 117 ans.fz = 1, ans.fm = 0; 118 for(i = 0; i < 4; ++ i) 119 scanf("%lld%lld%lld", &p[i].x, &p[i].y, &p[i].z); 120 if(!Paral(p[0], p[1], p[2], p[3]) && JudgeCZ(p[0], p[1], p[2], p[3])) 121 { 122 Point3 t = (p[1] - p[0]) * (p[3] - p[2]); 123 ans.fz = Sqr((p[2] - p[0]).dot(t)); 124 ans.fm = t.dot(t); 125 ans.yf(); 126 } 127 else 128 { 129 if(JudgeCZ(p[0], p[0], p[2], p[3])) 130 ans = min(ans, CalPtoL(p[0], p[2], p[3])); 131 if(JudgeCZ(p[1], p[1], p[2], p[3])) 132 ans = min(ans, CalPtoL(p[1], p[2], p[3])); 133 if(JudgeCZ(p[0], p[1], p[2], p[2])) 134 ans = min(ans, CalPtoL(p[2], p[0], p[1])); 135 if(JudgeCZ(p[0], p[1], p[3], p[3])) 136 ans = min(ans, CalPtoL(p[3], p[0], p[1])); 137 ans = min(ans, CalPtoP(p[0], p[2])); 138 ans = min(ans, CalPtoP(p[0], p[3])); 139 ans = min(ans, CalPtoP(p[1], p[2])); 140 ans = min(ans, CalPtoP(p[1], p[3])); 141 } 142 printf("%lld %lld\n", ans.fz, ans.fm); 143 } 144 return 0; 145 }
