この見た目で「幾何」になるの、面白い！

問題へのリンク

問題概要

高橋君は $N$ 種類の泳ぎ方ができる。 $i$ 種類目の泳ぎ方では、1 秒間に体力を $A_{i}$ だけ消費して、 $B_{i}$ メートル進むことができる。このとき、次の $Q$ 回のクエリに答えよ。

【クエリ】
正の実数 $C, D$ が与えられる。泳ぎ方を組み合わせることで、次の条件を満たせるかを判定し、満たせるならば満たすまでの最小時間を求めよ。

体力消費量が合計 $C$ 以下である
$D$ メートル以上進む

制約

$1 \le N, Q \le 2 \times 10^{5}$

考えたこと

まず、1 回のクエリについて考えることとする。

$i$ 種類目の泳ぎ方をした時間を $t_{i}$ 秒とすると、次の LP を解けば良いことになる。

Minimize:
$t_{0} + t_{1} + \dots + t_{N-1}$

s.t.

$A_{0}t_{0} + A_{1}t_{1} + \dots + A_{N-1}t_{N-1} \le C$
$B_{0}t_{0} + B_{1}t_{1} + \dots + B_{N-1}t_{N-1} \ge D$
$t_{0}, t_{1}, \dots, t_{N-1} \ge 0$

いかにも多角形を考えたくなるような問題設定である。多角形を考えるような問題にするために、

$\displaystyle T = t_{0} + t_{1} + \dots + t_{N-1}$
$x = \frac{A_{0}t_{0} + A_{1}t_{1} + \dots + A_{N-1}t_{N-1}}{T}$
$y = \frac{B_{0}t_{0} + B_{1}t_{1} + \dots + B_{N-1}t_{N-1}}{T}$

