LA 3890 半平面交

二分查询答案，判断每一个新形成的向量合在一块能否形成半平面交

  1 #include <iostream>
  2 #include <cstdio>
  3 #include <cstring>
  4 #include <cmath>
  5 #include <algorithm>
  6 using namespace std;
  7 #define N 110
  8 #define eps 1e-7
  9
 10 int dcmp(double x)
 11 {
 12     if(fabs(x)<eps) return 0;
 13     else return x<0?-1:1;
 14 }
 15
 16 struct Point{
 17     double x , y;
 18     Point(double x=0 , double y=0):x(x),y(y){}
 19 }po[N] , poly[N];
 20
 21 typedef Point Vector;
 22 Vector vec[N]; //记录每条边上对应的法向量
 23
 24 Vector operator+(Vector a , Vector b) { return Vector(a.x+b.x , a.y+b.y); }
 25 Vector operator-(Point a , Point b) { return Vector(a.x-b.x , a.y-b.y); }
 26 Vector operator*(Vector a , double b) { return Vector(a.x*b , a.y*b); }
 27 Vector operator/(Vector a , double b) { return Vector(a.x/b , a.y/b); }
 28 bool operator==(const Point &a , const Point &b) { return dcmp(a.x-b.x) == 0 && dcmp(a.y-b.y) == 0; }
 29
 30 double Dot(Vector a  , Vector b) { return a.x*b.x+a.y*b.y; }
 31 double Cross(Vector a , Vector b) { return a.x*b.y - b.x*a.y; }
 32 double Length(Vector a) { return sqrt(Dot(a,a)); }
 33 double Angle(Vector A , Vector B) { return acos(Dot(A,B)) / Length(A) / Length(B); }
 34 double Area2(Point A , Point B , Point C) { return Cross(B-A , C-A); }
 35
 36 Vector Rotate(Vector A , double rad) { return Vector(A.x*cos(rad)-A.y*sin(rad) , A.x*sin(rad)+A.y*cos(rad)); }
 37
 38 Vector Normal(Vector a)
 39 {
 40     double l = Length(a);
 41     return Vector(-a.y/l , a.x/l);
 42 }
 43
 44 struct Line{
 45     Point P;
 46     Vector v;
 47     double ang;
 48     Line(){}
 49     Line(Point P , Vector v):P(P),v(v){ang = atan2(v.y,v.x);}
 50     bool operator<(const Line &m)const {
 51         return ang<m.ang;
 52     }
 53 }line[N] , L[N];
 54
 55 bool OnLeft(Line L , Point P)
 56 {
 57     return Cross(L.v , P-L.P) > 0;
 58 }
 59
 60 Point GetIntersection(Line a , Line b)
 61 {
 62     Vector u = a.P-b.P;
 63     double t = Cross(b.v , u) / Cross(a.v , b.v);
 64     return a.P+a.v*t;
 65 }
 66
 67 /***半平面交的主过程，返回形成半平面交点的个数，无法形成就返回0***/
 68 int HalfplaneIntersection(Line *L , int n , Point *poly)
 69 {
 70     sort(L , L+n);
 71     int first , last;
 72     Point *p = new Point[n];
 73     Line *q = new Line[n];
 74     q[first=last=0] = L[0];
 75     for(int i=1 ; i<n ; i++)
 76     {
 77         while(first<last && !OnLeft(L[i] , p[last-1])) last--;
 78         while(first<last && !OnLeft(L[i] , p[first])) first++;
 79         q[++last] = L[i];
 80         if(fabs(Cross(q[last].v , q[last-1].v)) < eps)
 81         {
 82             last--;
 83             if(OnLeft(q[last] , L[i].P)) q[last]=L[i];
 84         }
 85         if(first<last) p[last-1] = GetIntersection(q[last-1] , q[last]);
 86     }
 87     while(first < last && !OnLeft(q[first] , p[last-1])) last--;
 88     //删除无用平面
 89     if(last-first<=1) return 0;
 90     p[last] = GetIntersection(q[last] , q[first]);
 91
 92     //从deque复制到输出中
 93     int m=0;
 94     for(int i=first ; i<=last ; i++) poly[m++] = p[i];
 95     return m;
 96 }
 97
 98 int main()
 99 {
100    // freopen("a.in" , "r" , stdin);
101     int n;
102     while(scanf("%d" , &n) , n)
103     {
104         for(int i=0 ; i<n ; i++)
105             scanf("%lf%lf" , &po[i].x , &po[i].y);
106
107         for(int i=0 ; i<n ; i++) vec[i] = Normal(po[(i+1)%n]-po[i]);
108         double l=0 , r=20000;
109         while(r-l>eps){
110             double m=(l+r)/2;
111             for(int i=0 ; i<n ; i++) L[i] = Line(po[i]+vec[i]*m , po[(i+1)%n]-po[i]);
112             if(HalfplaneIntersection(L , n , poly)>0) l=m;
113             else r=m;
114         }
115         printf("%.6f\n" , l);
116     }
117     return 0;
118 }