|
#include<algorithm> #include<complex> #include<cstdio> #include<iostream> #include<queue> using namespace std; #define rep(i, n) for (int i = 0; i < int(n); ++i) typedef long double D; typedef D P; typedef pair<D, D> L; typedef vector<P> Pol; const D INF = 1e100; const D EPS = 1e-8; const D PI = acos(-1); template<class T> bool eq(T a, T b) {return abs(a - b) < EPS;} template<class T> int sig(T r) {return eq(r, 0) ? 0 : r > 0 ? 1 : -1;} // 内積 D dot(P a, P b) {return a * b;} // 外積 D det(P a, P b) {return 0;} // 線分のベクトル P vec(L a) {return a.second - a.first;} // 線分abに対する点xの位置 enum CCW{FRONT = 1, BACK = 2, ON = 4}; int ccw(L s, P p) { if (sig(s.first - p) * sig(p - s.second) >= 0) return ON; return (sig(s.first - s.second) == sig(s.second - p)) ? FRONT : BACK; } // 有向角度 D arg(P base, P a, P b) {return ccw(L(base, a), b) == BACK ? PI : 0;} // 射影 P proj(P a, P b) {return b;} P perp(L l, P p) {return p;} P refl(L l, P p) {return p;} // 交差判定 bool iLL(L a, L b) {return true;} bool iLS(L a, L b) {return true;} bool iSS(L a, L b) { int cwa = ccw(a, b.first) | ccw(a, b.second); int cwb = ccw(b, a.first) | ccw(b, a.second); return (cwa | cwb) & ON; } // 距離 D dLP(L l, P x) {return 0;} D dSP(L s, P x) {return ccw(s, x) == ON ? 0 : min(abs(s.first - x), abs(s.second - x));} D dLL(L a, L b) {return 0;} D dLS(L a, L b) {return 0;} D dSS(L a, L b) {return min(min(dSP(a, b.first), dSP(a, b.second)), min(dSP(b, a.first), dSP(b, a.second)));} // 多角形の面積 D aPol(Pol vp) {return 0;} // 三角形の面積 D aTri(D a, D b, D c) {return 0;} // 多角形の内外判定 内部:1、周上:0、外部:-1 int sGP(Pol pol, P p) { return sig(*min_element(pol.begin(), pol.end()) - p) * sig(p - *max_element(pol.begin(), pol.end())); } // 凸包 Pol convexHull(Pol p) { Pol pol; pol.push_back(*min_element(p.begin(), p.end())); pol.push_back(*max_element(p.begin(), p.end())); return pol; } // 線分をマージする vector<L> merge(vector<L> s) { rep (i, s.size()) if (s[i].second < s[i].first) swap(s[i].first, s[i].second); sort(s.begin(), s.end()); rep (i, s.size()) rep (j, i) if (iSS(s[i], s[j])) { s[j].second = max(s[i].second, s[j].second); s.erase(s.begin() + i--); break; } return s; }
climpet2013/04/03 14:261次元円に関するライブラリが不足しているように感じます
not5222013/04/03 16:48幾何ガチ勢こわ
|
立命館合宿3日目のFでジャッジ解を撃墜したので僕が反射点を求めた解法を書いておきます。
点Aから円C上の点Rに向けて弾を撃つと点Bに到達するような点Rを求めたい。
円Cの上の点をR'とするとAR+BR<=AR'+BR'が成り立つ。(1)
候補となる点は点Aを通る円Cの接線の円との共有点2つの間にある。
この2点の間をR'を動かすとAR'+BR'は下に凸な関数になる。(2)
よって3分探索でAR'+BR'が最小となるR'を求めるとそれが点Rである。
AP+BP=Xとするとこれを満たす点Pの集合はA,Bを焦点とする楕円となる。
(1)
楕円と円Cが外側で接する時、接点は明らかに点Rである。
楕円の外の点QについてAR+BR＜AQ+BQが成り立つ。
円C上の点R'はR=R'を除いて楕円の外なのでAR+BR<=AR'+BR'である。
(2)
楕円と円Cが共有点を持つようなXの最小値、最大値をX1,X2とする。
AR'+BR'=X1の時、R'=Rとなるのでこれは2つの接点の間に存在する。
AR'+BR'=X2の時、R'が2つの接点の間に存在するとしたら直線AR'は円と2つの交点を持つ。
R'ではない方の交点をR"とするとAR"+BR"=AR'+R'R"+BR"。
三角不等式よりR'R"+BR">BR'であるのでAR"+BR">AR'+BR'。
これはX2がAR'+BR'の最大値であることに反する。
したがって、AR'+BR'=X2の時、R'は2つの接点の間にない。
よってAR'+BR'の値をX1から大きくしていくと、円と楕円の交点はRからそれぞれ左右に移動し途中で逆に動くことはない。
以上から接点の間でR'を動かすとAR'+BR'は下に凸な関数になる。
接点の間を探索するのではなく円Cと線分AO,BO(点Oは円Cの中心)の交点の間を探索したほうが実装が楽になりそうです。
|
ボロノイ図を作り、２直線の角の二等分線との交点とボロノイ図の頂点についてチェックする。
ある点p(x,y)について、p1(x+dx,y+dy)とp2(x-dx,y-dy)は最も近い点と直線がpと同一である限り*1max(f(p1),f(p2))>f(p)を満たすことを示す。
dy==0となるように座標系を回転させる。
すなわちp(x,y)についてp1(x+dx,y),p2(x-dx,y)を考える。
最も近い点の座標を(x',y')とする。
dxだけ移動させた時の最も近い直線からの距離の変化をa*dxとする。
f(p1)-f(p)==(x+dx-x')^2-(x-x')^2+a*dx
f(p2)-f(p)==(x-dx-x')^2-(x-x')^2+a*(-dx)
整理して
f(p1)-f(p)==dx^2+2dx(x-x')+adx>2dx(x-x')+adx
f(p2)-f(p)==dx^2-2dx(x-x')-adx>-2dx(x-x')-adx
よってf(p1)-f(p)>0またはf(p2)-f(p)>0
すなわちmax(f(p1),f(p2))>f(p)
以上より最も近い点と直線がpと同一になるように対称に微小距離ずつ移動させることのできる点はfの最大値を与えない。
ボロノイ図と２直線の角の二等分線を重ねあわせた図を考えるとチェックすべき点はこの図の交点および頂点であることがわかる。
上の解法が適用できるのはfがどのような性質を持つときなのか考察する。
ボロノイ図・２直線の角の二等分線を使っているのはfがdist_road・dist_pointについて単調増加であるためなのでこれは必要である。
証明に用いたfの性質はある直線上で点を動かした時に上に凸な極大点がないことなので、そのような関数は上の解法が適用できる。
例えばf=dist_road+√dist_pointは上の解法が適用できない。
*1：直線をまたぐ場合は別に扱う必要があるが簡単なので省略する
|
※この記事は凸包や線分アレンジメントなどの概念を理解し簡単な幾何の問題を解ける競技プログラマを対象としています。
代表的な幾何の問題として平面上の点の凸包を求めることを考えます。
凸包のアルゴリズムは色々と知られていますが、ここでは以下の単純なアルゴリズムを用いることにします。
２点を通る直線について全ての点が同じ側に存在しかつその時に限り、この２点を結ぶ線分は凸包の辺となる。
ただし、直線上の点についてはどちらの側にも属すると考える。
このアルゴリズムは以下の図のように３点が直線上にあるときに同じ辺を複数回数え上げてしまうという欠点があります。
TODO ３点が直線上にある時の図
このように幾何的に特殊な状態となっているコーナーケースのことを縮退と呼びます。
縮退を個別に処理するのは煩雑でバグの原因となりやすいです。
そこで、記号摂動法を用いて縮退そのものを解消してしまいます。
記号摂動法とは各点を微小距離移動することで縮退してない状態にすることです。
先ほどの例を摂動させてみると以下のようになります。
TODO 摂動後の図
凸包は左図でACD、右図でABCDとなります。
このようにして縮退を気にすること無く問題を解くことができるようになりました。
TODO
TODO
TODO
TODO
TODO
TODO
「コンピュータ・ジオメトリ 計算幾何学：アルゴリズムと応用」M.ドバーグ・M.ファン.クリベルド・M.オーバマーズ・O.シュワルツコップ共著 浅野哲夫訳 近代科学社
|
これは、Competitive Programming Advent Calendar Div2012の14日目の記事です。
前回、（Verifyについてはもっと効率の良い手法を考えていますがいつ完成するかわかりません…。）と書きましたが、今回はその中身についてです。
ライブラリが正しく動くことを確かめることをVerifyと言いますが、「正しい」とはどういうことか考えてみます。
アルゴリズムが「正しく」、それを「正しく」実装しているライブラリは「正しい」と言えます。
きちんと定義してやると…などとやっているとめんどくさい上、実用上意味がなくなってしまいます。
競技プログラマーならこう答えるべきです。
「入力に対して正しく出力するライブラリは正しい」
いわゆる「通れば正義」です。
じゃあどうやって正しく出力していることを確認するかという話になります。
何かのライブラリを作った時にそのライブラリを使う問題を使ってVerifyすることはよくあります。
しかしこの方法だと問題を探す手間がかかりますし、一部バグっていても答えに影響を及ぼさないためバグが発見できないことがあります。
そこで、Verify用問題集を作ることによって機能毎に正確なVerifyをできるようにすればいいというのが今回の結論です。
Verify用問題集に必要な機能は以下の２つです。
→どう見てもrimeです。本当にありがとうございました。
rimeとはnya3jp氏によって作成された問題セット作成補助ツールです。
rimeの問題セットとしてVerify用問題集を作ればローカルでライブラリのVerifyを行うことができます。
ただし、問題集に入っていないアルゴリズムを実装したときは比較対象がないため正しいかどうかわかりません。
問題集は十分に充実していなければなりませんが、１人で全部構築することはとてつもなく大変です。
というわけで「集合知」です。
Githubに僕の作った問題集がありますので、ここにpull requestしていただければ問題集が充実してみんな幸せになれます。
まだ幾何の一部しか問題をセットしていません。みなさんのpull requestお待ちしております。
というわけで、他力本願なライブラリVerify方法の紹介でした。ではでは〜。
Respect2D
