AOJ 0253 : Ant Nest

問題リンク : Ant Nest | Aizu Online Judge

問題概要 :
略

解法 :
重心求めてくるくるしながら二分探索して計算する

まず最初に、連続した３点以上が平行であるようなものが入力として与えられるので、それらを圧縮して１つの線分にする

各辺について、そこをケースの底辺として高さを計算する
ケースを回転するために、凸の重心を求め重心が原点にくるように移動する
現在底辺となっている辺の右隣の辺に垂直なベクトルを求め、原点からその方向に適当にのばし、その辺との交点を求める
その交点が原点の真下にくるように回転する

高さは二分探索で求める

角度計算で Runtime Error を出したがそれ以外にバグはなかったようで良かった

Sample Input 01:

8 1 1
0 0
1 0
2 0
2 1
2 2
1 2
0 2
0 1
0 0 0

Sample Output 01:

0.5000000000

Sample Input 02:

5 2 2
0 0
1 0
2 2
1 3
0 2
0 0 0

Sample Output 02:

0.8284271266

Sample Input 03:

4 1 1
0 0
1 0
1 1
0 1
0 0 0

Sample Output 03:

1.0000000000

Sample Input 04:

3 1 1
-1 -1
1 -1
0 1
0 0 0

Sample Output 04:

0.5857864320

コード :

#include<bits/stdc++.h>

#define REP(i,s,n) for(int i=s;i<n;i++)
#define rep(i,n) REP(i,0,n)
#define EPS (1e-8)
#define equals(a,b) (fabs((a)-(b)) < EPS)
#define COUNTER_CLOCKWISE 1
#define CLOCKWISE -1 
#define ONLINE_BACK 2
#define ONLINE_FRONT -2
#define ON_SEGMENT 0

using namespace std;

class Point{
public:
  double x,y;

  Point(double x = 0,double y = 0): x(x),y(y){}

  Point operator + (Point p){return Point(x+p.x,y+p.y);}
  Point operator - (Point p){return Point(x-p.x,y-p.y);}
  Point operator * (double a){return Point(a*x,a*y);}
  Point operator / (double a){return Point(x/a,y/a);}
  Point operator * (Point a){ return Point(x * a.x - y * a.y, x * a.y + y * a.x); }

  bool operator < (const Point& p) const{ return !equals(x,p.x)?x<p.x:(!equals(y,p.y)&&y<p.y); }

  bool operator == (const Point& p)const{ return fabs(x-p.x) < EPS && fabs(y-p.y) < EPS; }

};

struct Segment{
  Point p1,p2;
  Segment(Point p1 = Point(),Point p2 = Point()):p1(p1),p2(p2){}
  bool operator <  (const Segment& s) const { return ( p2 == s.p2 ) ? p1 < s.p1 : p2 < s.p2; }
  bool operator == (const Segment& s) const { return ( s.p1 == p1 && s.p2 == p2 ) || ( s.p1 == p2 && s.p2 == p1 ); }

};

typedef Point Vector;
typedef Segment Line;
typedef vector<Point> Polygon;

double dot(Point a,Point b){ return a.x*b.x + a.y*b.y; }

double cross(Point a,Point b){ return a.x*b.y - a.y*b.x; }

double norm(Point a){ return a.x*a.x+a.y*a.y; }

double abs(Point a){ return sqrt(norm(a)); }

Point rotate(Point a,double rad){ return Point(cos(rad)*a.x - sin(rad)*a.y,sin(rad)*a.x + cos(rad)*a.y); }

double toRad(double agl){ return agl*M_PI/180.0; }

double getArg(Point a,Point b,Point c){
  double arg1 = atan2(b.y-a.y,b.x-a.x);
  double arg2 = atan2(c.y-b.y,c.x-b.x);
  double arg = fabs( arg1 - arg2 );
  while( arg > M_PI ) arg -= 2.0 * M_PI;
  return fabs(arg);
}

int ccw(Point p0,Point p1,Point p2){
  Point a = p1-p0;
  Point b = p2-p0;
  if(cross(a,b) > EPS)return COUNTER_CLOCKWISE;
  if(cross(a,b) < -EPS)return CLOCKWISE;
  if(dot(a,b) < -EPS)return ONLINE_BACK;
  if(norm(a) < norm(b))return ONLINE_FRONT;
  return ON_SEGMENT;
}

bool intersectLL(Line l, Line m) {
  return abs(cross(l.p2-l.p1, m.p2-m.p1)) > EPS || // non-parallel
         abs(cross(l.p2-l.p1, m.p1-l.p1)) < EPS;   // same line
}

