AOJ 1183 鎖中経路 - mayoko’s diary

幾何むずい。

問題

解法

問題文中でご丁寧に示してあるように,最適解では円の交点を通りながらゴールに向かうことになります。

円の交点を求める公式はこちらを参考にしました。
円と円の交点を求める - Shogo Computing Laboratory

ということで,
dp[i][k] := (i番目の円とi+1番目の円のk番目の交点までの最短距離)
というdpを考えることが出来ます(kは0か1。それぞれの円の交点は2つ存在するので)。最初の円と最後の円は交点ではなく中心のみを考えるので適当に処理します。

少し難しいのは「どういう時にi番目の円からk番目の円まで一直線に移動できるか」を判定することです。今回はj = i+1〜 k-1についてj番目とj+1番目の円の交点から得られる線分を考えている直線が通過するかどうかでこれを判定しています。

double eps = 1e-9;

double add(double a, double b) {
    if (abs(a+b) < eps * (abs(a)+abs(b))) return 0;
    return a+b;
}

bool equal(double a, double b) {
    return add(a, -b) == 0;
}

struct P {
    double x, y;
    P() {}
    P(double x, double y) : x(x), y(y) {}
    P operator+(P p) const {return P(add(x, p.x), add(y, p.y));}
    P operator-(P p) const {return P(add(x, -p.x), add(y, -p.y));}
    P operator*(double d) const {return P(x*d, y*d);}
    double dot(P p) const {return add(x*p.x, y*p.y);} // 内積
    double det(P p) const {return add(x*p.y, -y*p.x);} // 外積
    double dist(P p) const {return sqrt((x-p.x)*(x-p.x) + (y-p.y)*(y-p.y));} // 距離
    bool operator<(const P& rhs) const {
        if (x != rhs.x) return x < rhs.x;
        return y < rhs.y;
    }
    bool operator==(const P& rhs) const {
        return equal(x, rhs.x) && equal(y, rhs.y);
    }
};

// 線分p1-p2上に点qがあるかを判定する
bool on_seg(P p1, P p2, P q) {
    return (p1-q).det(p2-q) == 0 && (p1-q).dot(p2-q) <= 0;
}

// 直線p1-p2と直線q1-q2の交点
P intersection(P p1, P p2, P q1, P q2) {
    return p1+(p2-p1)*((q2-q1).det(q1-p1)/(q2-q1).det(p2-p1));
}

double square(double x) {return x*x;}

const int MAXN = 101;
const double INF = 1e18;
int n;
double x[MAXN], y[MAXN], r[MAXN];
double dp[MAXN][2];
P ps[MAXN][2];

bool ok(int from, int fromi, int to, int toi) {
    for (int i = from+1; i < to; i++) {
        P inter = intersection(ps[from][fromi], ps[to][toi], ps[i][0], ps[i][1]);
        if (on_seg(ps[i][0], ps[i][1], inter)) continue;
        else return false;
    }
    return true;
}

void solve() {
    ps[0][0].x = ps[0][1].x = x[0];
    ps[0][0].y = ps[0][1].y = y[0];
    ps[n][0].x = ps[n][1].x = x[n-1];
    ps[n][0].y = ps[n][1].y = y[n-1];
    for (int i = 0; i <= n; i++) for (int j = 0; j < 2; j++) {
        dp[i][j] = INF;
    }
    dp[0][0] = dp[0][1] = 0;
    for (int i = 0; i < n-1; i++) {
        double x1 = x[i+1] - x[i];
        double y1 = y[i+1] - y[i];
        double r1 = r[i], r2 = r[i+1];
        double a = (square(x1) + square(y1) + square(r1) - square(r2)) / 2;
        ps[i+1][0].x = (a*x1 + y1 * sqrt((square(x1) + square(y1)) * square(r1) - square(a))) / (square(x1) + square(y1)) + x[i];
        ps[i+1][1].x = (a*x1 - y1 * sqrt((square(x1) + square(y1)) * square(r1) - square(a))) / (square(x1) + square(y1)) + x[i];
        ps[i+1][0].y = (a*y1 - x1 * sqrt((square(x1) + square(y1)) * square(r1) - square(a))) / (square(x1) + square(y1)) + y[i];
        ps[i+1][1].y = (a*y1 + x1 * sqrt((square(x1) + square(y1)) * square(r1) - square(a))) / (square(x1) + square(y1)) + y[i];
    }
    //for (int i = 0; i <= n; i++) {
    //    cout << i << endl;
    //    cout << ps[i][0].x << "  " << ps[i][0].y << endl;
    //    cout << ps[i][1].x << "  " << ps[i][1].y << endl;
    //}
    for (int i = 0; i < n; i++) for (int j = 0; j < 2; j++) {
        for (int k = i+1; k <= n; k++) for (int l = 0; l < 2; l++) {
            if (ok(i, j, k, l)) {
                dp[k][l] = min(dp[k][l], dp[i][j] + ps[i][j].dist(ps[k][l]));
            }
        }
    }
    printf("%.15lf\n", dp[n][0]);
}

int main() {
    cin.tie(0);
    ios::sync_with_stdio(false);
    while (cin >> n) {
        if (n == 0) break;
        for (int i = 0; i < n; i++) {
            cin >> x[i] >> y[i] >> r[i];
        }
        solve();
    }
    return 0;
}