POJ 3384 Feng Shui

http://poj.org/problem?id=3384

题意:给一个凸包，求往里面放两个圆(可重叠)的最大面积时的两个圆心坐标。

思路:先把凸包边往内推R，做半平面交，然后做旋转卡壳，此时得到最大距离的点对，就是圆心坐标。

PS:最大长度的初始值要设置为负数，因为距离有可能退化到0，就像这组数据

4 1

0 0

2 0

2 2

0 2

#include<cstdio>
#include<iostream>
#include<cmath>
#include<cstring>
#include<algorithm>
const double Pi=acos(-1);
double R;
int n,tot;
struct Point{
    double x,y;
    Point(){}
    Point(double x0,double y0):x(x0),y(y0){}
}p[200005];
struct Line{
    Point s,e;
    double slop;
    Line(){}
    Line(Point s0,Point e0):s(s0),e(e0){}
}L[200005],l[200005],c[200005];
int read(){
    int t=0,f=1;char ch=getchar();
    while (ch<'0'||ch>'9'){if (ch=='-')f=-1;ch=getchar();}
    while ('0'<=ch&&ch<='9'){t=t*10+ch-'0';ch=getchar();}
    return t*f;
}
Point operator /(Point p1,double x){
    return Point(p1.x/x,p1.y/x);
}
Point operator *(Point p,double x){
    return Point(p.x*x,p.y*x);
}
double operator *(Point p1,Point p2){
    return p1.x*p2.y-p1.y*p2.x;
}
Point operator -(Point p1,Point p2){
    return Point(p1.x-p2.x,p1.y-p2.y);
}
Point operator +(Point p1,Point p2){
    return Point(p1.x+p2.x,p1.y+p2.y);
}
bool cmp(Line p1,Line p2){
    if (p1.slop!=p2.slop) return p1.slop<p2.slop;
    else return (p1.e-p1.s)*(p2.e-p1.s)<=0;
}
Point inter(Line p1,Line p2){
    double k1=(p2.e-p1.s)*(p1.e-p1.s);
    double k2=(p1.e-p1.s)*(p2.s-p1.s);
    double t=(k2)/(k1+k2);
    double x=p2.s.x+(p2.e.x-p2.s.x)*t;
    double y=p2.s.y+(p2.e.y-p2.s.y)*t;
    return Point(x,y);
}
bool jud(Line p1,Line p2,Line p3){
    Point p=inter(p1,p2);
    return (p-p3.s)*(p3.e-p3.s)>0;
}
void phi(){
    std::sort(l+1,l+1+tot,cmp);
    int cnt=1;
    for (int i=2;i<=tot;i++)
     if (l[i].slop!=l[i-1].slop)
      l[++cnt]=l[i];
    int L=1,R=2;c[L]=l[1];c[R]=l[2];
    for (int i=3;i<=cnt;i++){
        while (L<R&&jud(c[R],c[R-1],l[i])) R--;
        while (L<R&&jud(c[L],c[L+1],l[i])) L++;
        c[++R]=l[i];
    }  
    while (L<R&&jud(c[R],c[R-1],c[L])) R--;
    while (L<R&&jud(c[L],c[L+1],c[R])) L++;
    tot=0;
    c[R+1]=c[L];
    for (int i=L;i<=R;i++)
     p[++tot]=inter(c[i],c[i+1]);
}
double sqr(double x){
    return x*x;
}
double dis(Point p){
    return sqrt(sqr(p.x)+sqr(p.y));
}
Point turn(Point p,double ang){
    double Cos=cos(ang),Sin=sin(ang);
    double x=Cos*p.x-Sin*p.y;
    double y=Cos*p.y+Sin*p.x;
    return Point(x,y);
}
Point e(Point p){
    double len=dis(p);p=p/len;return p;
}
double dis(Point p1,Point p2){
    return dis(p1-p2);
}
void rc(){
    p[tot+1]=p[1];
    int k=2;
    double mx=-1;
    Point ans1,ans2;
    for (int i=1;i<=tot;i++){
        while (fabs((p[i%tot+1]-p[i])*(p[k]-p[i]))<fabs((p[i%tot+1]-p[i])*(p[k%tot+1]-p[i]))) k=(k)%tot+1;
        if (mx<dis(p[i],p[k])){
            mx=dis(p[i],p[k]);
            ans1=p[i];
            ans2=p[k];
        }
    }
    printf("%.4f %.4f %.4f %.4f",ans1.x,ans1.y,ans2.x,ans2.y);
}
int main(){
    n=read(),R=read();
    for (int i=1;i<=n;i++) p[i].x=read(),p[i].y=read();
    for (int i=1;i<=n/2;i++) std::swap(p[i],p[n-i+1]);
    p[n+1]=p[1];
    for (int i=1;i<=n;i++)
     l[++tot]=Line(p[i],p[i+1]),l[tot].slop=atan2(l[tot].e.y-l[tot].s.y,l[tot].e.x-l[tot].s.x); 
    for (int i=1;i<=n;i++){
     Point p=e(turn((l[i].e-l[i].s),Pi/2.0))*R;
     l[i].s=l[i].s+p;
     l[i].e=l[i].e+p;
    }
    phi();
    rc();
}