AGC021 ~B問題 Holes ~ - banbooooのブログ

初投稿です。個人的にはアルゴリズムも実装も大変だったので書き残しておきます。
agc021.contest.atcoder.jp

乱数で決まる $(x,y)$ は半径 $R$ の円の内部。穴は1辺 $2 \times 10^6$ の正方形の内部の領域に存在する。
$R$ がめちゃめちゃ大きいので、穴の近くには滅多に落ちない。

f:id:banboooo:20180226021502j:image:w400 — 全体図

穴を頂点とする凸多角形で、かつ全ての穴を内部含むような図形(凸包というらしい)を作る。乱数で決まる点は主に、この凸多角形の外側の領域に落ちるので、凸多角形の頂点になってる穴のどれかに落ちる。

内部に落ちた時のみ、頂点に使われてない穴にも落ちる可能性があるが、確率は0とみなせる。(精度的に)つまり、凸多角形の頂点に対してのみ確率を考えればよい !!

隣同士の2つの穴について考えると、2つの穴を結んだ直線の垂直二等分線を基準に、どちらの穴に落ちるかがわかるので、これを凸多角形の全ての辺に対してやってみる。

すると、全部の垂直二等分線が中心に集まり、1点で交わるわけではないので小さい領域ができてしまうが、これらを無視すれば、下の図のように領域分割すれば良さそう。

f:id:banboooo:20180226020127j:image:w400 — 領域分割

中心角 = (垂直二等分線の交わる角度) = $\pi$ - (辺の交わる角度)
例: 赤い領域の中心角 = $\pi - \angle BAE$

アルゴリズムはここまでだが、実装でも苦労した点が2点ほど。

そもそも凸包はどうやって求めるのか。

凸包を求めるアルゴリズムで「グラハムスキャン」というのがあるらしい (蟻本p224)
2辺のベクトル( $\vec{v1},\vec{v2}$ )の符号付き角度が、 $det(v1,v2)$ の正負に一致するので、これが常に正になるように点を反時計回りに選ぶことがポイント。

長さの値をどうやって角度の値に直すか。

最初は $\arctan$ を使って角度に変換することを考えたが、場合分けが必要でめんどくさい。
そこでベクトルをおいて、 $\theta = \arccos{\dfrac{\vec{a} \vec{b}}{|a||b|}}$ で求めると楽だった。

ソースコード
まず二次元ベクトルの内積やノルムなどを簡単に求められるようにクラスを構築しておく。
なお、今回の問題では点ごとに確率を出力しないといけないので、indexの値を保持できるようにしておき、辺のベクトルを表すときは、indexの値は(-1)にしておいた。

typedef struct vec2D {
  double x, y;
  long ind;
  vec2D() {}
  vec2D(double x, double y,long index) : x(x), y(y) ,ind(index) {} 
  vec2D operator + (vec2D p){
	   return vec2D(x + p.x, y + p.y,-1);
  }
  vec2D operator - (vec2D p){
	   return vec2D(x - p.x, y - p.y,-1);
  }
  vec2D operator * (double d){
	   return vec2D(d*x, d*y,-1);
  }

  double norm(void){
  	return sqrt(x*x + y*y);
  }

  double twiceNorm(void){
  	return x*x + y*y;
  }

  double dot(vec2D p){ //pとの内積
	   return x * p.x + y * p.y;
  }
  double det(vec2D p){ // (x,p)の行列式
	   return x * p.y - p.x * y;
  }
  double dist(vec2D p){ //pとの距離の2乗
	   return (x-p.x)*(x-p.x) + (y-p.y)*(y-p.y);
  }
  bool operator <(const struct vec2D &e) const{
     return x == e.x? (y < e.y) : x < e.x;
  }
  bool operator >(const struct vec2D &e) const{
     return x == e.x? (y > e.y) : x > e.x;
  }
} pt;

次に、多角形クラスを定義しておく。これでまた凸包を求める問題が出ても怖くない !
なお、グラハムスキャンのところはほぼ蟻本写してます。

bool cmp_x(const pt &p1,const pt &p2){ //ソート用の比較関数
	if (p1.x != p2.x) return p1.x < p2.x;
    return p1.y < p2.y;
}
class Polygon{
  private:
    vector<pt> pts;
    size_t n;
  public:
    Polygon() : n(0) {}
    Polygon(vector<pt> v) : n(v.size()) {
      pts = v;
    }

    void addPoint(pt P){
      pts.push_back(P);
    }

    vector<pt> convexHull(void){ // 凸包を構成する点の配列を返す。
      vector<pt> qs(n * 2);
      sort(pts.begin(),pts.end(),cmp_x);
      long k = 0;

      REP(i,n){ //下側凸包の作成
        while(k > 1 && (qs[k-1] - qs[k-2]).det(pts[i] - qs[k-1]) <= 0) k--;
        qs[k++] = pts[i];
      }

      for(long i = n - 2, t = k; i >= 0; i--){ //上側凸包の作成
        while(k > t && (qs[k-1] - qs[k-2]).det(pts[i] - qs[k-1]) <= 0) k--;
        qs[k++] = pts[i];
      }
      qs.resize(k-1);
      n = k-1;
      pts = qs;
      return qs;
    }
};

最後に、main関数を書く

int main(void){
	short N;
	double x,y;
	cin >> N;
	vector<pt> point,convex;

	REP(i,N) {
		cin >> x >> y;
		point.push_back(pt(x,y,i));
	}

	if (N == 2) { // コーナーケース
		cout << 0.5 << "\n" << 0.5 << endl;
		exit(0);
	}

	Polygon P(point);
	convex = P.convexHull(); //凸包を求める
	vector<double> ans(N,0.0); //答えの確率を格納
	long m = convex.size();

	double alpha;
        // 角度を計算。
	REP(i,m-2){//凸包の点1 ~ (m-2)までの確率を計算
		pt v1 = convex[i]-convex[i+1],v2 = convex[i+2]-convex[i+1];
		alpha = acos((v1).dot(v2) / sqrt(v1.twiceNorm() * v2.twiceNorm()));
		ans[convex[i+1].ind] = (M_PI - alpha)/ (2 * M_PI);
	}
        // 凸包の点m-1の確率
	pt v1 = convex[m-2] - convex[m-1],v2 = convex[0] - convex[m-1];
	alpha = acos((v1).dot(v2) / sqrt(v1.twiceNorm() * v2.twiceNorm()));
	ans[convex[m-1].ind] = (M_PI - alpha)/ (2 * M_PI); 
        // 凸包の点0の確率
	v1 = convex[m-1] - convex[0],v2 = convex[1] - convex[0];
	alpha = acos((v1).dot(v2) / sqrt(v1.twiceNorm() * v2.twiceNorm()));
	ans[convex[0].ind] = (M_PI - alpha)/ (2 * M_PI);

	REP(i,N){
		cout << ans[i] << endl;
	}

	return 0;
}