Coverthe Polygon with Your Disk
题意：给定一个最多10个点的凸多边形，和一个圆，求圆与多边形相交的最大面积。
根据圆坐标的不同，面积交为一个二维凸函数。那么就有多种方法。
梯度下降，比较可靠的方法，过程相当于爬山。具体做法，对当前所在点求x,y的偏导，走的方向就是两个偏导组成的向量。这个向量称为梯度。梯度方向是最大的方向导数的方向，梯度的模是最大方向导数值。如果梯度为负，则往反方向走。
模拟退火，走的时候试探k个角度，取最大值。
这题还可以三分套三分，但是好像不方便处理。
#include<bits/stdc++.h> using namespace std; const double eps=1e-8; const double PI=acos(-1); inline int sgn(double x){return x<-eps?-1:x>eps;} inline double sqr(double x){return x*x;}//square struct Point{ void read(){scanf("%lf%lf",&x,&y);} void out(){printf("Point(%.7f,%.7f)",x,y);} Point(double x=0,double y=0):x(x),y(y){} bool operator ==(const Point &B)const{return !sgn(x-B.x)&&!sgn(y-B.y);} bool operator <(const Point &B)const{return sgn(x-B.x)<0||(sgn(x-B.x)==0&&sgn(y-B.y)<0);} Point operator +(Point b)const{return Point(x+b.x,y+b.y);} Point operator -(Point b)const{return Point(x-b.x,y-b.y);} Point operator *(double p)const{return Point(x*p,y*p);} Point operator /(double p)const{return Point(x/p,y/p);} double x,y; }; typedef Point Vector; //atan2(y,x)求向量极角(-PI,PI],x,y不同时等于0 double Dot(Vector a,Vector b){return a.x*b.x+a.y*b.y;} double Cross(Vector a,Vector b){return a.x*b.y-a.y*b.x;} double Length(Vector a){return sqrt(Dot(a,a));} double Dis(Point A,Point B){return sqrt(sqr(A.x-B.x)+sqr(A.y-B.y));} double Angle(Vector a,Vector b){return acos(Dot(a,b)/Length(a)/Length(b));}//[0~PI] Vector Rotate(Vector a,double rad){return Vector(a.x*cos(rad)-a.y*sin(rad),a.x*sin(rad)+a.y*cos(rad));}//逆时针旋转 Vector Normal(Vector a){double L=Length(a);return Vector(-a.y/L,a.x/L);}//单位法向量 Vector Trunc(Vector a,double L){return a/Length(a)*L;}//截成长度为L的向量 struct Line{ void read(){A.read();B.read();ang=atan2(B.y-A.y,B.x-A.x);} void out(){printf("Line(Point(%.7f,%.7f),Point(%.7f,%.7f))",A.x,A.y,B.x,B.y);} Line()=default; Line(Point A,Point B):A(A),B(B),ang(atan2(B.y-A.y,B.x-A.x)){} bool operator<(const Line &L)const{return ang<L.ang;} Point A,B; double ang; }; typedef Line Seg; void StdLine(Line L,double &A,double &B,double &C){ A=L.B.y-L.A.y; B=L.A.x-L.B.x; C=L.B.x*L.A.y-L.A.x*L.B.y; } Point IntersecLL(Point P,Vector v,Point Q,Vector w){//直线交点，参数方程 double t=Cross(w,P-Q)/Cross(v,w); return P+v*t; } Point IntersecLL(Line L1,Line L2){//直线交点，两点式 return IntersecLL(L1.A,L1.B-L1.A,L2.A,L2.B-L2.A); } Line PerpenLine(Line L){//Perpendicular中垂线 return Line(Point((L.A+L.B)/2),Point((L.A+L.B)/2)+Normal(L.B-L.A)); } double DisPL(Point P,Line L){ return fabs(Cross(L.B-L.A,P-L.A))/Dis(L.A,L.B); } double DisPS(Point P,Seg S){ if(S.A==S.B)return Dis(S.A,P); Vector v1=S.B-S.A,v2=P-S.A,v3=P-S.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); } Point ProjecPL(Point P,Line L){//投影 Vector v=L.B-L.A; return L.A+v*(Dot(v,P-L.A)/Dot(v,v)); } //1:left 0:on -1:right int RelationPL(Point P,Line L){ return sgn(Cross(L.B-L.A,P-L.A)); } //0:不在线段上,1:在线段非端点处,2:在端点上 int RelationPS(Point P,Seg S){ if(P==S.A||P==S.B)return 2; return !sgn(Cross(S.A-P,S.B-P))&&Dot(S.A-P,S.B-P)<0; } //0:不相交,1:规范相交,2:交于端点或重合 int RelationSS(Seg S1,Seg S2){ if(!(min(S1.A.x,S1.B.x)<=max(S2.A.x,S2.B.x)&&min(S2.A.x,S2.B.x)<=max(S1.A.x,S1.B.x)&&//快速排斥 min(S1.A.y,S1.B.y)<=max(S2.A.y,S2.B.y)&&min(S2.A.y,S2.B.y)<=max(S1.A.y,S1.B.y)))return 0; int c1=sgn(Cross(S1.B-S1.A,S2.A-S1.A)),c2=sgn(Cross(S1.B-S1.A,S2.B-S1.A)), c3=sgn(Cross(S2.B-S2.A,S1.A-S2.A)),c4=sgn(Cross(S2.B-S2.A,S1.B-S2.A)); if((c1^c2)==-2&&(c3^c4)==-2)return 1; if(c1*c2<=0&&c3*c4<=0)return 2; return 0; } struct Circle{ Circle()=default; Circle(Point C,double r):C(C),r(r){} void read(){C.read();scanf("%lf",&r);} Point point(double a){return Point(C.x+cos(a)*r,C.y+sin(a)*r);} Point C; double r; }; //0:在圆内1:在圆上2:在圆外 int RelationPC(Point P,Circle C){ return sgn(Dis(C.C,P)-C.r)+1; } //返回直线与圆交点个数 int RelationLC(Line L,Circle C){ return sgn(C.r-DisPL(C.C,L))+1; } //-1:重合0:内含1:内切2:相交3:外切4:相离 int RelationCC(Circle C1,Circle C2){ if(C1.C==C2.C&&C1.r==C2.r)return -1; double d=Dis(C1.C,C2.C); int cmp=sgn(d-fabs(C1.r-C2.r)); if(cmp<1)return cmp+1; cmp=sgn(d-C1.r-C2.r); return cmp+3; } //得到直线与圆的交点,返回交点个数,已测试 int IntersecLC(Line L,Circle C,Point &P1,Point &P2){ Point P=ProjecPL(C.C,L); double d=sqr(C.r)-sqr(Dis(C.C,P)); if(sgn(d)<0)return 0; if(sgn(d)==0){P1=P2=P;return 1;} d=sqrt(d);//防止相切时直线在圆外 Vector v=(L.B-L.A)/Dis(L.A,L.B); P1=P+v*d;P2=P-v*d; return P1==P2?1:2; } //得到两圆交点,返回两圆交点数,-1:重合,已测试 int IntersecCC(Circle C1,Circle C2,Point &P1,Point &P2){ if(C1.C==C2.C&&C1.r==C2.r)return -1; double d=Dis(C1.C,C2.C); if(sgn(C1.r+C2.r-d)<0||sgn(d-fabs(C1.r-C2.r))<0)return 0; double h=(d+(sqr(C1.r)-sqr(C2.r))/d)/2; Vector v=C2.C-C1.C; Point P=C1.C+Trunc(v,h); h=sqrt(fabs(sqr(C1.r)-sqr(h)));//fabs防止相切变nan P1=P+Normal(v)*h; P2=P-Normal(v)*h; return (sgn(C1.r+C2.r-d)==0||sgn(d-fabs(C1.r-C2.r))==0||P1==P2)?1:2; } //点向圆作切线,得到切点,返回切点个数,已测试 int TangentPC(Point P,Circle C,Point &P1,Point &P2){ int rel=RelationPC(P,C); if(rel==0)return 0; if(rel==1){P1=P2=P;return 1;} double d=Dis(C.C,P),l=C.r*C.r/d,h=sqrt(C.r*C.r-l*l); Vector v=P-C.C; P1=C.C+Trunc(v,l)+Normal(v)*h; P2=C.C+Trunc(v,l)-Normal(v)*h; return P1==P2?1:2; } //点向圆作切线,得到切线,返回切线条数 int TangentPC(Point P,Circle C,Line &L1,Line &L2){ int rel=RelationPC(P,C); if(rel==0)return 0; if(rel==1){ L1=L2=Line(P,P+Normal(P-C.C)); return 1; } double d=Dis(C.C,P),l=C.r*C.r/d,h=sqrt(C.r*C.r-l*l); Vector v=P-C.C; L1=Line(P,C.C+Trunc(v,l)+Normal(v)*h); L2=Line(P,C.C+Trunc(v,l)-Normal(v)*h); return 2; } //得到两圆公切点,返回切点对数,-1:两圆重合,已测试 int TangentCC(Circle C1,Circle C2,Point *p1,Point *p2){ int rel=RelationCC(C1,C2); if(rel<1)return rel; if(C1.r<C2.r){swap(C1,C2);swap(p1,p2);} int cnt=0; if(rel==1){//内切 p1[cnt]=p2[cnt]=C1.C+Trunc(C2.C-C1.C,C1.r); return ++cnt; } double d=Dis(C1.C,C2.C),ang=acos((C1.r-C2.r)/d); double base=atan2(C2.C.y-C1.C.y,C2.C.x-C1.C.x); p1[cnt]=C1.point(base+ang);p2[cnt++]=C2.point(base+ang); p1[cnt]=C1.point(base-ang);p2[cnt++]=C2.point(base-ang); if(rel==3){p1[cnt]=p2[cnt]=C1.C+Trunc(C2.C-C1.C,C1.r);++cnt;}//外切 if(rel!=4)return cnt; ang=acos((C1.r+C2.r)/d); p1[cnt]=C1.point(base+ang);p2[cnt++]=C2.point(base+ang+PI); p1[cnt]=C1.point(base-ang);p2[cnt++]=C2.point(base-ang+PI); return cnt; } //外接圆,已测试 Circle CircumscribedCircle(Point A,Point B,Point C){ Line L1=PerpenLine(Line(A,B)); Line L2=PerpenLine(Line(B,C)); Point P=IntersecLL(L1,L2); return Circle(P,DisPL(P,L1)); } //内切圆,已测试 Circle InscribedCircle(Point A,Point B,Point C){ double a=Dis(B,C),b=Dis(A,C),c=Dis(A,B); Point P((A*a+B*b+C*c)/(a+b+c)); return Circle(P,DisPL(P,Line(A,B))); } typedef vector<Point> Poly; double Area2(Point A,Point B,Point C){return Cross(B-A,C-A);}//三角形有向面积的二倍 //简单多边形面积，已测试 double AreaPoly(Poly &p){ int n=p.size(); double area=0; for(int i=1;i<n-1;++i) area+=Cross(p[i]-p[0],p[i+1]-p[0]); return area/2; } //圆与三角形相交的面积，要求三角形一个点是圆心，hdu2892测 double AreaIntersecTriangleC(Point A,Point B,Circle C){ if(sgn(Cross(A-C.C,B-C.C))==0)return 0; int rela=RelationPC(A,C),relb=RelationPC(B,C); if(rela<=1&&relb<=1)return fabs(Cross(A-C.C,B-C.C))/2;//两个点都在圆内 if(rela>relb){swap(A,B);swap(rela,relb);} Point a,b; if(IntersecLC(Line(A,B),C,a,b)<=1)return Angle(A-C.C,B-C.C)/2*C.r*C.r;//小于两个交点 if(sgn(Cross(a-C.C,b-C.C))*sgn(Cross(a-C.C,B-C.C))<0)swap(a,b); if(rela<=1)return Angle(B-C.C,b-C.C)/2*C.r*C.r+fabs(Cross(b-C.C,A-C.C))/2;//两个点一内一外 if(sgn(Cross(a-C.C,A-C.C))*sgn(Cross(b-C.C,B-C.C))>0) return Angle(A-C.C,B-C.C)/2*C.r*C.r;//两个交点在一侧，两个点都在圆外 return (Angle(A-C.C,a-C.C)+Angle(B-C.C,b-C.C))/2*C.r*C.r+fabs(Cross(a-C.C,b-C.C))/2;//两个点都在圆外 } //简单多边形与圆相交的面积,拆成三角形计算,hdu2892测 double AreaIntersecPolyC(Poly &p,Circle C){ double area=0;p.push_back(p[0]); for(int i=0;i<p.size()-1;++i) area+=AreaIntersecTriangleC(p[i],p[i+1],C)*sgn(Cross(p[i]-C.C,p[i+1]-C.C)); p.pop_back(); return fabs(area); } //点与简单多边形位置关系-1:在边界，0:外部，1:内部,O(N),已测试 int RelationPPoly(Point A,Poly &p){ int n=p.size(),wn=0;//winding number for(int i=0;i<n;++i){ if(RelationPS(A,Seg(p[i],p[(i+1)%n]))!=0)return -1; int k=sgn(Cross(p[(i+1)%n]-p[i],A-p[i]));//-1:右侧1:左侧 int d1=sgn(p[i].y-A.y),d2=sgn(p[(i+1)%n].y-A.y); if(k>0&&d1<=0&&d2>0)++wn; if(k<0&&d2<=0&&d1>0)--wn; } return wn!=0;//wn==0说明在外部 } //点与凸多边形位置关系O(logN)0:outside,1:in or on,已测试 int RelationPConvex(Point A,Poly &p){ int n=p.size()-1; if(sgn(Cross(p[1]-p[0],A-p[0]))<0||sgn(Cross(p[n]-p[0],A-p[0]))>0) return 0; int l=1,r=n; while(r-l>1){ int mid=(l+r)>>1; if(Cross(p[mid]-p[0],A-p[0])>0)l=mid; else r=mid; } return sgn(Cross(p[r]-p[l],A-p[l]))>=0; } //求凸包，已测试 Poly ConvexHull(Poly &p){ int n=p.size(),m=0; Poly ch(n+1); sort(p.begin(),p.end()); for(int i=0;i<n;++i){ while(m>1&&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&&Cross(ch[m-1]-ch[m-2],p[i]-ch[m-2])<=0)--m; ch[m++]=p[i]; } ch.resize(m>1?m-1:m);//起点包含两次 return ch; } //直线切割凸多边形，退化返回点或线段，已测试 Poly LineCutPoly(Line L,Poly p){ Poly newpoly; int n=p.size(); Point A=L.A,B=L.B,C,D,ip; for(int i=0;i<n;++i){ C=p[i];D=p[(i+1)%n]; if(sgn(Cross(B-A,C-A))>=0)newpoly.push_back(C);//在左侧 if(!sgn(Cross(B-A,D-C)))continue;//不平行 ip=IntersecLL(Line(A,B),Line(C,D)); if(RelationPS(ip,Seg(C,D))==1)//交点在线段上去且不在端点 newpoly.push_back(ip); } return newpoly; } //半平面交,需要保证得到的不是无穷区域，无穷区域加一个大矩形。 //半平面交不存在返回空集。无穷大返回空集或求最后的交点时出错。 Poly Halfplane(Line *L,int n){ sort(L,L+n); Point *p=new Point[n]; Line *q=new Line[n]; int s=0,t=0; q[0]=L[0]; for(int i=1;i<n;++i){ while(s<t&&RelationPL(p[t-1],L[i])!=1)--t; while(s<t&&RelationPL(p[s],L[i])!=1)++s; q[++t]=L[i]; if(sgn(q[t].ang-q[t-1].ang)==0)if(RelationPL(L[i].A,q[--t])==1)q[t]=L[i];//平行取内侧 if(s<t)p[t-1]=IntersecLL(q[t],q[t-1]); } while(s<t&&RelationPL(p[t-1],q[s])!=1)--t;//最后一个交点与第一条线 //while(s<t&&RelationPL(p[s],q[t])!=1)++s; if(t-s>1)p[t++]=IntersecLL(q[s],q[t]);//首尾交点，不存在则无穷 else s=t=0; Poly poly(p+s,p+t); delete[] p;delete[] q; return poly; } //凸包闵可夫斯基和，要求pa,pb都是凸包且起点均为最左下 Poly Minkowski(Poly &pa,Poly &pb){ int na=pa.size(),nb=pb.size(),tot=0,ia=0,ib=0; Vector *va=new Vector[na],*vb=new Vector[nb]; for(int i=0;i<na;++i)va[i]=pa[(i+1)%na]-pa[i]; for(int i=0;i<nb;++i)vb[i]=pb[(i+1)%nb]-pb[i]; Poly p(na+nb+1); p[tot++]=pa[0]+pb[0]; while(ia<na&&ib<nb)p[tot]=p[tot-1]+(Cross(va[ia],vb[ib])>0?va[ia++]:vb[ib++]),++tot; while(ia<na)p[tot]=p[tot-1]+va[ia++],++tot; while(ib<nb)p[tot]=p[tot-1]+vb[ib++],++tot; delete[] va;delete[] vb; //删除在边上的点 { int x=1; for(int i=2;i<tot;++i) if(sgn(Cross(p[x]-p[x-1],p[i]-p[i-1]))==0)p[x]=p[i]; else p[++x]=p[i]; tot=x+1; } p.resize(tot-1);//删除首尾重复 return p; } double r; double calc(Point o,Poly &p){ return AreaIntersecPolyC(p,Circle(o,r)); } int main(){ int n; scanf("%d%lf",&n,&r); Poly p(n); Point o; for(int i=0;i<n;++i){ p[i].read(); o=o+p[i]; } o=o/n; double T=100,D=0.9; double EPS=1e-6; while(T>EPS){ double f=calc(o,p); double fx=calc(o+Point(EPS,0),p); double fy=calc(o+Point(0,EPS),p); Point o2=o+Trunc(Vector((fx-f)/EPS,(fy-f)/EPS),T); if(calc(o2,p)>calc(o,p))o=o2; else T*=D; } printf("%.6f",calc(o,p)); return 0; } /* 4 4 0 0 6 0 6 6 0 6 35.759506 3 1 0 0 2 1 1 3 2.113100 3 1 0 0 100 1 99 1 0.019798 4 1 0 0 100 10 100 12 0 1 3.137569 10 10 0 0 10 0 20 1 30 3 40 6 50 10 60 15 70 21 80 28 90 36 177.728187 */
#include<bits/stdc++.h> using namespace std; const double eps=1e-8; const double PI=acos(-1); inline int sgn(double x){return x<-eps?-1:x>eps;} inline double sqr(double x){return x*x;}//square struct Point{ void read(){scanf("%lf%lf",&x,&y);} void out(){printf("Point(%.7f,%.7f)",x,y);} Point(double x=0,double y=0):x(x),y(y){} bool operator ==(const Point &B)const{return !sgn(x-B.x)&&!sgn(y-B.y);} bool operator <(const Point &B)const{return sgn(x-B.x)<0||(sgn(x-B.x)==0&&sgn(y-B.y)<0);} Point operator +(Point b)const{return Point(x+b.x,y+b.y);} Point operator -(Point b)const{return Point(x-b.x,y-b.y);} Point operator *(double p)const{return Point(x*p,y*p);} Point operator /(double p)const{return Point(x/p,y/p);} double x,y; }; typedef Point Vector; //atan2(y,x)求向量极角(-PI,PI],x,y不同时等于0 double Dot(Vector a,Vector b){return a.x*b.x+a.y*b.y;} double Cross(Vector a,Vector b){return a.x*b.y-a.y*b.x;} double Length(Vector a){return sqrt(Dot(a,a));} double Dis(Point A,Point B){return sqrt(sqr(A.x-B.x)+sqr(A.y-B.y));} double Angle(Vector a,Vector b){return acos(Dot(a,b)/Length(a)/Length(b));}//[0~PI] Vector Rotate(Vector a,double rad){return Vector(a.x*cos(rad)-a.y*sin(rad),a.x*sin(rad)+a.y*cos(rad));}//逆时针旋转 Vector Normal(Vector a){double L=Length(a);return Vector(-a.y/L,a.x/L);}//单位法向量 Vector Trunc(Vector a,double L){return a/Length(a)*L;}//截成长度为L的向量 struct Line{ void read(){A.read();B.read();ang=atan2(B.y-A.y,B.x-A.x);} void out(){printf("Line(Point(%.7f,%.7f),Point(%.7f,%.7f))",A.x,A.y,B.x,B.y);} Line()=default; Line(Point A,Point B):A(A),B(B),ang(atan2(B.y-A.y,B.x-A.x)){} bool operator<(const Line &L)const{return ang<L.ang;} Point A,B; double ang; }; typedef Line Seg; void StdLine(Line L,double &A,double &B,double &C){ A=L.B.y-L.A.y; B=L.A.x-L.B.x; C=L.B.x*L.A.y-L.A.x*L.B.y; } Point IntersecLL(Point P,Vector v,Point Q,Vector w){//直线交点，参数方程 double t=Cross(w,P-Q)/Cross(v,w); return P+v*t; } Point IntersecLL(Line L1,Line L2){//直线交点，两点式 return IntersecLL(L1.A,L1.B-L1.A,L2.A,L2.B-L2.A); } Line PerpenLine(Line L){//Perpendicular中垂线 return Line(Point((L.A+L.B)/2),Point((L.A+L.B)/2)+Normal(L.B-L.A)); } double DisPL(Point P,Line L){ return fabs(Cross(L.B-L.A,P-L.A))/Dis(L.A,L.B); } double DisPS(Point P,Seg S){ if(S.A==S.B)return Dis(S.A,P); Vector v1=S.B-S.A,v2=P-S.A,v3=P-S.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); } Point ProjecPL(Point P,Line L){//投影 Vector v=L.B-L.A; return L.A+v*(Dot(v,P-L.A)/Dot(v,v)); } //1:left 0:on -1:right int RelationPL(Point P,Line L){ return sgn(Cross(L.B-L.A,P-L.A)); } //0:不在线段上,1:在线段非端点处,2:在端点上 int RelationPS(Point P,Seg S){ if(P==S.A||P==S.B)return 2; return !sgn(Cross(S.A-P,S.B-P))&&Dot(S.A-P,S.B-P)<0; } //0:不相交,1:规范相交,2:交于端点或重合 int RelationSS(Seg S1,Seg S2){ if(!(min(S1.A.x,S1.B.x)<=max(S2.A.x,S2.B.x)&&min(S2.A.x,S2.B.x)<=max(S1.A.x,S1.B.x)&&//快速排斥 min(S1.A.y,S1.B.y)<=max(S2.A.y,S2.B.y)&&min(S2.A.y,S2.B.y)<=max(S1.A.y,S1.B.y)))return 0; int c1=sgn(Cross(S1.B-S1.A,S2.A-S1.A)),c2=sgn(Cross(S1.B-S1.A,S2.B-S1.A)), c3=sgn(Cross(S2.B-S2.A,S1.A-S2.A)),c4=sgn(Cross(S2.B-S2.A,S1.B-S2.A)); if((c1^c2)==-2&&(c3^c4)==-2)return 1; if(c1*c2<=0&&c3*c4<=0)return 2; return 0; } struct Circle{ Circle()=default; Circle(Point C,double r):C(C),r(r){} void read(){C.read();scanf("%lf",&r);} Point point(double a){return Point(C.x+cos(a)*r,C.y+sin(a)*r);} Point C; double r; }; //0:在圆内1:在圆上2:在圆外 int RelationPC(Point P,Circle C){ return sgn(Dis(C.C,P)-C.r)+1; } //返回直线与圆交点个数 int RelationLC(Line L,Circle C){ return sgn(C.r-DisPL(C.C,L))+1; } //-1:重合0:内含1:内切2:相交3:外切4:相离 int RelationCC(Circle C1,Circle C2){ if(C1.C==C2.C&&C1.r==C2.r)return -1; double d=Dis(C1.C,C2.C); int cmp=sgn(d-fabs(C1.r-C2.r)); if(cmp<1)return cmp+1; cmp=sgn(d-C1.r-C2.r); return cmp+3; } //得到直线与圆的交点,返回交点个数,已测试 int IntersecLC(Line L,Circle C,Point &P1,Point &P2){ Point P=ProjecPL(C.C,L); double d=sqr(C.r)-sqr(Dis(C.C,P)); if(sgn(d)<0)return 0; if(sgn(d)==0){P1=P2=P;return 1;} d=sqrt(d);//防止相切时直线在圆外 Vector v=(L.B-L.A)/Dis(L.A,L.B); P1=P+v*d;P2=P-v*d; return P1==P2?1:2; } //得到两圆交点,返回两圆交点数,-1:重合,已测试 int IntersecCC(Circle C1,Circle C2,Point &P1,Point &P2){ if(C1.C==C2.C&&C1.r==C2.r)return -1; double d=Dis(C1.C,C2.C); if(sgn(C1.r+C2.r-d)<0||sgn(d-fabs(C1.r-C2.r))<0)return 0; double h=(d+(sqr(C1.r)-sqr(C2.r))/d)/2; Vector v=C2.C-C1.C; Point P=C1.C+Trunc(v,h); h=sqrt(fabs(sqr(C1.r)-sqr(h)));//fabs防止相切变nan P1=P+Normal(v)*h; P2=P-Normal(v)*h; return (sgn(C1.r+C2.r-d)==0||sgn(d-fabs(C1.r-C2.r))==0||P1==P2)?1:2; } //点向圆作切线,得到切点,返回切点个数,已测试 int TangentPC(Point P,Circle C,Point &P1,Point &P2){ int rel=RelationPC(P,C); if(rel==0)return 0; if(rel==1){P1=P2=P;return 1;} double d=Dis(C.C,P),l=C.r*C.r/d,h=sqrt(C.r*C.r-l*l); Vector v=P-C.C; P1=C.C+Trunc(v,l)+Normal(v)*h; P2=C.C+Trunc(v,l)-Normal(v)*h; return P1==P2?1:2; } //点向圆作切线,得到切线,返回切线条数 int TangentPC(Point P,Circle C,Line &L1,Line &L2){ int rel=RelationPC(P,C); if(rel==0)return 0; if(rel==1){ L1=L2=Line(P,P+Normal(P-C.C)); return 1; } double d=Dis(C.C,P),l=C.r*C.r/d,h=sqrt(C.r*C.r-l*l); Vector v=P-C.C; L1=Line(P,C.C+Trunc(v,l)+Normal(v)*h); L2=Line(P,C.C+Trunc(v,l)-Normal(v)*h); return 2; } //得到两圆公切点,返回切点对数,-1:两圆重合,已测试 int TangentCC(Circle C1,Circle C2,Point *p1,Point *p2){ int rel=RelationCC(C1,C2); if(rel<1)return rel; if(C1.r<C2.r){swap(C1,C2);swap(p1,p2);} int cnt=0; if(rel==1){//内切 p1[cnt]=p2[cnt]=C1.C+Trunc(C2.C-C1.C,C1.r); return ++cnt; } double d=Dis(C1.C,C2.C),ang=acos((C1.r-C2.r)/d); double base=atan2(C2.C.y-C1.C.y,C2.C.x-C1.C.x); p1[cnt]=C1.point(base+ang);p2[cnt++]=C2.point(base+ang); p1[cnt]=C1.point(base-ang);p2[cnt++]=C2.point(base-ang); if(rel==3){p1[cnt]=p2[cnt]=C1.C+Trunc(C2.C-C1.C,C1.r);++cnt;}//外切 if(rel!=4)return cnt; ang=acos((C1.r+C2.r)/d); p1[cnt]=C1.point(base+ang);p2[cnt++]=C2.point(base+ang+PI); p1[cnt]=C1.point(base-ang);p2[cnt++]=C2.point(base-ang+PI); return cnt; } //外接圆,已测试 Circle CircumscribedCircle(Point A,Point B,Point C){ Line L1=PerpenLine(Line(A,B)); Line L2=PerpenLine(Line(B,C)); Point P=IntersecLL(L1,L2); return Circle(P,DisPL(P,L1)); } //内切圆,已测试 Circle InscribedCircle(Point A,Point B,Point C){ double a=Dis(B,C),b=Dis(A,C),c=Dis(A,B); Point P((A*a+B*b+C*c)/(a+b+c)); return Circle(P,DisPL(P,Line(A,B))); } typedef vector<Point> Poly; double Area2(Point A,Point B,Point C){return Cross(B-A,C-A);}//三角形有向面积的二倍 //简单多边形面积，已测试 double AreaPoly(Poly &p){ int n=p.size(); double area=0; for(int i=1;i<n-1;++i) area+=Cross(p[i]-p[0],p[i+1]-p[0]); return area/2; } //圆与三角形相交的面积，要求三角形一个点是圆心，hdu2892测 double AreaIntersecTriangleC(Point A,Point B,Circle C){ if(sgn(Cross(A-C.C,B-C.C))==0)return 0; int rela=RelationPC(A,C),relb=RelationPC(B,C); if(rela<=1&&relb<=1)return fabs(Cross(A-C.C,B-C.C))/2;//两个点都在圆内 if(rela>relb){swap(A,B);swap(rela,relb);} Point a,b; if(IntersecLC(Line(A,B),C,a,b)<=1)return Angle(A-C.C,B-C.C)/2*C.r*C.r;//小于两个交点 if(sgn(Cross(a-C.C,b-C.C))*sgn(Cross(a-C.C,B-C.C))<0)swap(a,b); if(rela<=1)return Angle(B-C.C,b-C.C)/2*C.r*C.r+fabs(Cross(b-C.C,A-C.C))/2;//两个点一内一外 if(sgn(Cross(a-C.C,A-C.C))*sgn(Cross(b-C.C,B-C.C))>0) return Angle(A-C.C,B-C.C)/2*C.r*C.r;//两个交点在一侧，两个点都在圆外 return (Angle(A-C.C,a-C.C)+Angle(B-C.C,b-C.C))/2*C.r*C.r+fabs(Cross(a-C.C,b-C.C))/2;//两个点都在圆外 } //简单多边形与圆相交的面积,拆成三角形计算,hdu2892测 double AreaIntersecPolyC(Poly &p,Circle C){ double area=0;p.push_back(p[0]); for(int i=0;i<p.size()-1;++i) area+=AreaIntersecTriangleC(p[i],p[i+1],C)*sgn(Cross(p[i]-C.C,p[i+1]-C.C)); p.pop_back(); return fabs(area); } //点与简单多边形位置关系-1:在边界，0:外部，1:内部,O(N),已测试 int RelationPPoly(Point A,Poly &p){ int n=p.size(),wn=0;//winding number for(int i=0;i<n;++i){ if(RelationPS(A,Seg(p[i],p[(i+1)%n]))!=0)return -1; int k=sgn(Cross(p[(i+1)%n]-p[i],A-p[i]));//-1:右侧1:左侧 int d1=sgn(p[i].y-A.y),d2=sgn(p[(i+1)%n].y-A.y); if(k>0&&d1<=0&&d2>0)++wn; if(k<0&&d2<=0&&d1>0)--wn; } return wn!=0;//wn==0说明在外部 } //点与凸多边形位置关系O(logN)0:outside,1:in or on,已测试 int RelationPConvex(Point A,Poly &p){ int n=p.size()-1; if(sgn(Cross(p[1]-p[0],A-p[0]))<0||sgn(Cross(p[n]-p[0],A-p[0]))>0) return 0; int l=1,r=n; while(r-l>1){ int mid=(l+r)>>1; if(Cross(p[mid]-p[0],A-p[0])>0)l=mid; else r=mid; } return sgn(Cross(p[r]-p[l],A-p[l]))>=0; } //求凸包，已测试 Poly ConvexHull(Poly &p){ int n=p.size(),m=0; Poly ch(n+1); sort(p.begin(),p.end()); for(int i=0;i<n;++i){ while(m>1&&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&&Cross(ch[m-1]-ch[m-2],p[i]-ch[m-2])<=0)--m; ch[m++]=p[i]; } ch.resize(m>1?m-1:m);//起点包含两次 return ch; } //直线切割凸多边形，退化返回点或线段，已测试 Poly LineCutPoly(Line L,Poly p){ Poly newpoly; int n=p.size(); Point A=L.A,B=L.B,C,D,ip; for(int i=0;i<n;++i){ C=p[i];D=p[(i+1)%n]; if(sgn(Cross(B-A,C-A))>=0)newpoly.push_back(C);//在左侧 if(!sgn(Cross(B-A,D-C)))continue;//不平行 ip=IntersecLL(Line(A,B),Line(C,D)); if(RelationPS(ip,Seg(C,D))==1)//交点在线段上去且不在端点 newpoly.push_back(ip); } return newpoly; } //半平面交,需要保证得到的不是无穷区域，无穷区域加一个大矩形。 //半平面交不存在返回空集。无穷大返回空集或求最后的交点时出错。 Poly Halfplane(Line *L,int n){ sort(L,L+n); Point *p=new Point[n]; Line *q=new Line[n]; int s=0,t=0; q[0]=L[0]; for(int i=1;i<n;++i){ while(s<t&&RelationPL(p[t-1],L[i])!=1)--t; while(s<t&&RelationPL(p[s],L[i])!=1)++s; q[++t]=L[i]; if(sgn(q[t].ang-q[t-1].ang)==0)if(RelationPL(L[i].A,q[--t])==1)q[t]=L[i];//平行取内侧 if(s<t)p[t-1]=IntersecLL(q[t],q[t-1]); } while(s<t&&RelationPL(p[t-1],q[s])!=1)--t;//最后一个交点与第一条线 //while(s<t&&RelationPL(p[s],q[t])!=1)++s; if(t-s>1)p[t++]=IntersecLL(q[s],q[t]);//首尾交点，不存在则无穷 else s=t=0; Poly poly(p+s,p+t); delete[] p;delete[] q; return poly; } //凸包闵可夫斯基和，要求pa,pb都是凸包且起点均为最左下 Poly Minkowski(Poly &pa,Poly &pb){ int na=pa.size(),nb=pb.size(),tot=0,ia=0,ib=0; Vector *va=new Vector[na],*vb=new Vector[nb]; for(int i=0;i<na;++i)va[i]=pa[(i+1)%na]-pa[i]; for(int i=0;i<nb;++i)vb[i]=pb[(i+1)%nb]-pb[i]; Poly p(na+nb+1); p[tot++]=pa[0]+pb[0]; while(ia<na&&ib<nb)p[tot]=p[tot-1]+(Cross(va[ia],vb[ib])>0?va[ia++]:vb[ib++]),++tot; while(ia<na)p[tot]=p[tot-1]+va[ia++],++tot; while(ib<nb)p[tot]=p[tot-1]+vb[ib++],++tot; delete[] va;delete[] vb; //删除在边上的点 { int x=1; for(int i=2;i<tot;++i) if(sgn(Cross(p[x]-p[x-1],p[i]-p[i-1]))==0)p[x]=p[i]; else p[++x]=p[i]; tot=x+1; } p.resize(tot-1);//删除首尾重复 return p; } double r; double calc(Point o,Poly &p){ return AreaIntersecPolyC(p,Circle(o,r)); } int main(){ int n; Point o; scanf("%d%lf",&n,&r); Poly p(n); for(int i=0;i<n;++i){ p[i].read(); o=o+p[i]; } o=o/n; double ans=calc(o,p); default_random_engine e(time(0)); uniform_real_distribution<double> u(0,1); bool flag; double T=30; int cnt=2048; while(cnt--){ Point o2=o; double res=ans; flag=0; double ang=2*PI*u(e); double step=2*PI/50; for(int i=0;i<50;++i){ Point o3=o2+Point(T*cos(ang),T*sin(ang)); double tmp=calc(o3,p); if(tmp>res){ res=tmp; flag=1; o2=o3; } ang+=step; } if(!flag)T*=0.9; else{ ans=res; o=o2; T*=1.2; } } printf("%.6f",ans); return 0; } /* 4 4 0 0 6 0 6 6 0 6 35.759506 3 1 0 0 2 1 1 3 2.113100 3 1 0 0 100 1 99 1 0.019798 4 1 0 0 100 10 100 12 0 1 3.137569 10 10 0 0 10 0 20 1 30 3 40 6 50 10 60 15 70 21 80 28 90 36 177.728187 */