• 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轴的,然后就除零爆炸了……又脑补了很久,后来才想起来两个向量的叉积必定同时垂直于这两个向量……这人没救了

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

  • 相关阅读:
    vim复制
    嵌入式Linux学习(二)
    (Java实现) 洛谷 P1042 乒乓球
    (Java实现) 洛谷 P1042 乒乓球
    (Java实现) 洛谷 P1071 潜伏者
    (Java实现) 洛谷 P1071 潜伏者
    (Java实现) 洛谷 P1025 数的划分
    (Java实现)洛谷 P1093 奖学金
    (Java实现)洛谷 P1093 奖学金
    Java实现 洛谷 P1064 金明的预算方案
  • 原文地址:https://www.cnblogs.com/hzoier/p/6418036.html
Copyright © 2020-2023  润新知