SRM 201 DIV1 Hard DogWoods
問題
ある性質をもって走る犬がいる。
- 木にぶつかるまでの間、点(0, 0)を中心とした円周上を時計回りに走る
- 木にぶつかった場合、上記の動きができるようになる位置まで、木にそって反時計回りに走る
犬が点(startx, starty)から走り始める時、点(0, 0)からの距離が10以下になるまでの、走る距離を求めよ。
- 犬は点と考える
- 木は円として複数(<=50)与えられる
解答
がんばって式を解いてシミュレーションする。
二つの円周の交差する点は問題文に与えられているのでそれを使った。
「木にぶつかってから、時計回りの動きができるようになる点」はWolframAlphaで適当に解いて、適当に角度のminを取ったりした
コメント
普通に式を解けばいいんだから、WolframAlphaさんにやってもらえば普通にできるよね。
頑張れば書けるんだからがんばろう！
いちいち式を解くのに複雑になっていて、複素数とかでもっと簡単な記述ができるものもあるだろうから、幾何とか微分とか数学の基本をやっぱりきちんと勉強したほうがいいと思う
コード
#include <string> #include <vector> #include <algorithm> #include <numeric> #include <set> #include <map> #include <queue> #include <iostream> #include <sstream> #include <cstdio> #include <cmath> #include <ctime> #include <cstring> #include <cctype> #include <complex> #define rep(i,n) for(int (i)=0;(i)<(int)(n);++(i)) #define each(it,o) for(__typeof((o).begin()) it = (o).begin(); it != (o).end(); ++ it) #define pb(x) push_back(x) #define mp(x,y) make_pair((x),(y)) #define EPS 1e-8 using namespace std; //ばっど #define double long double struct DogWoods { typedef complex<double> P; double sq(double x) { return x*x; } //剰余だが、%演算子やfmodはmodでなくremのほう(Haskellの場合)なので、modを実装する //QuickCheckでこの定義でいいらしいとわかった //角度はmodがつきものなので便利 double fmod2(double x, double y) { if(x > 0) return fmod(x, y); else return fmod(y - fmod(abs(x), y), y); } double howFar(vector <int> xs, vector <int> ys, vector <int> ds, int startx, int starty) { const double pi = acos(-1.); int n = xs.size(); P c = P(startx, starty); double total = 0; if(abs(c) < 10+EPS) return 0; rep(ii, 100) { double r1 = abs(c); double ca = -arg(c); vector<pair<int,P> > v; //当たりそうな木を列挙する rep(i, n) { double x0 = xs[i], y0 = ys[i], r0 = ds[i] / 2.; double d = sq(x0) + sq(y0); double pp = (sq(r0 + r1) - d)*(d - sq(r1 - r0)); //sqrtの中身がマイナス<->解無し if(pp < -EPS) { continue; } if(pp < 0) pp = 0; //計算式、(見た目が)ひどいなあ！ double p = sqrt(pp); double tx0 = x0 / 2 - x0*(sq(r0) - sq(r1)) / (2*d) + (y0*p) / (2*d), ty0 = y0 / 2 - y0*(sq(r0) - sq(r1)) / (2*d) - (x0*p) / (2*d), tx1 = x0 / 2 - x0*(sq(r0) - sq(r1)) / (2*d) - (y0*p) / (2*d), ty1 = y0 / 2 - y0*(sq(r0) - sq(r1)) / (2*d) + (x0*p) / (2*d); v.pb(mp(i, P(tx0, ty0))); if(p >= EPS) v.pb(mp(i, P(tx1, ty1))); } double ma = 999; P mp; int mi = -1; each(i, v) { P p = i->second; //どれくらいの角度走るかを計算 double a = (-arg(p)) - ca; a = fmod2(a, 2*pi); //今いたところっぽかったら無視 if(a < EPS || 2*pi-a < EPS) continue; if(ma > a) ma = a, mp = p, mi = i->first; } if(mi == -1) return -1.0; total += ma * r1; c = mp; //atan(x-x0, y-y0) = atan(x, y) { double x0 = xs[mi], y0 = ys[mi], r = ds[mi] / 2.; double b = arg(c - P(x0, y0)); double d = sq(x0) + sq(y0); double pp = (sq(r + 10) - d)*(d - sq(r - 10)); if(abs(pp) < EPS) pp = 0; if(pp >= 0) { double p = sqrt(pp); double ux0 = x0 / 2 - x0*(sq(r) - sq(10)) / (2*d) + (y0*p) / (2*d), uy0 = y0 / 2 - y0*(sq(r) - sq(10)) / (2*d) - (x0*p) / (2*d), ux1 = x0 / 2 - x0*(sq(r) - sq(10)) / (2*d) - (y0*p) / (2*d), uy1 = y0 / 2 - y0*(sq(r) - sq(10)) / (2*d) + (x0*p) / (2*d); double a0 = arg(P(ux0, uy0) - P(x0, y0)) - b, a1 = arg(P(ux1, uy1) - P(x0, y0)) - b; a0 = fmod2(a0, 2*pi); a1 = fmod2(a1, 2*pi); double aa = min(a0, a1); total += aa * r; break; } //これも見た目がやばい。式は適当だけど角度のminを取ったらうまくいった double t = sqrt(sq(r)*sq(x0)*sq(y0) + sq(r)*sq(sq(y0))); double tx0 = -((x0*t) / (y0*d)) + (x0*sq(y0)) / d + (x0*x0*x0) / d, ty0 = -(t / d) + (sq(x0)*y0) / d + (y0*y0*y0) / d, tx1 = ((x0*t) / (y0*d)) + (x0*sq(y0)) / d + (x0*x0*x0) / d, ty1 = (t / d) + (sq(x0)*y0) / d + (y0*y0*y0) / d; P tp0 = P(tx0, ty0), tp1 = P(tx1, ty1); double a0 = arg(tp0 - P(x0, y0)) - b, a1 = arg(tp1 - P(x0, y0)) - b; a0 = fmod2(a0, 2*pi); a1 = fmod2(a1, 2*pi); double a; P tp; if(a0 < a1) a = a0, tp = tp0; else a = a1, tp = tp1; total += a * r; c = tp; } } return total; } }; #undef double
