Hnoi2004 金属包裹

 1 #include<cstdio>
 2 #include<cstring>
 3 #include<cmath>
 4 #include<algorithm>
 5 #include<iostream>
 6 #include<fstream>
 7 #include<iomanip>
 8 using namespace std;
 9 const int maxn=110;
10 const long double eps=1e-15;
11 struct Point{
12     long double x,y,z;
13     Point(long double x=0.0,long double y=0.0,long double z=0.0):x(x),y(y),z(z){}
14     Point operator-(const Point &a)const{return Point(a.x-x,a.y-y,a.z-z);}//A - B = B到A的位移
15     Point operator/(const long double &a)const{return Point(x/a,y/a,z/a);}
16 }a[maxn];
17 typedef Point Vector;
18 long double noise();
19 long double Dot(const Vector&,const Vector&);
20 Vector Cross(const Vector&,const Vector&);
21 long double Length(const Vector&);
22 long double Area2(const Point&,const Point&,const Point&);
23 Vector LawVector(const Vector&,const Vector&);
24 long double Distance(const Point&,const Point&,const Vector&);
25 ifstream fin("enwrap.in");
26 ofstream fout("enwrap.out");
27 long double ans=0.0;
28 int n,t;
29 int main(){
30     fin>>n;
31     for(int i=1;i<=n;i++){
32         fin>>a[i].x>>a[i].y>>a[i].z;
33         a[i].x+=noise();
34         a[i].y+=noise();
35         a[i].z+=noise();
36     }
37     Vector A;
38     long double d;
39     bool bad;
40     for(int i=1;i<=n;i++)for(int j=1;j<i;j++)for(int k=1;k<j;k++){
41         bad=false;
42         t=0;
43         A=LawVector(a[j]-a[i],a[k]-a[i]);
44         for(int l=1;l<=n;l++)if(i!=l&&j!=l&&k!=l){
45             d=Distance(a[l],a[i],A);
46             if(!t)t=(d<0?-1:1);
47             else if(t*d<0){
48                 bad=true;
49                 break;
50             }
51         }
52         if(!bad)ans+=Area2(a[i],a[j],a[k]);
53     }
54     ans/=2.0;
55     fout<<fixed<<ans;
56     return 0;
57 }
58 long double noise(){
59     static int a=1213,b=678424137,p=998244353,x=753497501;
60     x=a*x+b;x%=p;
61     if(x<0)x+=p;
62     return (long double)(x-p)/p*5e-10;
63 }
64 inline long double Dot(const Vector &A,const Vector &B){return A.x*B.x+A.y*B.y+A.z*B.z;}
65 inline Vector Cross(const Vector &A,const Vector &B){return Vector(A.y*B.z-B.y*A.z,A.z*B.x-B.z*A.x,A.x*B.y-B.x*A.y);}
66 inline long double Length(const Vector &A){return sqrt(Dot(A,A));}
67 inline long double Area2(const Point &A,const Point &B,const Point &C){return Length(Cross(B-A,C-A));}
68 inline Vector LawVector(const Vector &A,const Vector &B){
69     Vector n=Cross(A,B);
70     return n/Length(n);
71 }//两个向量的叉积一定同时垂直于这两个向量
72 inline long double Distance(const Point &A,const Point &P,const Vector &n){return Dot(A-P,n);}//(P,n)是平面的点法式,n是单位向量
73 /*
74 三维凸包——暴力法
75 暴力枚举三个点组成的三角形，
76 判断其他点是否都在这个三角形的同侧，
77 是则说明这个三角形一定是凸包的一个面。
78 */