土下座しながら探索中

主に競技プログラミング

AOJ 2246 : Alice and Bomb

問題リンク : Alice and Bomb | Aizu Online Judge

問題概要:
建物、人、爆弾があり、建物は2次元平面上の多角形として、人と爆弾は2次元平面上の点として与えられる
爆弾は爆発すると全方位に爆風が広がる
建物に爆風が当たるとそこでそこの爆風は遮られる
人が爆風に当たらない場所まで移動する際にかかる移動距離の最小値を求めよ
ただし、爆弾が建物を表す多角形の辺の延長線上にあることはない

解法 :
幾何 + dijkstra
爆風が当たらない場所の候補は公式の解説スライドにある通り
あとは候補点を含めてdijkstraして答えを得る
次の2つの関数を作っておくと良い
1. 点を引数とし、その点が爆風を遮っている ( その点がゴールになっている ) か判定する関数
2. 2つの点を引数とし、その2点間を直線で移動することができるか判定する関数

2.の関数では点の内包判定をすることになるが、角度で点の内包判定をするととても重い (実装にもよるだろうが)
というのもacosが非常に重い、これだけで10secは違った
自分の持っているライブラリは角度のものだけだったのでspagettie sourceのものを利用した

dijkstraをする際には、今いる点がゴールならばすぐにそこで処理を終えること
でないとTLEする
自分の場合は、dijkstraは各頂点間だけでその後に他の候補店を試したのだが
dijkstraで人から全ての点の最短距離をもとめていたらTLEした
(これが原因で数ヶ月解けなかった)


コード:

#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;

ostream& operator << (ostream& os,const Point& a){ os << "(" << a.x << "," << a.y << ")"; }

ostream& operator << (ostream& os,const Segment& a){ os << "( " << a.p1 << " , " << a.p2 << " )"; }

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)); }

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);
}


//cross product of pq and pr
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); }
  
//returns true if point r is on the same line as the line pq
bool collinear(Point p,Point q,Point r) { return fabs(cross3p(p,q,r)) < EPS; }
  
//returns true if point t is on the left side of line pq
bool ccwtest(Point p,Point q,Point r){
  return cross3p(p,q,r) > 0; //can be modified to accept collinear points
}
 
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)) ) ;
}
  

double angle(Point a,Point b,Point c) {
  double ux = b.x - a.x, uy = b.y - a.y;
  double vx = c.x - a.x, vy = c.y - a.y;
  return acos((ux*vx + uy*vy)/sqrt((ux*ux + uy*uy) * (vx*vx + vy*vy)));
}  
  
//多角形poly内(線分上も含む)に点pが存在するかどうは判定する  
bool inPolygon(Polygon poly,Point p,bool flag=true){
  if((int)poly.size() == 0)return false;
  if( flag ) rep(i,poly.size()) if(onSegment(poly[i],poly[(i+1)%poly.size()],p))return true;
  double sum = 0;
  for(int i=0; i < (int)poly.size() ;i++) {
    if( equals(cross(poly[i]-p,poly[(i+1)%poly.size()]-p),0.0) ) continue; // 3点が平行だとangleがnanを返しsumがnanになり死ぬ
    if( cross3p(p,poly[i],poly[(i+1)%poly.size()]) < 0 ) sum -= angle(p,poly[i],poly[(i+1)%poly.size()]);
    else                                                 sum += angle(p,poly[i],poly[(i+1)%poly.size()]);
  }
  // あまり誤差を厳しくしすぎると良くないので以下のほうが良い 
  const double eps = 1e-5;
  return (fabs(sum - 2*M_PI) < eps || fabs(sum + 2*M_PI) < eps);
}  






enum { OUT, ON, IN };
int contains(Polygon& P, Point& p) {
  bool in = false;
  for (int i = 0; i < P.size(); ++i) {
    Point a = P[i] - p, b = P[(i+1)%P.size()] - p;
    if (a.y > b.y) swap(a, b);
    if (a.y <= 0 && 0 < b.y)
      if (cross(a, b) < 0) in = !in;
    if (cross(a, b) == 0 && dot(a, b) <= 0) return ON;
  }
  return in ? IN : OUT;
}
// -------------------

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

bool DEBUG = false;

struct Data {
  Point cur;
  double weight;
  int a,b;
  bool operator < ( const Data &data ) const { return LT(data.weight,weight); }
};

const double DINF = 1e30;
int N;
Point bomb;
vector<Polygon> polies;
vector<vector<double> > mindist;
vector<vector<bool> > goals;

