1 条题解
-
0
C++ :
#include <cstdio> #include <cstring> #include <cstdlib> #include <cmath> #include <algorithm> #include <vector> using namespace std; #define sqr(x) ((x)*(x)) const double eps=1e-6; const double inf=8; const double pi=3.14159265358979323846; inline bool cmpDouble(const double &a,const double &b) {return fabs(a-b)<eps;} struct Tpoint { double x,y; Tpoint(){} Tpoint(double a,double b){x=a;y=b;} inline double norm(){ return sqrt( sqr(x)+sqr(y) ); } }; inline Tpoint operator +(const Tpoint &a,const Tpoint &b){return Tpoint(a.x+b.x,a.y+b.y);} inline Tpoint operator -(const Tpoint &a,const Tpoint &b){return Tpoint(a.x-b.x,a.y-b.y);} inline Tpoint operator *(const double &a,const Tpoint &b){return Tpoint(a*b.x,a*b.y);} inline Tpoint operator *(const Tpoint &b,const double &a){return Tpoint(a*b.x,a*b.y);} inline Tpoint operator /(const Tpoint &a,const double &b){return Tpoint(a.x/b,a.y/b);} inline bool operator <(const Tpoint &a,const Tpoint &b){return a.x+eps<b.x || fabs(a.x-b.x)<eps && a.y+eps<b.y;} inline bool operator ==(const Tpoint &a,const Tpoint &b){return fabs(a.x-b.x)<eps && fabs(a.y-b.y)<eps;} inline double det(const Tpoint &a,const Tpoint &b){return a.x*b.y-a.y*b.x;} inline double dot(const Tpoint &a,const Tpoint &b){return a.x*b.x+a.y*b.y;} struct P3 { double x,y,z; P3(){} P3(double a,double b,double c){x=a;y=b;z=c;} inline void read(){scanf("%lf%lf%lf",&x,&y,&z);} }; inline P3 operator +(const P3 &a,const P3 &b){return P3(a.x+b.x,a.y+b.y,a.z+b.z);} inline P3 operator -(const P3 &a,const P3 &b){return P3(a.x-b.x,a.y-b.y,a.z-b.z);} inline P3 operator *(const double &a,const P3 &b){return P3(a*b.x,a*b.y,a*b.z);} inline P3 operator *(const P3 &b,const double &a){return P3(a*b.x,a*b.y,a*b.z);} inline P3 operator /(const P3 &a,const double &b){return P3(a.x/b,a.y/b,a.z/b);} struct Tcir { double r; Tpoint o; Tcir(){} Tcir(Tpoint x,double y){o=x,r=y;} }; vector <Tcir> Circle; typedef vector <Tpoint> TP; vector <TP> Poly; struct Tinter { double x,y,Area,mid; int delta; Tinter(){} Tinter(double xx,double yy,double mm,int dd,double aa) { x=xx;y=yy;mid=mm; delta=dd;Area=aa; } }; inline bool operator <(const Tinter &a,const Tinter &b){return a.mid>b.mid+eps;} inline bool operator ==(const Tinter &a,const Tinter &b){return fabs(a.mid-b.mid)<eps;} vector <Tinter> inter; vector <double> bak; inline double dist(const Tpoint &a,const Tpoint &b) { return sqr(a.x-b.x)+sqr(a.y-b.y); } inline void Add(double x) { bak.push_back(x); } inline void CircleIntersectCircle(const Tcir &a,const Tcir &b) { double l=dist(a.o,b.o); double s=((a.r-b.r)*(a.r+b.r)/l+1)/2; double t=sqrt(-(l-sqr(a.r+b.r))*(l-sqr(a.r-b.r))/(l*l*4)); double ux=b.o.x-a.o.x,uy=b.o.y-a.o.y; double ix=a.o.x+s*ux+t*uy,iy=a.o.y+s*uy-t*ux; double jx=a.o.x+s*ux-t*uy,jy=a.o.y+s*uy+t*ux; Add(ix); Add(jx); } inline bool InLine(const Tpoint &a,const Tpoint &b,const Tpoint &c) { return fabs(det(b-a,c-a))<eps && dot(a-c,b-c)<eps; } inline void LineToLine(const Tpoint &a,const Tpoint &b,const Tpoint &c,const Tpoint &d) { double s1=det(c-a,b-a),s2=det(b-a,d-a); if (s1*s2<-eps) return; Tpoint e=c+(d-c)*s1/(s1+s2); if (InLine(a,b,e) && InLine(c,d,e)) { Add(e.x); } } inline void LineToCircle(const Tpoint &a,const Tpoint &b,const Tcir &c) { double h=fabs(det(c.o-a,b-a))/(b-a).norm(); if (h>c.r+eps) return; double lamda=dot(c.o-a,b-a); lamda/=dist(a,b); Tpoint x=a+(b-a)*lamda; double d=sqrt( sqr(c.r)-sqr(h) ); d/=(b-a).norm(); Tpoint e=x+(b-a)*d; Tpoint f=x-(b-a)*d; if (InLine(a,b,e)) Add(e.x); if (InLine(a,b,f)) Add(f.x); return; } inline void CircleToCircle() { for (int i=0;i<Circle.size();++i) { Add(Circle[i].o.x-Circle[i].r); Add(Circle[i].o.x+Circle[i].r); Add(Circle[i].o.x); for (int j=i+1;j<Circle.size();++j) if (dist(Circle[i].o,Circle[j].o)<=sqr(Circle[i].r+Circle[j].r)) if (dist(Circle[i].o,Circle[j].o)>=sqr(Circle[i].r-Circle[j].r)) CircleIntersectCircle(Circle[i],Circle[j]); } } inline void CircleToPoly() { for (int i=0;i<Circle.size();++i) for (int j=0;j<Poly.size();++j) for (int v=0;v<Poly[j].size();++v) LineToCircle(Poly[j][v],Poly[j][(v+1)%Poly[j].size()],Circle[i]); } inline void PolyToPoly() { for (int i=0;i<Poly.size();++i) for (int u=0;u<Poly[i].size();++u) Add(Poly[i][u].x); for (int i=0;i<Poly.size();++i) for (int j=i+1;j<Poly.size();++j) for (int u=0;u<Poly[i].size();++u) for (int v=0;v<Poly[j].size();++v) LineToLine(Poly[i][u],Poly[i][(u+1)%Poly[i].size()],Poly[j][v],Poly[j][(v+1)%Poly[j].size()]); } inline void Get(const Tcir &c,double x,double &l,double &r) { double y=fabs(c.o.x-x); double d=sqrt(fabs( sqr(c.r)-sqr(y) )); l=c.o.y+d; r=c.o.y-d; } inline double arcArea(const Tcir &a,double l,double x,double r,double y) { double len=sqrt(sqr(l-r) + sqr(x-y)); double d=sqrt(sqr(a.r)-sqr(len)/4.0); double angle=atan(len/2.0/d); return fabs(angle*sqr(a.r)-d*len/2.0); } inline void Get_Interval(const Tcir &a,double l,double r) { double L1,L2,R1,R2,M1,M2; Get(a,l,L1,L2); Get(a,r,R1,R2); Get(a,(l+r)/2.0,M1,M2); int D1=1,D2=-1; double A1=arcArea(a,l,L1,r,R1),A2=arcArea(a,l,L2,r,R2); inter.push_back( Tinter(L1,R1,M1,D1,A1) ); inter.push_back( Tinter(L2,R2,M2,D2,A2) ); } inline double calcSlice(double xl,double xr) { inter.clear(); double lmost=-inf,rmost=inf; for (int i=0;i<Poly.size();++i) { int cc=0; Tinter I[5]; for (int u=0;u<Poly[i].size();++u) { Tpoint x=Poly[i][u]; Tpoint y=Poly[i][(u+1)%Poly[i].size()]; double l=min(x.x,y.x),r=max(x.x,y.x); if (l<=xl+eps && xr<=r+eps) { if (fabs(l-r)<eps) continue; Tpoint d=y-x; Tpoint Left=x+d/d.x*(xl-x.x); Tpoint Right=x+d/d.x*(xr-x.x); Tpoint Mid=(Left+Right)/2; I[cc++]=Tinter(Left.y,Right.y,Mid.y,1,0); } } sort(I,I+cc); if (cc==2) { I[1].delta=-1; inter.push_back(I[0]); inter.push_back(I[1]); lmost=max(lmost,I[1].mid); rmost=min(rmost,I[0].mid); } } for (int i=0;i<Circle.size();++i) if (fabs(Circle[i].o.x-xl)<Circle[i].r+eps && fabs(Circle[i].o.x-xr)<Circle[i].r+eps) Get_Interval(Circle[i],xl,xr); if (!inter.size()) return 0; double ans=0; sort(inter.begin(),inter.end()); int cnt=0; for (int i=0;i<inter.size();++i) { if (cnt>0) { ans+=(fabs(inter[i-1].x-inter[i].x)+fabs(inter[i-1].y-inter[i].y))*(xr-xl)/2.0; ans+=inter[i-1].delta*inter[i-1].Area; ans-=inter[i].delta*inter[i].Area; } cnt+=inter[i].delta; } return ans; } #define maxn 105 int n,m; struct Poly4 { P3 p[4]; }a[maxn]; struct Sphere { P3 o; double r; inline void read(){o.read();scanf("%lf",&r);} }b[maxn]; inline bool equal(double a,double b) { return fabs(a-b)<eps; } inline bool InterSect(const Tpoint &a,const Tpoint &b,const Tpoint &c,const Tpoint &d) { double s1=det(b-a,c-a),s2=det(b-a,d-a); if (s1*s2>-eps) return false; s1=det(d-c,a-c),s2=det(d-c,b-c); if (s1*s2>-eps) return false; return true; } inline void ToHull(vector <Tpoint> &a) { sort(a.begin(),a.end()); int hull[10],len,limit=1; hull[len=1]=0; for (int i=1;i<4;++i) { while (len>limit && det(a[hull[len]]-a[hull[len-1]],a[i]-a[hull[len]])>=0) --len; hull[++len]=i; } limit=len; for (int i=2;i>=0;--i) { while (len>limit && det(a[hull[len]]-a[hull[len-1]],a[i]-a[hull[len]])>=0) --len; hull[++len]=i; } vector <Tpoint> b=a; a.resize(len-1); for (int i=0;i<len-1;++i) a[i]=b[hull[i+1]]; } inline double calcArea(double z) { Poly.clear(); Circle.clear(); bak.clear(); for (int i=0;i<n;++i) { vector <Tpoint> cross; for (int j=0;j<4;++j) for (int k=j+1;k<4;++k) if (!equal(a[i].p[j].z,a[i].p[k].z)) { double l=min(a[i].p[j].z,a[i].p[k].z),r=max(a[i].p[j].z,a[i].p[k].z); if (l<=z+eps && z<=r+eps) { P3 d=a[i].p[k]-a[i].p[j]; d=d/d.z; d=d*(z-a[i].p[j].z); d=d+a[i].p[j]; Tpoint x(d.x,d.y); cross.push_back(x); } } sort(cross.begin(),cross.end()); cross.erase(unique(cross.begin(),cross.end()),cross.end()); if (cross.size()>2) { if (cross.size()>4){puts("Too Many Points!!!");while (1);} if (cross.size()==4) ToHull(cross); Poly.push_back(cross); } } for (int i=0;i<m;++i) if (fabs(z-b[i].o.z)+eps<b[i].r) { Tpoint o(b[i].o.x,b[i].o.y); double r=sqrt( sqr(b[i].r)-sqr(z-b[i].o.z) ); Circle.push_back(Tcir(o,r)); } CircleToCircle(); CircleToPoly(); PolyToPoly(); sort(bak.begin(),bak.end()); double res=0; for (int i=0;i+1<bak.size();++i) if (fabs(bak[i+1]-bak[i])>eps) res+=calcSlice(bak[i],bak[i+1]); return res; } int main() { for (;scanf("%d%d",&n,&m)!=EOF && (n+m);) { for (int i=0;i<n;++i) for (int j=0;j<4;++j) a[i].p[j].read(); for (int i=0;i<m;++i) b[i].read(); double last=0,Ans=calcArea(-inf)+calcArea(inf); const int Block=4000; double h=(inf+inf)/(double)Block; for (int i=1;i<Block;i+=2) Ans+=4*calcArea(-inf+i*h); for (int i=2;i<Block;i+=2) Ans+=2*calcArea(-inf+i*h); Ans=Ans*h/3.0; printf("%.3f\n",Ans); } return 0; }
信息
- ID
- 191
- 时间
- 1000ms
- 内存
- 128MiB
- 难度
- (无)
- 标签
- 递交数
- 0
- 已通过
- 0
- 上传者