[HNOI2004]最佳包裹

C++14(g++5.4) 解法, 执行用时: 3ms, 内存消耗: 552K, 提交时间: 2020-09-12 14:19:52

#include<bits/stdc++.h>
#define eps 1e-11
using namespace std;

struct Point3
{
    double x,y,z;
    Point3 (double x=0,double y=0,double z=0):x(x),y(y),z(z){}
};

typedef Point3 Vector3;
const int N=510;

Vector3 operator + (Vector3 A, Vector3 B){return Vector3{A.x+B.x,A.y+B.y,A.z+B.z};}
Vector3 operator - (Vector3 A, Vector3 B){return Vector3{A.x-B.x,A.y-B.y,A.z-B.z};}
Vector3 operator * (Vector3 A, double B){return Vector3{A.x*B,A.y*B,A.z*B};}
Vector3 operator / (Vector3 A, double B){return Vector3{A.x/B,A.y/B,A.z/B};}

double Dot(Vector3 A, Vector3 B){return A.x*B.x+A.y*B.y+A.z*B.z;}
Vector3 Cross(Vector3 A,Vector3 B) {return Vector3(A.y*B.z-A.z*B.y,A.z*B.x-A.x*B.z,A.x*B.y-A.y*B.x);}
double len(Vector3 A){return sqrt(Dot(A,A));}
double angle(Vector3 A, Vector3 B){return acos(Dot(A,B)/len(A)/len(B));}

double rand01(){return rand()/(double)RAND_MAX;}
double randeps(){return (rand01()-0.5)*eps;}
Point3 add_noise(Point3 p){return Point3{p.x+randeps(),p.y+randeps(),p.z+randeps()};}

struct Face{
    int v[3];
    Vector3 normal(Point3 *P)const{
        return Cross(P[v[1]]-P[v[0]],P[v[2]]-P[v[0]]);
    }
    int cansee(Point3 *P,int i)const{
        return Dot(P[i]-P[v[0]],normal(P))>0?1:0;
    }
};

vector<Face>F;
Point3 P[N],Q[N];

vector<Face> CH3D(Point3 *P,int n){
    vector<Face> cur;
    int vis[N][N];
    cur.push_back((Face){0,1,2});
    cur.push_back((Face){2,1,0});
    for(int i=3;i<n;i++){
        vector<Face> next;
        for(int j=0;j<cur.size();j++){
            Face &f=cur[j];
            int res=f.cansee(P,i);
            if(!res) next.push_back(f);
            for(int k=0;k<3;k++) vis[f.v[k]][f.v[(k+1)%3]]=res;
        }
        for(int j=0;j<cur.size();j++){
            for(int k=0;k<3;k++){
                int a=cur[j].v[k],b=cur[j].v[(k+1)%3];
                if(vis[a][b]!=vis[b][a]&&vis[a][b]) next.push_back((Face){a,b,i});
            }
        }
        cur=next;
    }
    return cur;
}

//三角形面积*2
double S(Point3 a,Point3 b,Point3 c) {
    return len(Cross(b-a,c-a));    //注意要取lenth
}
//四面体有向体积*6  混合积
double V(Point3 a,Point3 b,Point3 c,Point3 d) {
    return Dot(d-a,Cross(b-a,c-a));
}
//表面积
double Ssum(Point3 *P,int n) {
    double ans=0;
    if (n==3) {
        ans=S(P[0],P[1],P[2]);
        return ans/2.0;
    }
    for (int i=0;i<F.size();i++)
        ans=ans+S(P[F[i].v[0]],P[F[i].v[1]],P[F[i].v[2]])/2.0;
    return ans;
}
//体积
double Vsum(Point3 *P) {
    double ans=0;
    Point3 tmp(0,0,0);
    for (int i=0;i<F.size();i++)
        ans+=V(tmp,P[F[i].v[0]],P[F[i].v[1]],P[F[i].v[2]]);
    return fabs(ans/6.0);
}
bool same(int s,int t,Point3 *P) {
    Point3 &a=P[F[s].v[0]];
    Point3 &b=P[F[s].v[1]];
    Point3 &c=P[F[s].v[2]];
    return fabs(V(a,b,c,P[F[t].v[0]]))<eps &&
            fabs(V(a,b,c,P[F[t].v[1]]))<eps &&
            fabs(V(a,b,c,P[F[t].v[2]]))<eps;
}
//表面多边形个数
int facecnt(Point3 *P) {
    int flag,ans=0;
    for (int i=0;i<F.size();i++) {
        flag=1;
        for (int j=0;j<i;j++)
            if (same(i,j,P)) {
                flag=0;
                break;
            }
        ans+=flag;
    }
    return ans;
}
//三维凸包重心
Point3 centre(Point3 *P) {
    Point3 ans(0,0,0),o(0,0,0);
    double all=0;
    for (int i=0;i<F.size();i++) {
        double vol=V(o,P[F[i].v[0]],P[F[i].v[1]],P[F[i].v[2]]);
        ans=ans+(o+P[F[i].v[0]]+P[F[i].v[1]]+P[F[i].v[2]])/4.0*vol;
        all+=vol;
    }
    ans=ans/all;
    return ans;
}
//点到面的距离
double dis_face(Point3 p,Point3 *P,int i) {
    return fabs(V(P[F[i].v[0]],P[F[i].v[1]],P[F[i].v[2]],p)/len(Cross(P[F[i].v[1]]-P[F[i].v[0]],P[F[i].v[2]]-P[F[i].v[0]])));
}