とおく。すると、 $N$ 点 $(A_{0}, B_{0}), (A_{1}, B_{1}, \dots, (A_{N-1}, B_{N-1})$ の凸包を $S$ として、次の 2 つは等価となる。

$t_{0}, t_{1}, \dots, t_{N-1}$ が正の実数全体を動く
二次元座標平面上で点 $(x, y)$ は $S$ の周上または内部の点である

ここで、さらに $U = \frac{1}{T}$ とおくと、上の最適化問題は、次の最適化問題と等価となることがわかる。

Maximize:
$U$

s.t.

$x \le CU$
$y \ge DU$
$(x, y)$ は凸多角形 $S$ の周上または内部の点

この答えは、以下の条件を満たすような点のうちの Y 座標の最大値となる。

凸多角形 $S$ の周上または内部の点
直線 $y = \frac{D}{C}x$ 上、または、直線 $y = \frac{D}{C}x$ よりも上側にある点

残りは、そのような点を高速に見つける方法を考える。

クエリ対応

結局、次の 2 つを試せばよいこととなる。これらのうち、Y 座標が大きい方 (で $D$ を割った値) を答えればよい。

凸多角形 $S$ の頂点であって、Y 座標が最大の点 (ただし、直線 $y = \frac{D}{C}x$ よりも下側にある場合は除外)
凸多角形 $S$ の周囲と、直線 $y = \frac{D}{C}x$ との交点

後者を求めるために、二分探索による解法を考える。まず、下図のように、 $S$ の周のうち "右上側" にあるもののみを抽出した。

そして、この折れ線をなす線分のうち、直線 $y = \frac{D}{C}x$ と交わるものを二分探索によって特定し、その交点を求めた。

以上により、全体の計算量は $O((N + Q) \log N)$ となった。

メモ

一般に、凸多角形と直線の交点を $O( \log N)$ の計算量で求める方法もあるらしい。

tjkendev.github.io

コード

#include <bits/stdc++.h>
using namespace std;

using pint = pair<int, int>;
using pll = pair<long long, long long>;
template<class T> inline bool chmax(T& a, T b) { if (a < b) { a = b; return 1; } return 0; }
template<class T> inline bool chmin(T& a, T b) { if (a > b) { a = b; return 1; } return 0; }


//------------------------------//
// 幾何の基本要素 (点, 線分, 円)
//------------------------------//

// basic settings
long double EPS = 1e-10;  // to be set appropriately
constexpr long double PI = 3.141592653589793238462643383279502884L;
long double torad(long double deg) {return (long double)(deg) * PI / 180;}
long double todeg(long double ang) {return ang * 180 / PI;}

// Point or Vector
template<class DD> struct Point {
    // inner value
    DD x, y;
    
    // constructor
    constexpr Point() : x(0), y(0) {}
    constexpr Point(DD x, DD y) : x(x), y(y) {}
    
    // various functions
    constexpr Point conj() const {return Point(x, -y);}
    constexpr DD dot(const Point &r) const {return x * r.x + y * r.y;}
    constexpr DD cross(const Point &r) const {return x * r.y - y * r.x;}
    constexpr DD norm() const {return dot(*this);}
    constexpr long double abs() const {return sqrt(norm());}
    constexpr long double arg() const {
        if (this->eq(Point(0, 0))) return 0L;
        else if (x < -EPS && this->eq(Point(x, 0))) return PI;
        else return atan2((long double)(y), (long double)(x));
    }
    constexpr bool eq(const Point &r) const {return (*this - r).abs() <= EPS;}
    constexpr Point rot90() const {return Point(-y, x);}
    constexpr Point rot(long double ang) const {
        return Point(cos(ang) * x - sin(ang) * y, sin(ang) * x + cos(ang) * y);
    }
    
    // arithmetic operators
    constexpr Point operator - () const {return Point(-x, -y);}
    constexpr Point operator + (const Point &r) const {return Point(*this) += r;}
    constexpr Point operator - (const Point &r) const {return Point(*this) -= r;}
    constexpr Point operator * (const Point &r) const {return Point(*this) *= r;}
    constexpr Point operator / (const Point &r) const {return Point(*this) /= r;}
    constexpr Point operator * (DD r) const {return Point(*this) *= r;}
    constexpr Point operator / (DD r) const {return Point(*this) /= r;}
    constexpr Point& operator += (const Point &r) {
        x += r.x, y += r.y;
        return *this;
    }
    constexpr Point& operator -= (const Point &r) {
        x -= r.x, y -= r.y;
        return *this;
    }
    constexpr Point& operator *= (const Point &r) {
        DD tx = x, ty = y;
        x = tx * r.x - ty * r.y;
        y = tx * r.y + ty * r.x;
        return *this;
    }
    constexpr Point& operator /= (const Point &r) {
        return *this *= r.conj() / r.norm();
    }
    constexpr Point& operator *= (DD r) {
        x *= r, y *= r;
        return *this;
    }
    constexpr Point& operator /= (DD r) {
        x /= r, y /= r;
        return *this;
    }

    // friend functions
    friend ostream& operator << (ostream &s, const Point &p) {
        return s << '(' << p.x << ", " << p.y << ')';
    }
    friend constexpr Point conj(const Point &p) {return p.conj();}
    friend constexpr DD dot(const Point &p, const Point &q) {return p.dot(q);}
    friend constexpr DD cross(const Point &p, const Point &q) {return p.cross(q);}
    friend constexpr DD norm(const Point &p) {return p.norm();}
    friend constexpr long double abs(const Point &p) {return p.abs();}
    friend constexpr long double arg(const Point &p) {return p.arg();}
    friend constexpr bool eq(const Point &p, const Point &q) {return p.eq(q);}
    friend constexpr Point rot90(const Point &p) {return p.rot90();}
    friend constexpr Point rot(const Point &p, long long ang) {return p.rot(ang);}
};

// necessary for some functions
template<class DD> constexpr bool operator < (const Point<DD> &p, const Point<DD> &q) {
    return (abs(p.x - q.x) > EPS ? p.x < q.x : p.y < q.y);
}

// Line
template<class DD> struct Line : vector<Point<DD>> {
    Line(Point<DD> a = Point<DD>(0, 0), Point<DD> b = Point<DD>(0, 0)) {
        this->push_back(a);
        this->push_back(b);
    }
    friend ostream& operator << (ostream &s, const Line<DD> &l) {
        return s << '{' << l[0] << ", " << l[1] << '}';
    }
};

// 交点
template<class DD> Point<DD> proj_for_crosspoint(const Point<DD> &p, const Line<DD> &l) {
    DD t = dot(p - l[0], l[1] - l[0]) / norm(l[1] - l[0]);
    return l[0] + (l[1] - l[0]) * t;
}
template<class DD> vector<Point<DD>> crosspoint(const Line<DD> &l, const Line<DD> &m) {
    vector<Point<DD>> res;
    DD d = cross(m[1] - m[0], l[1] - l[0]);
    if (abs(d) < EPS) return vector<Point<DD>>();
    res.push_back(l[0] + (l[1] - l[0]) * cross(m[1] - m[0], m[1] - l[0]) / d);
    return res;
}

// 凸包 (一直線上の3点を含めない) の上側だけとってくる
template<class DD> vector<Point<DD>> convex_hull(vector<Point<DD>> &ps) {
    int n = (int)ps.size();
    if (n == 1) return ps;
    vector<Point<DD>> res(2*n);
    auto cmp = [&](const Point<DD> &p, const Point<DD> &q) -> bool {
        return (abs(p.x - q.x) > EPS ? p.x < q.x : p.y < q.y);
    };
    sort(ps.begin(), ps.end(), cmp);
    int k = 0;
    for (int i = 0; i < n; ++i) {
        if (k >= 2) {
            while (cross(res[k-1] - res[k-2], ps[i] - res[k-2]) < EPS) {
                --k;
                if (k < 2) break;
            }
        }
        res[k] = ps[i]; ++k;
    }
    int t = k+1;
    for (int i = n-2; i >= 0; --i) {
        if (k >= t) {
            while (cross(res[k-1] - res[k-2], ps[i] - res[k-2]) < EPS) {
                --k;
                if (k < t) break;
            }
        }
        res[k] = ps[i]; ++k;
    }
    res.resize(k-1);
    return res;
}

using DD = long double;

int main() {
    int N, Q;
    cin >> N;
    vector<Point<DD>> ab(N);
    for (int i = 0; i < N; ++i) cin >> ab[i].x >> ab[i].y;
    
    // 凸包を求めて、原点から見て外側の折れ線カーブのみを取り出す
    auto ch = convex_hull(ab);
    DD min_arg = 1LL<<60;
    int min_p = -1;
    for (int i = 0; i < ch.size(); ++i) {
        if (chmin(min_arg, ch[i].y / ch[i].x)) min_p = i;
    }
    vector<Point<DD>> ps({ch[min_p]});
    for (int i = min_p + 1; i < ch.size() * 2; ++i) {
        if (ch[i % ch.size()].y / ch[i % ch.size()].x
            < ch[(i-1) % ch.size()].y / ch[(i-1) % ch.size()].x) {
            break;
        }
        ps.push_back(ch[i % ch.size()]);
    }
    
    // 凸包全体で y 座標が最大の点を求める
    DD max_y = 0;
    Point<DD> max_point;
    for (int i = 0; i < ch.size(); ++i) {
        if (chmax(max_y, ch[i].y)) max_point = ch[i];
    }

    // クエリ処理
    cin >> Q;
    while (Q--) {
        DD C, D;
        cin >> C >> D;
        Line<DD> l(Point<DD>(0, 0), Point<DD>(C, D));
        
        if (D / C > ps.back().y / ps.back().x + EPS) {
            cout << -1 << endl;
            continue;
        }
        if (max_point.y / max_point.x >= D / C - EPS) {
            cout << fixed << setprecision(10) << D / max_y << endl;
            continue;
        }
        
        int left = 0, right = ps.size();
        while (right - left > 1) {
            int mid = (left + right) / 2;
            if (D / C > ps[mid].y / ps[mid].x) left = mid;
            else right = mid;
        }
        Line l2(ps[left], ps[left+1]);
        
        auto crs = crosspoint(l, l2);
        DD res = 1;
        if (max_point.y / max_point.x > D / C) chmax(res, max_y);
        if (!crs.empty()) chmax(res, crs[0].y);

        cout << fixed << setprecision(10) << D / res << endl;
    }
}

けんちょんの競プロ精進記録

競プロの精進記録や小ネタを書いていきます

AtCoder ABC 356 G - Freestyle (5D, 赤色, 575 点)

問題概要

制約

考えたこと

クエリ対応

メモ

コード