[Gym-101158J]Coverthe Polygon with Your Disk——梯度下降,模拟退火

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

*/
模拟退火

 

posted @ 2020-04-17 16:27  _ZPENG  阅读(263)  评论(0编辑  收藏  举报