bool isValidMove(Point src,Point dst){

  Segment ururu_beam = Segment(src,dst);
  Point mp = ( src + dst ) * 0.5; // $
  rep(i,N) {
    int M = polies[i].size();
    bool in = false;
    bool on = false;
    rep(j,M) {
      Segment seg = Segment(polies[i][j],polies[i][(j+1)%M]);

      Point a = polies[i][j] - mp, b = polies[i][(j+1)%M] - mp;
      if( a.y > b.y ) swap(a,b);
      if( a.y <= 0 && 0 < b.y )  if( cross(a,b) < 0 ) in = !in;
      if( cross(a,b) == 0 && dot(a,b) <= 0 ) on = true;

      // うるるビームと平行なのであればどの様な状況であっても正しい移動
      if( equals(cross(ururu_beam.p1-ururu_beam.p2,seg.p1-seg.p2),0.0) ) continue;
      seg = Segment(polies[i][j],polies[i][(j-1+M)%M]);
      if( equals(cross(ururu_beam.p1-ururu_beam.p2,seg.p1-seg.p2),0.0) ) continue;
      seg = Segment(polies[i][j],polies[i][(j+1)%M]);

      if( !intersectSS(ururu_beam,seg) ) continue;

      Point cp = crosspoint(ururu_beam,seg);

      //交点が src or dst である場合必ず正しい移動
      if( cp == ururu_beam.p1 || cp == ururu_beam.p2 ) continue;

      if( cp == seg.p1 || cp == seg.p2 ) {

	// p0 -> p1 ->p2, p0 -> p1 -> p3 が平行であることはない ( 前でcontinueしてるから )
	Point p0 = src;
	Point p1 = cp;
	Point p2 = polies[i][(j+1)%M];
	Point p3 = polies[i][(j-1+M)%M];
	if( cp == p2 || cp == p3 ) continue;
	int ccw_res1 = ccw(p0,p1,p2);
	int ccw_res2 = ccw(p0,p1,p3);
	if( ccw_res1 != ccw_res2 ) return false;
      } else return false;
    }
    if( !( src == bomb || dst == bomb ) && !on && in ) return false;
  }


  /*
  if( !( src == bomb || dst == bomb ) ) {
    Point mp = ( src + dst ) / 2.0;
    //rep(i,N) if( inPolygon(polies[i],mp,false) ) return false;
    rep(i,N) if( contains(polies[i],mp) == IN ) return false;
  }
  */

  return true;
}

double temper;
void dijkstra(){
  temper = DINF;
  priority_queue<Data> Q;
  Q.push((Data){(Point){0,0},0,-1,-1});
  while( !Q.empty() ){  
    Data data = Q.top(); Q.pop();
    if( data.a != -1 && LT(mindist[data.a][data.b],data.weight) ) continue;
    if( data.a != -1 && goals[data.a][data.b] ) {
      temper = data.weight;
      return;
    }
    rep(i,N) rep(j,polies[i].size()){
      Point next = polies[i][j];
      if( data.cur == next ) continue;
      if( !isValidMove(data.cur,next) ) continue;
      double next_weight = data.weight + abs(next-data.cur);
      if( LT(next_weight,mindist[i][j]) ) {
	mindist[i][j] = next_weight;
	Q.push((Data){next,next_weight,i,j});
      }
    }
  }
}

//問題の制約より、ururu_beamとsegが平行であることはない
bool isGoal(Point p){
  Segment ururu_beam = Segment(bomb,p);
  rep(i,N){
    int M = polies[i].size();
    rep(j,M){
      Segment seg = Segment(polies[i][j],polies[i][(j+1)%M]);
      if( !intersectSS(ururu_beam,seg) ) continue;
      Point cp = crosspoint(ururu_beam,seg);

      if( !( cp == p ) ) return true;
      else if( cp == seg.p1 || cp == seg.p2 ) {
	int k = (cp==seg.p1?j:(j+1)%M);
	int c1 = ccw(bomb,p,polies[i][(k-1+M)%M]);
	int c2 = ccw(bomb,p,polies[i][(k+1)%M]);
	assert( !( p == polies[i][(k-1+M)%M]));
	assert( !( p == polies[i][(k+1)%M]));
	assert( c1 == COUNTER_CLOCKWISE || c1 == CLOCKWISE );
	assert( c2 == COUNTER_CLOCKWISE || c2 == CLOCKWISE );
	if( c1 == c2 ) return true;
      }
    }
  }
  return false;
}

//爆弾と線分は平行ではないことが保証されている
Point getOnSegmentPoint(Point p){
  Segment ururu_beam = Segment(bomb,p);
  Vector e = (ururu_beam.p2-ururu_beam.p1) / abs(ururu_beam.p2-ururu_beam.p1);
  ururu_beam.p2 = bomb + e * 1000000;
  Point ret = Point(-DINF,-DINF);
  rep(i,N) {
    int M = polies[i].size();
    rep(j,M){
      Segment seg = Segment(polies[i][j],polies[i][(j+1)%M]);
      if( !intersectSS(ururu_beam,seg) ) continue;
      Point cp = crosspoint(ururu_beam,seg);
      if( cp == p ) continue;
      if( ret == Point(-DINF,-DINF) ) ret = cp;
      else if( LT(abs(cp-bomb),abs(ret-bomb)) ) ret = cp;
    }
  }
  return ret;
}

