The 46th World Championship ICPC T問題：Carl’s Vacation

問題. Carl’s Vacation

2 つのピラミッドがある。そのピラミッドの頂上間の最短距離を求めよ。

ピラミッドの形状は $(x_1, y_1, x_2, y_2, h)$ で表される。ピラミッドの形状は真上から見ると正方形で、正方形のある一辺をピラミッドを左にしながら点 $(x_1, y_1)$ から点 $(x_2, y_2)$ へ向かう線分で表す。ピラミッドの頂点は正方形の中心で高さ $h$ である。
2 つのピラミッドは交差しない。

制約： $-10^5 \le x_1, x_2, y_1, y_2 \le 10^5$ $\left( (x_1, y_1) \neq (x_2, y_2) \right)$ , $1 \le h \le 10^5$

解法．三分探索

初めに、ピラミッドの頂上から地上との境界である辺上の一点 $P$ への最短距離を考える。
ピラミッドを真上から見た図を下に載せる。対称性から青色の部分についてのみ考える。 $P$ の座標はパラメータ $t \in [0, S / 2]$ を用いて $P = P_1 + t \times (P_2 - P_1)$ となる。このときの距離は $d(t) = \sqrt{t^2 - S t + \frac{S^2}{2} + h^2}$ である。ただし、 $S$ は正方形の一辺の長さとする。 $d(t)$ は $t$ に関する凸関数である。

それぞれのピラミッドの頂上から辺上までの距離はそれぞれパラメータ $t_1 \in [0, S_1 / 2 ]$ 、 $t_2 \in [0, S_2 / 2]$ として $d(t_1)$ と $d(t_2)$ となる（ $S_1, S_2$ はそれぞれ一番目と二番目のピラミッドの一辺の長さ）。また、地上でのピラミッド間の距離はユークリッド距離で、 $t_1, t_2$ の凸関数として表される。したがって、2 つのピラミッドの頂点間の距離は凸関数の和となり、全体で凸関数なので三分探索で答えを求めることができる（三分探索 - 忘れても大丈夫）。

ソースコードを表示

#include <iostream>
#include <iomanip>
#include <cmath>
#include <limits>

using ld = long double;

struct Point {
    ld x, y;

    Point() : x(0.0), y(0.0) {}
    Point(ld x, ld y) : x(x), y(y) {}

    // 反時計回りに90度回転
    void rotate90() {
        *this = Point(-y, x);
    }
};

ld distance(const Point &p1, const Point &p2) {
    return std::hypot(p1.x - p2.x, p1.y - p2.y);
}

Point advance(const Point &p1, const Point &p2, ld t) {
    Point dir(p2.x - p1.x, p2.y - p1.y);
    ld len = std::hypot(dir.x, dir.y);
    dir.x /= len;
    dir.y /= len;
    return Point(p1.x + t * dir.x, p1.y + t * dir.y);    
}

struct Pyramid {
    Point e[2];
    ld S, H;

    void next() {
        Point v(e[1].x - e[0].x, e[1].y - e[0].y);
        v.rotate90();
        e[0] = e[1];
        e[1] = Point(e[1].x + v.x, e[1].y + v.y);
    }

    void swap() {
        std::swap(e[0], e[1]);
    }
};

struct Solver {
    Solver() {}

    Pyramid org_pm1, org_pm2;

    void Input() {
        std::cin >> org_pm1.e[0].x >> org_pm1.e[0].y >> org_pm1.e[1].x >> org_pm1.e[1].y >> org_pm1.H;
        std::cin >> org_pm2.e[0].x >> org_pm2.e[0].y >> org_pm2.e[1].x >> org_pm2.e[1].y >> org_pm2.H;
        org_pm1.S = distance(org_pm1.e[0], org_pm1.e[1]);
        org_pm2.S = distance(org_pm2.e[0], org_pm2.e[1]);
    }

    ld MinimumDistance() {
        ld min_d = std::numeric_limits<ld>::max();
        Pyramid pm1 = org_pm1;

        for (int i = 0; i < 4; ++i) {
            for (int i1 = 0; i1 < 2; ++i1) {
                pm1.swap();
                Pyramid pm2 = org_pm2;

                for (int j = 0; j < 4; ++j) {
                    for (int j1 = 0; j1 < 2; ++j1) {
                        pm2.swap();
                        min_d = std::min(min_d, MinimumDistance(pm1, pm2));
                    }
                    pm2.next();
                }
            }
            pm1.next();
        }

        return min_d;
    }

    ld MinimumDistance(Pyramid pm1, Pyramid pm2) {
        ld lb = 0.0, ub = 0.5 * pm1.S;
        for (int i = 0; i < 148; ++i) {
            ld t1 = (2.0 * lb + ub) / 3.0;
            ld t2 = (lb + 2.0 * ub) / 3.0;

            if (f(t1, pm1, pm2) > f(t2, pm1, pm2)) lb = t1;
            else ub = t2;
        }

        return f(lb, pm1, pm2);
    }

    ld f(const ld t, const Pyramid &pm1, const Pyramid &pm2) {
        const Point p1 = advance(pm1.e[0], pm1.e[1], t);
        const ld dist = std::sqrt(t * t - pm1.S * t + 0.5 * pm1.S * pm1.S + pm1.H * pm1.H);

        ld lb = 0.0, ub = 0.5 * pm2.S;
        for (int i = 0; i < 148; ++i) {
            ld t1 = (2.0 * lb + ub) / 3.0;
            ld t2 = (lb + 2.0 * ub) / 3.0;

            if (g(t1, p1, pm2) > g(t2, p1, pm2)) lb = t1;
            else ub = t2;
        }

        return g(lb, p1, pm2) + dist;
    }

    ld g(const ld t, const Point& p1, const Pyramid &pm2) {
        const Point p2 = advance(pm2.e[0], pm2.e[1], t);
        const ld dist1 = std::sqrt(t * t - pm2.S * t + 0.5 * pm2.S * pm2.S + pm2.H * pm2.H);
        const ld dist2 = distance(p1, p2);
        return dist1 + dist2;
    }
};

int main() {
    std::cout << std::fixed << std::setprecision(10);

    Solver solver;
    solver.Input();
    std::cout << solver.MinimumDistance() << std::endl;

    return 0;
}