POJ3608

这题真神了，同样的代码用G++就WA用C++就AC……

题意是说，给两个凸多边形，求它们边界之间距离的最小值。做法类似于求凸包的旋转卡壳法，先取一个凸包的最高点和另一个凸包的最低点，画出方向相反的两条射线。然后同时逆时针旋转两条射线至其中一条与对应多边形的一条边相交，然后继续更新射线旋转。因为可以证明说最短距离是点到线的距离，所以不断更新就好了……

#include<stdio.h>#include<vector>#include<map>#include<set>#include<string.h>#include<math.h>#include<algorithm>#include<iostream>using namespace std;typedef long long ll;const double eps=1e-8;const double pi=acos(-1.0);int sgn(double x){if(fabs(x)<eps) return 0;return x>0?1:-1;}struct point{double x,y;point(double xt=0,double yt=0) :x(xt),y(yt){}};double dot(point a,point b){return a.x*b.x+a.y*b.y;}double cross(point a,point b){return a.x*b.y-a.y*b.x;}double length(point a){return sqrt(dot(a,a));}point operator +(point a,point b){return point(a.x+b.x,a.y+b.y);}point operator -(point a,point b){return point(a.x-b.x,a.y-b.y);}point operator *(point a,double b){return point(a.x*b,a.y*b);}point operator /(point a,double b){return point(a.x/b,a.y/b);}point unit(point a){return a/length(a)*100.0;}bool cmp(point a,point b){if(fabs(a.x-b.x)>eps) return a.x<b.x;return a.y<b.y;}bool operator==(point a,point b){return !sgn(a.x-b.x)&&!sgn(a.y-b.y);}struct line{point x,v;line(point xt=point(),point vt=point(1,0)) :x(xt),v(vt){}};line rotate(line te,double a){line re=te;re.v=point(te.v.x*cos(a)-te.v.y*sin(a),te.v.x*sin(a)+te.v.y*cos(a));return re;}double rotang(point a,point b){double ang=acos(dot(a,b)/length(a)/length(b));return fabs(ang);}int convexhull(point *p,int n,point *ch){sort(p,p+n,cmp);int m=0;for(int i=0;i<n;i++){while(m>1&&sgn(cross(ch[m-1]-ch[m-2],p[i]-ch[m-2]))<0) --m;ch[m++]=p[i];}int k=m;for(int i=n-2;i>=0;i--){while(m>k&&sgn(cross(ch[m-1]-ch[m-2],p[i]-ch[m-2]))<0) --m;ch[m++]=p[i];}if(n>1) --m;return m;}double ptos(point p,point a,point b){if(a==b) return length(p-a);point v1=b-a,v2=p-a,v3=p-b;if(sgn(dot(v1,v2))<0) return length(v2);if(sgn(dot(v1,v3))>0) return length(v3);return fabs(cross(v1,v2))/length(v1);}double stos(point a,point b,point c,point d){double tet=min(ptos(a,c,d),ptos(b,c,d));tet=min(tet,min(ptos(c,a,b),ptos(d,a,b)));return tet;}point pq[100005],qp[100005],p[100005],q[100005];int n,m,ni,mi;void solve(){int i,j;double ret=0;double mindist=length(p[ni]-q[mi]);line pt(p[ni],point(-100,0)),qt(q[mi],point(100,0));if(n<=2||m<=2){for(i=0;i<n;i++) for(j=0;j<m;j++) mindist=min(mindist,stos(p[(i+1)%n],p[i],q[(j+1)%m],q[j]));printf("%.5lf\n",mindist);return;}int nt=ni,mt=mi,t=0;while(ret<3*pi+eps){double ag1=rotang(pt.v,unit(p[(nt+1)%n]-p[nt])),ag2=rotang(qt.v,unit(q[(mt+1)%m]-q[mt]));//cout<<ag1<<" "<<ag2<<" !!"<<endl;if(fabs(ag1-ag2)<eps){double tet=stos(p[(nt+1)%n],p[nt],q[(mt+1)%m],q[mt]);mindist=min(tet,mindist);pt=line(p[(nt+1)%n],unit(p[(nt+1)%n]-p[nt]));qt=line(q[(mt+1)%m],unit(q[(mt+1)%m]-q[mt]));nt++;nt%=n;mt++;mt%=m;ret+=ag1;}else if(ag1>ag2){double tet=stos(p[(nt+1)%n],p[nt],q[(mt+1)%m],q[mt]);mindist=min(tet,mindist);//cout<<pt.v.x<<" "<<pt.v.y<<endl;pt=rotate(pt,ag2);//cout<<pt.v.x<<" "<<pt.v.y<<endl;qt=line(q[(mt+1)%m],unit(q[(mt+1)%m]-q[mt]));mt++;mt%=m;ret+=ag2;}else{double tet=stos(p[(nt+1)%n],p[nt],q[(mt+1)%m],q[mt]);mindist=min(tet,mindist);//cout<<qt.v.x<<" "<<qt.v.y<<endl;qt=rotate(qt,ag1);//cout<<qt.v.x<<" "<<qt.v.y<<endl;pt=line(p[(nt+1)%n],unit(p[(nt+1)%n]-p[nt]));nt++;nt%=n;ret+=ag1;}//printf("%lf %lf %lf %lf %lf %lf %lf %lf\n",pt.x.x,pt.x.y,pt.v.x,pt.v.y,qt.x.x,qt.x.y,qt.v.x,qt.v.y);}printf("%.5lf\n",mindist);}int main(){int i;while(scanf("%d %d",&n,&m)!=EOF&&n&&m){for(i=0,ni=0;i<n;i++) scanf("%lf %lf",&pq[i].x,&pq[i].y);for(i=0,mi=0;i<m;i++) scanf("%lf %lf",&qp[i].x,&qp[i].y);n=convexhull(pq,n,p);m=convexhull(qp,m,q);for(i=0;i<n;i++) if(p[i].y-p[ni].y>eps) ni=i;for(i=0;i<m;i++) if(q[i].y-q[mi].y<-eps) mi=i;solve();}return 0;}