Point vir;
bool comp_dist(const Point &p,const Point &q) {
  return LT(abs(vir-p),abs(vir-q));
}

void compute(){


  if( isGoal(Point(0,0)) ) {
    printf("%.10f\n",0.0);
    return;
  }

  rep(i,N) rep(j,polies[i].size()) goals[i][j] = isGoal(polies[i][j]);

  // dijkstra で各頂点間の最短距離を求める
  dijkstra();


  //rep(i,N) rep(j,polies[i].size()) goals[i][j] = isGoal(polies[i][j]);
  //double answer = DINF;
  double answer = temper;
  rep(i,N) rep(j,polies[i].size()) if( goals[i][j] && LT(mindist[i][j],answer) ) {
    answer = mindist[i][j];
  }

  // ゴールとなる候補店の列挙
  // 各点について、ゴールか、ゴールであれば最小値を更新する
  rep(i,N) {
    int M = polies[i].size();
    rep(j,M){
      Point p0 = bomb;
      Point p1 = polies[i][j];
      Point p2 = polies[i][(j+1)%M];
      Point p3 = polies[i][(j-1+M)%M];
      int ccw_res1 = ccw(p0,p1,p2);
      int ccw_res2 = ccw(p0,p1,p3);
      if( ( ccw_res1 == COUNTER_CLOCKWISE || ccw_res1 == CLOCKWISE ) && ( ccw_res2 == COUNTER_CLOCKWISE || ccw_res2 == CLOCKWISE ) && ccw_res1 != ccw_res2 ) {
	continue;
      }

      //爆弾は多角形を構成する線分の延長線上にはない
      //線分上
      Point p = getOnSegmentPoint(polies[i][j]);
      if( p.x != -DINF && p.y != -DINF ) {
	//if( !( LT(abs(p),answer) ) )
	if( isGoal(p) )	{
	  if( LT(abs(p),answer) && isValidMove(Point(0,0),p) ) {
	    answer = min(answer,abs(p));
	  } else {
	    double mini = DINF;
	    rep(k,N) rep(l,polies[k].size()) if( LT(mindist[k][l]+abs(polies[k][l]-p),min(answer,mini)) && isValidMove(polies[k][l],p) ) {
	      mini = min(mini,mindist[k][l]+abs(polies[k][l]-p));
	    }
	    answer = min(answer,mini);
	  }
	}
      }

      //垂直
      p = projection(Segment(bomb,polies[i][j]),Point(0,0));

      //if( !( LT(abs(p),answer) ) )
      if( isGoal(p) )	{
	if( isValidMove(Point(0,0),p) ) {
	  answer = min(answer,abs(p));
	} else {
	  bool success = true;
	  Segment ururu_beam = Segment(p,Point(0,0));
	  rep(k,N) {
	    int M2 = polies[k].size();
	    rep(l,M2){
	      Segment seg = Segment(polies[k][l],polies[k][(l+1)%M2]);
	      if( !intersectSS(ururu_beam,seg) ) continue;
	      if( equals(cross(ururu_beam.p1-ururu_beam.p2,seg.p1-seg.p2),0.0) ) continue;
	      Point cp = crosspoint(ururu_beam,seg);
	      if( cp == seg.p1 || cp == seg.p2 || cp == ururu_beam.p1 || cp == ururu_beam.p2 ) continue;
	      success = false;
	      break;
	    }
	    if( !success ) break;
	  }
	  if( success ) {
	    double mini = DINF;
	    rep(k,N) rep(l,polies[k].size()) if( goals[k][l] ) {
	    mini = min(mini,mindist[k][l]+abs(ururu_beam.p1-ururu_beam.p2));
	    }
	  }
	}
      }
    }
  }

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

int main(){
  while( scanf("%d",&N) , N ) {
    polies.clear(), mindist.clear(); goals.clear();
    polies.resize(N), mindist.resize(N), goals.resize(N);
    //cin >> bomb.x >> bomb.y;
    scanf("%lf %lf",&bomb.x,&bomb.y);
    rep(i,N) {
      int m;
      scanf("%d",&m);
      polies[i].resize(m), mindist[i].resize(m,DINF), goals[i].resize(m,false);
      rep(j,m) scanf("%lf %lf",&polies[i][j].x,&polies[i][j].y);//cin >> polies[i][j].x >> polies[i][j].y;
    }
    compute();
  }
  return 0;
}