int main()
{
    int n;
    scanf("%d",&n);
    for (int i=0;i<n;i++)
        scanf("%lf%lf%lf",&P[i].x,&P[i].y,&P[i].z);
    for(int i=0;i<n;i++) Q[i]=add_noise(P[i]);
    F=CH3D(Q,n);
    printf("%.6f\n",fabs(Ssum(P,n)));
    return 0;
}

C++ 解法, 执行用时: 3ms, 内存消耗: 436K, 提交时间: 2021-09-09 21:37:00

#include<bits/stdc++.h>

using namespace std;
typedef pair<double,double>pdd;
#define x first
#define y second
const int N = 250;
const double eps = 1e-12;
const double pi = acos(-1);

bool g[N][N];
int n,m;

double rand_eps(){
    return ((double)rand() / RAND_MAX - 0.5) * eps;
}

struct Point{
    double x,y,z;
    void skank(){
        x+=rand_eps();y+=rand_eps();z+=rand_eps();
    }
    Point operator + (Point b){
        return {x+b.x,y+b.y,z+b.z};
    }
    Point operator - (Point b){
        return {x-b.x,y-b.y,z-b.z};
    }
    Point operator * (Point b){
        return {y*b.z - b.y*z,z*b.x - b.z*x, x*b.y - y*b.x}; 
    }
    double operator & (Point t){
        return x*t.x+y*t.y+z*t.z;
    }
    double len(){
        return sqrt(x*x+y*y+z*z);
    }
}q[N];

struct Plane{
    int v[3];
    Point norm(){
        return (q[v[1]] - q[v[0]]) * (q[v[2]] - q[v[0]]);
    }
    double area(){
        return norm().len() / 2;
    }
    bool above(Point a){
        return ((a- q[v[0]]) & norm()) >= 0;
    }
}plane[N], np[N]; 

void get_conx_3d(){
    plane[m++] = {0,1,2};
    plane[m++] = {2,1,0};
    for(int i=3;i<n;i++){
        int cnt=0;
        for(int j=0;j<m;j++){
            bool t = plane[j].above(q[i]);
            if(!t)np[cnt++] = plane[j];
            for(int k=0;k<3;k++){
                g[plane[j].v[k]][plane[j].v[(k+1)%3]] = t;
            }
        }
        for(int j=0;j<m;j++){
            for(int k=0;k<3;k++){
                int a = plane[j].v[k], b = plane[j].v[(k+1)%3];
                if(g[a][b] && !g[b][a]){
                    np[cnt++] = {a,b,i};
                }
            }
        }
        m = cnt;
        for(int j=0;j<m;j++)plane[j] = np[j];
    }
}

int main(){
    cin>>n;
    for(int i=0;i<n;i++){
        cin>>q[i].x>>q[i].y>>q[i].z;
        q[i].skank();
    }
    get_conx_3d();

    double ans=0;
    for(int i=0;i<m;i++){
        ans+=plane[i].area();
    }
    printf("%lf",ans);
    return 0;
}

详情