• Hnoi2004 金属包裹


    传送门

    三维凸包模板题……只是听了听计算几何的课之后心血来潮想写的……

    我的做法很无脑是吧……暴力枚举三个点组成的三角形,然后枚举剩下的点,判断其余点是否都在这个三角形的同一侧,是的话则说明这个三角形是凸包的一个面。

    理论复杂度应该是$O(n^4)$,不过看上去跑得飞快?人帅自带小常数哈哈

    这个故事告诉我们:大力出奇迹……

     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 */
    View Code

    话说一开始忘了怎么求平面的法向量,手推了个公式然后发现搞出NaN了……因为我推的公式默认法向量在z维的长度不为0,但其实有很多平面的法向量是垂直于z轴的,然后就除零爆炸了……又脑补了很久,后来才想起来两个向量的叉积必定同时垂直于这两个向量……这人没救了

    还有一件事,这题卡精度,随机扰乱搞得太大会炸精度……计算几何毁我青春

  • 相关阅读:
    Codeforces 543E. Listening to Music
    UOJ #138. 【UER #3】开学前的涂鸦
    bzoj 3569: DZY Loves Chinese II
    bzoj 2428: [HAOI2006]均分数据
    bzoj 4589: Hard Nim
    UOJ #119. 【UR #8】决战圆锥曲线
    spoj5973
    codeforces555E
    poj1275
    bzoj4152
  • 原文地址:https://www.cnblogs.com/hzoier/p/6418036.html
Copyright © 2020-2023  润新知