Gym-101158J&&CSU2300 Cover the Polygon with Your Disk 三分套三分

题目连接：http://acm.csu.edu.cn:20080/csuoj/problemset/problem?pid=2300

题目大意：给定一个凸多边形和一个圆，移动圆的位置使得和多边形的交最大，求交的最大面积。

此题显然是个凸函数，所以可以三分套三分，也可模拟退火，感觉还是三分套三分稳一点。

关于如何求多边形和圆的交，这里有一篇博客：http://www.cnblogs.com/lxglbk/archive/2012/08/12/2634192.html
代码里也写了点注释，可以参考一下。

AC代码：

#include<bits/stdc++.h>
using namespace std;
const int N=15;
const double pi=acos(-1.0);
const double eps=1e-9;
const double inf=1e9;
struct Point{
    double x,y;
    Point(){};
    Point(double x,double y):x(x),y(y){}
    Point operator-(const Point p)const{
        return Point(x-p.x,y-p.y);
    }
    Point operator+(const Point p)const{
        return Point(x+p.x,y+p.y);
    }
    double operator*(const Point p)const{
        return x*p.x+y*p.y;
    }
    double operator^(const Point p)const{
        return x*p.y-y*p.x;
    }
}a[N],p[N];
inline double dis(Point a,Point b){
    return sqrt((b-a)*(b-a));
}
inline double cross(Point a,Point b,Point c){//叉积,关于a
    return (b-a)^(c-a);
}
inline double dot(Point a,Point b,Point c){//点积，关于a
    return (b-a)*(c-a);
}
inline bool online(Point a,Point b,Point c){//c是否在线段ab上
    return cross(a,b,c)==0&&c.x>=min(a.x,b.x)&&\
    c.y>=min(a.y,b.y)&&c.x<=max(a.x,b.x)&&c.y<=max(a.y,b.y);
}
double area(Point b,Point c,double r){//圆和一个顶点为圆心的三角形的交
    Point a(0.0,0.0);
    if(dis(b,c)<eps)
        return 0.0;
    double h=fabs(cross(a,b,c))/dis(b,c);
    if(dis(a,b)>r-eps&&dis(a,c)>r-eps){//两个端点都在圆的外面
        double angle=acos(dot(a,b,c)/dis(a,b)/dis(a,c));
        if(h>r-eps){
            return 0.5*r*r*angle;
        }
        else if(dot(b,a,c)>0&&dot(c,a,b)){
            double angle1=2*acos(h/r);
            return 0.5*r*r*fabs(angle-angle1)+0.5*r*r*sin(angle1);
        }
        else{
            return 0.5*r*r*angle;
        }
    }
    else if(dis(a,b)<r+eps&&dis(a,c)<r+eps){   //两个端点都在圆内的情况
        return 0.5*fabs(cross(a,b,c));
    }
    else{//一个端点在圆外，一个端点在圆内的情况
        if(dis(a,b)>dis(a,b)){//默认b在圆内
            swap(b,c);
        }
        if(fabs(dis(a,b)<eps)){//ab距离为零直接返回
            return 0.0;
        }
        if(dot(b,a,c)<eps){//角abc为钝角
            double angle1=acos(h/dis(a,b));
            double angle2=acos(h/r)-angle1;
            double angle3=acos(h/dis(a,c))-angle2;
            return 0.5*r*dis(a,b)*sin(angle1+angle2)+0.5*r*r*angle3;
        }
        else{           //角abc为锐角
            double angle1=acos(h/dis(a,b));
            double angle2=acos(h/r);
            double angle3=acos(h/dis(a,c))-angle2;
            return 0.5*r*dis(a,b)*sin(angle1+angle2)+0.5*r*r*angle3;
        }
    }
}
int n;
double r;
double gets(Point o){//圆和多边形的交
    for(int i=0;i<n;i++)
        p[i]=a[i]-o;
    o=Point(0.0,0.0);
    double sum=0;
    for(int i=0;i<n;i++){
        int nxt=(i+1)%n;
        double s=area(p[i],p[nxt],r);
        if(cross(o,p[i],p[nxt])>0)
            sum+=s;
        else
            sum-=s;
    }
    return fabs(sum);
}
double cal(double x){
    double ly=inf,ry=-inf;
    for(int i=0,j=1;i<n;i++,j=(j+1)%n){
        if((a[i].x-x)*(a[j].x-x)<=0){//仅当边和X=x相交时才要更新ly,ry
            double tmp=a[i].y+(a[j].y-a[i].y)/(a[j].x-a[i].x)*(x-a[i].x);
            ly=min(ly,tmp);
            ry=max(ry,tmp);
        }
    }
    double s1,s2;
    while(ry-ly>eps){
        double m1=(ly+ly+ry)/3;
        double m2=(ly+ry+ry)/3;
        s1=gets(Point(x,m1)),s2=gets(Point(x,m2));
        if(s1>s2)
            ry=m2;
        else
            ly=m1;
    }
    return s1;
}
int main(){
    double lx,rx;
    scanf("%d%lf",&n,&r);
    lx=inf,rx=-inf;
    for(int i=0;i<n;i++){
        scanf("%lf%lf",&a[i].x,&a[i].y);
        lx=min(lx,a[i].x);
        rx=max(rx,a[i].x);
    }
    double s1,s2;
    while(rx-lx>eps){
        double m1=(lx+lx+rx)/3;
        double m2=(lx+rx+rx)/3;
        s1=cal(m1),s2=cal(m2);
        if(s1>s2)
            rx=m2;
        else
            lx=m1;
    }
    printf("%.6f\n",s1);
    return 0;
}