poj 3384 Feng Shui(半平面交)

#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
using namespace std;
const double epsi=1e-10;
const int maxn=1010;
const double pi=acos(-1.0);
inline int sign(const double &x){
    if(x>epsi) return 1;
    if(x<-epsi) return -1;
    return 0;
}
struct point{
    double x,y;
    point(double xx=0,double yy=0):x(xx),y(yy){}
    point operator +(const point &op2) const{
        return point(x+op2.x,y+op2.y);
    }
    point operator -(const point &op2) const{
        return point(x-op2.x,y-op2.y);
    }
    double operator *(const point &op2) const{
        return x*op2.x+y*op2.y;
    }
    point operator *(const double &d) const{
        return point(x*d,y*d);
    }
    point operator /(const double &d) const{
        return point(x/d,y/d);
    }
    double operator ^(const point &op2) const{
        return x*op2.y-y*op2.x;
    }
    bool operator ==(const point &op2) const{
        return sign(x-op2.x)==0&&sign(y-op2.y)==0;
    }
};
int n;
double r;
point p[maxn];
int t[2];
point plane[2][maxn],q1,q2;
inline double sqr(const double &x){
    return x*x;
}
inline double mul(const point &p0,const point &p1,const point &p2){
    return (p1-p0)^(p2-p0);
}
inline double dot(const point &p0,const point &p1,const point &p2){
    return (p1-p0)*(p2-p0);        //p1p0与p2p0的点积
}
inline double dis2(const point &p0,const point &p1){
    return sqr(p0.x-p1.x)+sqr(p0.y-p1.y);
}
inline double dis(const point &p0,const point &p1){
    return sqrt(dis2(p0,p1));
}
inline double dis(const point &p0,const point &p1,const point &p2){
    if(sign(dot(p1,p0,p2)<0)) return dis(p0,p1); //若p1p0与p1p2的夹角超过90度，返回p1p0
    if(sign(dot(p2,p0,p1))<0) return dis(p0,p2); //若p2p0与p2p1的夹角超过90度，返回p2p0
    return fabs(mul(p0,p1,p2)/dis(p1,p2));  //返回p0至p1p2的垂线长度
}
inline point rotate(const point &p,const double &ang){  //点p旋转ang角度后的点
    return point(p.x*cos(ang)-p.y*sin(ang),p.x*sin(ang)+p.y*cos(ang));
}
inline void translation(const point &p1,const point &p2,const double &d,point &q1,point &q2){
    q1=p1+rotate(p2-p1,pi/2)*d/dis(p1,p2);   //p2p1向内推进d后形成直线q2q1
    q2=q1+p2-p1;
}
inline void cross(const point &p1,const point &p2,const point &p3,const point &p4,point &p){
    double a1=mul(p1,p3,p4),a2=mul(p2,p3,p4);
    p.x=(a1*p2.x-a2*p1.x)/(a1-a2);
    p.y=(a1*p2.y-a2*p1.y)/(a1-a2);    //利用公式计算交点
}
inline int half_plane_cross(point *a,int n,point *b,const point &p1,const point &p2){
    int newn=0;
    for(int i=0,j;i<n;++i){
        if(sign(mul(a[i],p1,p2))>=0){
            b[newn++]=a[i];
            continue;
        }
        j=i-1;
        if(j==-1) j=n-1;   //计算i的左相邻点j
        if(sign(mul(a[j],p1,p2))>0){
            cross(p1,p2,a[j],a[i],b[newn++]);
        }
        j=i+1;
        if(j==n) j=0;      //计算i右相邻点j
        if(sign(mul(a[j],p1,p2))>0)
            cross(p1,p2,a[j],a[i],b[newn++]);
    }
    return newn;
}
int main()
{
   scanf("%d%lf",&n,&r);
   for(int i=0;i<n;i++)
    scanf("%lf%lf",&p[i].x,&p[i].y);
   p[n]=p[0];
   int o1=0,o2;
   t[0]=4;
   plane[0][0]=point(-1e3,-1e3);
   plane[0][1]=point(1e3,-1e3);
   plane[0][2]=point(1e3,1e3);
   plane[0][3]=point(-1e3,1e3);
   for(int i=0;i<n;i++){
       o2=o1^1;
       translation(p[i+1],p[i],r,q1,q2);
       t[o2]=half_plane_cross(plane[o1],t[o1],plane[o2],q1,q2);
       o1=o2;
   }
   double maxd=-1,curd;
   for(int i=0;i<t[o1];i++)
   for(int j=i;j<t[o1];j++){
      curd=dis2(plane[o1][i],plane[o1][j]);
      if(sign(curd-maxd)>0){
        maxd=curd;
        q1=plane[o1][i];q2=plane[o1][j];
      }
   }
   printf("%.4f %.4f %.4f %.4f\n",q1.x,q1.y,q2.x,q2.y);
   return 0;
}