tetrahedron
题意：
求解一个四面体的内切球。
解法：
首先假设内切球球心为$(x0,x1,x2)$，可以用$r = \frac{3V}{S_1+S_2+S_3+S_4}$得出半径，
这样对于四个平面列出三个方程，解得
$x_n = \sum_{i=0}^3{Ai_{x_n} \cdot S_i } / (S_1 + S_2 + S_3 + S_4)$
这样，即可得出内切球。
时间复杂度$O(1)$。
1 #include <iostream> 2 #include <cstdio> 3 #include <cstring> 4 #include <cmath> 5 6 #define LD double 7 #define sqr(x) ((x)*(x)) 8 #define eps 1e-13 9 10 using namespace std; 11 12 struct node 13 { 14 LD x,y,z; 15 void scan() 16 { 17 scanf("%lf%lf%lf",&x,&y,&z); 18 } 19 void print() 20 { 21 printf("%.4lf %.4lf %.4lf\n",x,y,z); 22 } 23 LD length() 24 { 25 return sqrt(sqr(x)+sqr(y)+sqr(z)); 26 } 27 node operator+(const node &tmp) 28 { 29 return (node){x+tmp.x,y+tmp.y,z+tmp.z}; 30 } 31 node operator-(const node &tmp) 32 { 33 return (node){x-tmp.x,y-tmp.y,z-tmp.z}; 34 } 35 node operator/(LD tmp) 36 { 37 return (node){x/tmp,y/tmp,z/tmp}; 38 } 39 node operator*(LD tmp) 40 { 41 return (node){x*tmp,y*tmp,z*tmp}; 42 } 43 }; 44 45 node cross(node a,node b) 46 { 47 node ans; 48 ans.x = a.y*b.z - b.y*a.z; 49 ans.y = b.x*a.z - a.x*b.z; 50 ans.z = a.x*b.y - b.x*a.y; 51 return ans; 52 } 53 54 LD dist(node a,node b) 55 { 56 return (b-a).length(); 57 } 58 59 LD dot(node a,node b) 60 { 61 return a.x*b.x + a.y*b.y + a.z*b.z; 62 } 63 64 LD get_angle(node a,node b) 65 { 66 LD tmp = dot(a,b)/a.length()/b.length(); 67 return acos(tmp); 68 } 69 70 node get_node(node A,node B,node C) 71 { 72 LD Lth = (B-A).length() + (C-A).length() + (C-B).length(); 73 cout << sqr(Lth-4) << endl; 74 LD r = fabs(cross(B-A,C-A).length()) / Lth; 75 cout << r*r << endl; 76 node v1 = C-A; 77 node v2 = B-A; 78 node v = (v1+v2)/(v1+v2).length(); 79 LD d = (C-A).length()/2; 80 LD L = sqrt(sqr(d)+sqr(r)); 81 v = v*L; 82 return A+v; 83 } 84 85 int main() 86 { 87 node A,B,C,D; 88 while(~scanf("%lf%lf%lf",&A.x,&A.y,&A.z)) 89 { 90 B.scan(); 91 C.scan(); 92 D.scan(); 93 if(fabs(dot(cross(B-A,C-A),D-A)) < eps) 94 { 95 puts("O O O O"); 96 continue; 97 } 98 LD S1 = fabs(cross(B-D,C-D).length())/2; 99 LD S2 = fabs(cross(D-A,C-A).length())/2; 100 LD S3 = fabs(cross(B-A,D-A).length())/2; 101 LD S4 = fabs(cross(B-A,C-A).length())/2; 102 LD Ve = fabs(dot(cross(B-A,C-A),D-A))/6; 103 LD R = 3*Ve / ((S1+S2+S3+S4)); 104 node ans; 105 ans.x = (S1*A.x + S2*B.x + S3*C.x + S4*D.x)/(S1+S2+S3+S4); 106 ans.y = (S1*A.y + S2*B.y + S3*C.y + S4*D.y)/(S1+S2+S3+S4); 107 ans.z = (S1*A.z + S2*B.z + S3*C.z + S4*D.z)/(S1+S2+S3+S4); 108 printf("%.4lf %.4lf %.4lf %.4lf\n",ans.x,ans.y,ans.z,R); 109 } 110 return 0; 111 } 112 /* 113 0 0 0 2 0 0 0 0 2 0 2 0 114 0 0 0 2 0 0 3 0 0 4 0 0 115 */
