SRM 215 DIV1 Hard ShortCut

問題 Editorial

問題

平面上に一本の道がある。これは終点が次の始点となる有向線分(roadX,roadY)をつなげた形で、交差・接触はしない。
この道の始点から終点までを最小時間で移動したい。
道にそって走る時は単位距離に対して単位時間かかる。ただし道の向きの通りに走らなければいけない。
また、道をずれて走ることもできる。この場合は道を通る時と比べて2倍の時間がかかってしまう(1/2の速度で走ることが出来る)。
このとき(最小時間/道の総距離)を求めよ

解答

Editorial・rng_58さんの解答を読んで参考にした。
まず、移動する際に曲がることがある点は限られており、列挙できる。
これは詳しくはEditorialを読む・最後の図を見る方がいいが、
全ての(s,t)∈(道の線分の両端)と頂点p∈(全ての線分の端点)に対して(
線分(s,t)だけを使えてpからtに最短時間で向かう時、その移動経路が線分(s,t)に最初に入る点 ∪
線分(s,t)だけを使えてsからpに最短時間で向かう時、その移動経路で線分(s,t)から最初に出る点 )
となる。
線分が一つの時の最短経路で線分に(入る・出る)点は凸(の逆さ)関数になるので、その点は微分しても求められるし黄金探索しても良い。
その点たちを求めたらあとはダイクストラとかをすれば良い。
(同じ線分上であって向きが合っているなら1・else 2)*点の距離を重みとする。
しかし重みを2500*2500はMLEしてしまう。しかしこれを解決する上手く美しい方法がある。
まず点たちを列挙する時に、対象となる線分上の距離((s,t)に対してsからの距離)でソートする。
またそれを、使っている線分のインデックスでソートする。
こうすると必ず、本来の道の通りの順番ということになる。
そして重みを求める際、二つの(v,u)があるとして、((v+1 == u)なら1・else 2)*点の距離、というふうに求められる。
これは考えればわかるが、そのような二点は必ず道の通りに進んでいけるからである

コメント

色々とよくわからないところがあってかなり時間を使ってしまった。しかしそういうのが(解答を知ってても)コードを書けるといいね。そのおかげで今は理解できているし

コード

#include <vector>
#include <algorithm>
#include <queue>
#include <cmath>
#include <cstring>
#include <complex>
#include <limits>
#define rep(i,n) for(int (i)=0;(i)<(int)(n);++(i))
#define all(o) (o).begin(), (o).end()
#define pb(x) push_back(x)
#define mp(x,y) make_pair((x),(y))
#define mset(m,v) memset(m,v,sizeof(m))
using namespace std;
typedef complex<double> P;
struct ShortCut {
	double suvTime(vector <int> roadX, vector <int> roadY) {
		int n = roadX.size();
		double standardDist = 0;
		rep(i, n-1) standardDist += hypot(roadX[i]-roadX[i+1], roadY[i]-roadY[i+1]);
		const double inf = numeric_limits<double>::infinity();
		vector<double> pointsx, pointsy;
		rep(i, n-1) {
			P s0 = P(roadX[i], roadY[i]), s1 = P(roadX[i+1], roadY[i+1]);
			vector<double> mids;
			double len = abs(s1 - s0);
			P cs;
			rep(dir, 2) {
			swap(s0, s1);
			cs = P(cos(arg(s1-s0)), sin(arg(s1-s0)));
			rep(j, n) {
				P p = P(roadX[j], roadY[j]);
				double a = 0, b = len;
				#define DIST(x) 2*abs((s0 + x*cs) - p) + (len - x)
				const double ratio = 2 / (3 + sqrt(5));
				double c = a + ratio * (b - a), d = b - ratio * (b - a);
				double fc = DIST(c), fd = DIST(d);
				rep(ii, 100) {
					if (fc > fd) { // '<': maximum, '>': minimum
						a = c; c = d; d = b - ratio * (b - a);
						fc = fd; fd = DIST(d);
					} else {
						b = d; d = c; c = a + ratio * (b - a);
						fd = fc; fc = DIST(c);
					}
				}
				mids.pb(!dir ? len - c : c);
			}
			}
			sort(all(mids));
			rep(j, mids.size()) {
				P p = s0 + mids[j]*cs;
				pointsx.pb(real(p)), pointsy.pb(imag(p));
			}
		}
		pointsx.pb(roadX[n-1]), pointsy.pb(roadY[n-1]);
		int N = pointsx.size();
		static double dist[250000]; static bool vis[250000];
		fill(dist, dist+N, inf); dist[0] = 0; mset(vis, 0);
		priority_queue<pair<double,int>,vector<pair<double,int> >,greater<pair<double,int> > > q;
		for (q.push(mp(0., 0)); !q.empty(); ) {
			double w = q.top().first; int v = q.top().second; q.pop();
			if (vis[v]) continue; vis[v] = 1;
			rep(u, N) {
				double d = (v+1==u?1:2)*hypot(pointsx[u] - pointsx[v], pointsy[u] - pointsy[v]);
				if (dist[u] > w + d)
					q.push(mp(dist[u] = w + d, u));
				}
		}
		return dist[N-1] / standardDist;
	}
};