bool intersectLS(Line l, Line s) {
  return cross(l.p2-l.p1, s.p1-l.p1)*       // s[0] is left of l
         cross(l.p2-l.p1, s.p2-l.p1) < EPS; // s[1] is right of l
}
bool intersectLP(Line l,Point p) {
  return abs(cross(l.p2-p, l.p1-p)) < EPS;
}
bool intersectSS(Line s, Line t) {
  return ccw(s.p1,s.p2,t.p1)*ccw(s.p1,s.p2,t.p2) <= 0 &&
         ccw(t.p1,t.p2,s.p1)*ccw(t.p1,t.p2,s.p2) <= 0;
}
bool intersectSP(Line s, Point p) {
  return abs(s.p1-p)+abs(s.p2-p)-abs(s.p2-s.p1) < EPS; // triangle inequality
}

Point projection(Line l,Point p) {
  double t = dot(p-l.p1, l.p1-l.p2) / norm(l.p1-l.p2);
  return l.p1 + (l.p1-l.p2)*t;
}
Point reflection(Line l,Point p) {
  return p + (projection(l, p) - p) * 2;
}
double distanceLP(Line l, Point p) {
  return abs(p - projection(l, p));
}
double distanceLL(Line l, Line m) {
  return intersectLL(l, m) ? 0 : distanceLP(l, m.p1);
}

double distanceLS(Line l, Line s) {
  if (intersectLS(l, s)) return 0;
  return min(distanceLP(l, s.p1), distanceLP(l, s.p2));
}
double distanceSP(Line s, Point p) {
  Point r = projection(s, p);
  if (intersectSP(s, r)) return abs(r - p);
  return min(abs(s.p1 - p), abs(s.p2 - p));
}

double distanceSS(Line s, Line t) {
  if (intersectSS(s, t)) return 0;
  return min(min(distanceSP(s, t.p1), distanceSP(s, t.p2)),
             min(distanceSP(t, s.p1), distanceSP(t, s.p2)));
}

Point crosspoint(Line l,Line m){
  double A = cross(l.p2-l.p1,m.p2-m.p1);
  double B = cross(l.p2-l.p1,l.p2-m.p1);
  if(abs(A) < EPS && abs(B) < EPS){
    vector<Point> vec;
    vec.push_back(l.p1),vec.push_back(l.p2),vec.push_back(m.p1),vec.push_back(m.p2);
    sort(vec.begin(),vec.end());
    assert(vec[1] == vec[2]); //同じセグメントかもよ
    return vec[1];
    //return m.p1;
  }
  if(abs(A) < EPS)assert(false);
  return m.p1 + (m.p2-m.p1)*(B/A);
}

double getArea(vector<Point>& vec) {
  double sum = 0;
  for(int i=0;i<vec.size();i++)
    sum += cross(vec[i],vec[(i+1)%vec.size()]);
  return fabs(sum)/2.0;
}

double cross3p(Point p,Point q,Point r) { return (r.x-q.x) * (p.y -q.y) - (r.y - q.y) * (p.x - q.x); }

bool collinear(Point p,Point q,Point r) { return fabs(cross3p(p,q,r)) < EPS; }

bool ccwtest(Point p,Point q,Point r){ return cross3p(p,q,r) > 0; }
 
bool onSegment(Point p,Point q,Point r){
  return collinear(p,q,r) && equals(sqrt(pow(p.x-r.x,2)+pow(p.y-r.y,2)) + sqrt(pow(r.x-q.x,2) + pow(r.y-q.y,2) ),sqrt(pow(p.x-q.x,2)+pow(p.y-q.y,2)) ) ;
}


Point getCentroidOfConvex(Polygon& poly){
  double area = getArea(poly);
  int V = poly.size();
  assert( !equals(area,0.0) );
  double x = 0, y = 0;
  rep(i,(int)poly.size()) {
    x += ( poly[i].x + poly[(i+1)%V].x ) * ( poly[i].x*poly[(i+1)%V].y - poly[(i+1)%V].x*poly[i].y );
    y += ( poly[i].y + poly[(i+1)%V].y ) * ( poly[i].x*poly[(i+1)%V].y - poly[(i+1)%V].x*poly[i].y );
  }
  return Point(x/(6.0*area),y/(6.0*area));
}

bool LT(double a,double b) { return !equals(a,b) && a < b; }
bool LTE(double a,double b) { return equals(a,b) || a < b; }

