AOJ 1151 Twirl Around(くるくる)

AOJ 幾何

問題

日本語なので本文参照（http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=1151&lang=jp）


多角形の内側にそって線分を回転させて、2*PI*Rラジアンまわした後の座標を答える。
棒がつっかかったらそこで終了。

方針

シミュレートする。
線分が次にぶつかる場所を求め、回転してそこまでぶつけて、
そこを中心にしてまた次にぶつかる場所を……と次々にやっていけばよい。


0距離でぶつかるようになってしまったら、つっかえているということなので終了。
ただし、離れる向きで最初から接触している分には問題なくて、
めり込む向きに最初から接しているのだけがダメ。


次にぶつかる場所はどうやって求めるかというと、

  • 回している線分の端点が多角形の辺にぶつかる場合
  • 回している線分の途中が多角形の頂点にぶつかる場合

の両方の場合が考えられて、前者は円と線分の交点で求まり、後者は点と直線の距離を使えばよい。


交点がもとまったら、ぶつかる角度の順にソートして、
一番小さいものが次にぶつかる点。同じ角度が複数あったら一番遠いところ。
で、最初からくっついている場合があるが、離れる向きのものはセーフ（離れるだけなので無視してよい）で、めりこむ向きのものがあったら、つっかかるのでそこでループ打ち切り。


離れる向きかめり込む向きかは、上の前者の場合と後者の場合で判定方法が違って、

  • 前者なら、(回転中心)-(接点)-(多角形の接点の次の頂点)が時計回りならめりこむ
  • 後者なら、微小な角度だけ回転させて、内部から出てしまったらアウト。


バグったら頑張って半自動で図示するビジュアライザー作って確認する（した）。

ソースコード

long double len, R;
int n;
 
inline double arg_(const P &p){
	double a = arg(conj(p));
	while(a < -EPS) a += 2 * PI;
	return a;
}
inline void rot(P &p, const P &c, double a){
	p = P(cos(a), sin(a)) * (p - c) + c;
}
 
int main(){
	while(cin >> len >> R >> n, n){
		R *= 2 * PI;
		 
		G g;
		rep(i, n){
			double x, y; cin >> x >> y;
			g.pb(P(x, y));
		}
		L bar(P(0, 0), P(0, len));
		P c(0, 0);
		 
		 
		int loop = 0;
		 
		while(R > EPS){
			//dbg2(bar[0], bar[1]);
			 
			long double r0 = abs(bar[0] - c);
			long double r1 = abs(bar[1] - c);
			 
			vector<pair<long double, P> > p0;
			vector<pair<long double, P> > p1;
			P next(c);
			long double nd = 0;
			 
			//頂点との交わり
			rep(i, n){
				long double d = abs(g[i] - c);
				if(d < EPS) continue;
				if(d - EPS < r0){
					 
					long double a = arg_((g[i] - c) / (bar[0] - c));
					if(a > EPS) p0.pb(mp(a, g[i]));
					else if(ccw(c, g[i], g[(i + 1) % n]) == -1){
						R = -1; break;
					}
				}
				if(d - EPS < r1){
					 
					long double a = arg_((g[i] - c) / (bar[1] - c));
					if(a > EPS) p1.pb(mp(a, g[i]));
					else if(ccw(c, g[i], g[(i + 1) % n]) == -1){
						R = -1; break;
					}
				}
			}
			
			//辺との交わり
			rep(i, n){
				const P &p = g[i];
				const P &q = g[(i + 1) % n];
				G g0 = isCL(c, r0, p, q);
				G g1 = isCL(c, r1, p, q);
				 
				if(r0 > EPS) rep(j, g0.size()){
					if(ccw(g0[j], p, q) != 2) continue;
					 
					long double a = arg_((g0[j] - c) / (bar[0] - c));
					if(a > EPS) p0.pb(mp(a, g0[j]));
					else{
						rot(g0[j], c, -.00005);
						if(contains(g, g0[j]) != IN){
							R = -1; break;
						}
					}
					 
				}
				if(r1 > EPS) rep(j, g1.size()){
					if(ccw(g1[j], p, q) != 2) continue;
					 
					long double a = arg_((g1[j] - c) / (bar[1] - c));
					if(a > EPS) p1.pb(mp(a, g1[j]));
					else{
						rot(g1[j], c, -.00005);
						if(contains(g, g1[j]) != IN){
							R = -1; break;
						}
					}
				}
			}
			if(R < EPS) break;
			if(!p0.size() && !p1.size()) break; //??????
			sort(all(p0)); sort(all(p1));
			 
			p0.insert(p0.end(), all(p1));
			sort(all(p0));
			 
			rep(i, p0.size()) if(abs(p0[i].first - p0[0].first) < EPS){
				long double d = abs(p0[i].second - c);
				if(d > nd + EPS){
					nd = d;
					next = p0[i].second;
				}
			}
			 
			long double ang = min(p0[0].first, R);
			R -= ang;
			rot(bar[0], c, -ang);
			rot(bar[1], c, -ang);
			c = next;
		}
		printf("%.9f %.9f\n", (double)bar[0].real(), (double)bar[0].imag());
	}
	 
	return 0;
}