int V;
double d,vol;
Polygon poly;
Point ps[110];

double lets_compute_arg(int L,int R){
  assert( L != R );
  Point lp = poly[L], rp = poly[R];
  rp = rp - lp;
  lp = Point(0,0);
  Vector e = rp / abs(rp);
  e = e * 100000;
  e = rotate(e,toRad(90.0));
  Line line = Line(e,e * -1);
  Point cp = crosspoint(line,Line(poly[L],poly[R]));
  double ret = getArg(Point(0,-1000),Point(0,0),cp);
  return ret + toRad(180.0);
}


double getNewVolume(double h){
  Segment cut_line = Segment(Point(-20000,h),Point(20000,h));
  int nV = 0;
  
  int index = 0;
  int cross_ps[2];

  rep(i,V) {
    Segment seg = Segment(poly[i],poly[(i+1)%V]);
    if( intersectSS(cut_line,seg) ) { // 平行な場合はない
      Point cp = crosspoint(cut_line,seg);
      if( !( cp == seg.p1 || cp == seg.p2 ) ) {
        ps[nV++] = cp;
      } 
    } 
    ps[nV++] = poly[(i+1)%V];
  }

  rep(i,nV) if( equals(cut_line.p2.y,ps[i].y) ) { cross_ps[index++] = i; }
  
  int sp = -1;
  double maxi_h = -1e30;
  rep(i,nV) maxi_h = max(maxi_h,ps[i].y);
  rep(i,nV) if( equals(maxi_h,ps[i].y) ) { sp = i; break; }

  bool found = false;
  int cur = sp;
  vector<Point> npoly;
  rep(i,nV) {
    bool same_point = ( cross_ps[0] == cur || cross_ps[1] == cur );
    if( same_point && found ) { npoly.push_back(ps[cur]); break; }
    if( same_point ) found = true;
    if( found ) npoly.push_back(ps[cur]);
    cur = ( cur + 1 ) % nV;
  }

  double area = getArea(npoly);
  return area * d;
}

void compute(){

  bool full = false;
  {
    double area = getArea(poly) * d;
    if( equals(area,vol) ) full = true;
  }

  double maxi = 0;
  Point centroid = getCentroidOfConvex(poly);
  rep(i,V) poly[i] = poly[i] - centroid;
  centroid = Point(0,0);

  rep(_,V){

    double maxi_y = -1e30, mini_y = 1e30;
    rep(i,V) maxi_y = max(maxi_y,poly[i].y), mini_y = min(mini_y,poly[i].y);
    double low_y = mini_y, high_y = maxi_y;

    if( !full ) {
      rep(__,60){
        double mid_y = ( low_y + high_y ) / 2.0;

        double new_volume = getNewVolume(mid_y);

        if( equals(new_volume,vol) ) { maxi = max(maxi,fabs(mid_y-mini_y)); break; }
        if( LT(new_volume,vol) ) low_y  = mid_y;
        else                     high_y = mid_y;
      }
    } else {
      rep(j,V) maxi = max(maxi,abs(mini_y-poly[j].y));
    }

    
    int left_p = -1, right_p = -1;
    rep(i,V) if( equals(mini_y,poly[i].y) ) {
      if( left_p == -1 ) left_p = i;
      if( right_p == -1 ) right_p = i;
      if( LT(poly[i].x,poly[left_p].x)  ) left_p  = i;
      if( LT(poly[right_p].x,poly[i].x) ) right_p = i;
    }

    ( left_p  += 1 ) %= V;
    ( right_p += 1 ) %= V;
    double rad = lets_compute_arg(left_p,right_p);
    rep(i,V) poly[i] = rotate(poly[i],rad);

  }
  printf("%.10f\n",maxi);
}

void fix(){
  bool update = true;
  while( update ){
    update = false;
    rep(i,poly.size()){
      int j = ( i + 1 ) % poly.size();
      int k = ( i + 2 ) % poly.size();
      if( equals(cross(poly[j]-poly[i],poly[k]-poly[i]),0.0) ) {
        update = true;
        poly.erase(poly.begin()+j);
        goto Skip;
      }
    }
  Skip:;
  }
  V = poly.size();
}

int main(){
  int a,b,c;
  while( scanf("%d %d %d",&a,&b,&c), a|b|c ){
    V = a, d = b, vol = c;
    poly.clear();
    poly.resize(V);
    rep(i,V) scanf("%lf %lf",&poly[i].x,&poly[i].y);
    fix();
    compute();
  }
  return 0;
}