有理数さん。。。
誤差がとんでもなく厳しいらしく、有理数ライブラリを使うということをついに覚えた！！！

あと、問題の解法自体は「区間串刺し」といういわゆる「区間スケジューリング問題」の亜種。

問題概要

とある水族館に住むイルカ君は、ジャンプをして $N$ 個のリングをくぐり抜けるとご褒美がもらえます。

イルカ君は座標 $(0, 0)$ から飛び、 $(T, 0)$ で着水する。
ジャンプの軌道は放物線である。
$i$ 番目のリングは、ジャンプの軌道が $(X_i, L_i)$ と $(X_i, H_i)$ を結ぶ線分と交わると、くぐり抜けたと判定される。

1 回のジャンプには初速と同じだけの体力が必要である。

イルカ君は、必要であれば何度でもジャンプをすることができます。重力加速度ベクトルを $(0, −1)$ として、イルカ君が全てのリングを通り抜けるために必要な体力の合計の最小値を求めてください。

絶対誤差または相対誤差が $10^{-9}$ まで許容。

制約

$1 \le N \le 10^{5}$
$1 \le X_i$ < $T \le 10^{6}$
$1 \le L_i$ < $H_i \le 10^{6}$

考えたこと

各リングは、「必要な初速度が $L$ 以上 $R$ 以下」といった制約を表現する区間と考えることができる。

そうするとこの問題は

すべての区間を串刺しにしたい
各串にはコストがある
串のコストの総和を最小にしたい

というタイプの問題になる。もし串の重みが 1 であれば、これは実は区間スケジューリング問題の双対問題になっている (AtCoder ABC 103 D - Islands War (400 点) を参照)。

重みがあると解きづらいが、今回は打ち出し 45 度がコスト最小になっていて、

45 度で倒せず 45 度より小さい角度帯にあるリングはなるべく大きな角度で倒したい
45 度で倒せず 45 度より大きい角度帯にあるリングはなるべく小さな角度で倒したい

という感じになっている。とりあえず例によって区間全体を終端でソートしておいて、

まだ串刺しされていないリングのうち、終端が 45 度より手前にあるなら終端で串刺しして行く (終端でソートしているので、そもそも後ろで刺せば刺すほど他の区間にとっても良い)
終端が 45 度を超えたら急に「終端で刺すのがいい」が崩れてしまい、かといって 45 度で刺すことが残りの区間にとってよいかも定かでない状態になる。
そこで残ったやつを (というかむしろ最初から) 始点が遅い順に区間をソートして、まだ串刺しされていないリングたちを始点で串刺しにして行く
同様に始点が 45 度よりも小さくなったら止める。

これらの操作をやってもし最後にいくつかの区間が残ったら、それらはすべて「45 度を含む」という区間なので、残りまとめて 45 度で刺して終わりになる。

以上をまとめると

区間全体を終端の早い順にソートして、終端が 45 度を超えない範囲を区間スケジューリングする
区間全体を始点の遅い順にソートして、始点が 45 度より小さくならない範囲を区間スケジューリングする
それでもし何も余らないなら終わりで、余ったら最後まとめて 45 度を刺す

という感じになる。注意点として、「45 度を含む区間」があったとしても、いきなりそれを刺してはいけない。上記の 1. 2. の作業ですべての区間を刺せるかもしれないからだ。その場合、45 度を先に刺してしまうと解が悪化してしまう場合もある。

さて、残った巨大な注意点として、この問題は実は誤差が非常にシビアらしい！！！！！！！！！！！！！！！！！！！

今回は放物線だから明らかにそれよりヤバそうで、10^6だと3乗までしか耐えられない(long doubleは精度10^19なので)ことを考えると、やばそうな気がしてくる
— よすぽ (@yosupot) 2019年3月7日

そこで有理数を用いることが推奨される。特に区間の終端を扱う部分がかなりシビアな精度が求められるので。。。

ただ有理数として分母分子を long long 型で表してオーバーフローしない設定にするのも、それはそれで大変だった。。。

#include <iostream>
#include <vector>
#include <cmath>
#include <iomanip>
#include <algorithm>
using namespace std;

// chmax, chmin
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; }

// debug stream of pair, vector 
#define COUT(x) cout << #x << " = " << (x) << " (L" << __LINE__ << ")" << endl
template<class T1, class T2> ostream& operator << (ostream &s, pair<T1,T2> P)
{ return s << '<' << P.first << ", " << P.second << '>'; }
template<class T> ostream& operator << (ostream &s, vector<T> P)
{ for (int i = 0; i < P.size(); ++i) { if (i > 0) { s << " "; } s << P[i]; } return s; }


// 有理数
long long calc_gcd(long long a, long long b) {return b ? calc_gcd(b, a % b) : a;}
struct frac {
    long long first, second;

    using D = long double;
    inline frac normalize() {
        if (second < 0) {first = -first; second = -second;}
        long long d = calc_gcd(first, second);
        if (d == 0) {first = 0; second = 1;}
        else {first /= d; second /= d;}
        return *this;
    }
    frac(long long f = 0, long long s = 1) : first(f), second(s) { normalize(); }
    inline D to_d() const { return D(first) / second; }
    inline frac operator - () { (*this).first *= -1; return (*this); }
    inline const frac& operator = (long long a) { *this = frac(a, 1); return *this; }
    inline const frac& operator += (const frac& a);
    inline const frac& operator += (long long a);
    inline const frac& operator -= (const frac& a);
    inline const frac& operator -= (long long a);
    inline const frac& operator *= (const frac& a);
    inline const frac& operator *= (long long a);
    inline const frac& operator /= (const frac& a);
    inline const frac& operator /= (long long a);
    inline friend ostream& operator << (ostream& s, const frac& f) { 
        s << f.first; if (f.second != 1) s << "/" << f.second; return s;
    }
};
inline bool operator == (const frac &a, const frac&b) {
    return a.first * b.second == a.second * b.first;
}
inline bool operator != (const frac &a, const frac &b) { return !(a == b); }
inline bool operator < (const frac& a, const frac& b) {
    return a.first * b.second < a.second * b.first;
}
inline bool operator > (const frac& a, const frac& b) { return b < a; }
inline bool operator <= (const frac& a, const frac& b) {
    return a.first * b.second <= a.second * b.first;
}
inline bool operator >= (const frac& a, const frac& b) { return b <= a; }
inline frac operator + (const frac& a, const frac& b) {
    frac res;
    res.first = a.first * b.second + a.second * b.first;
    res.second = a.second * b.second;
    res.normalize();
    return res;
}
inline frac operator - (const frac& a, const frac& b) {
    frac res;
    res.first = a.first * b.second - a.second * b.first;
    res.second = a.second * b.second;
    res.normalize();
    return res;
}
inline frac operator * (const frac& a, const frac& b) {
    frac res;
    res.first = a.first * b.first;
    res.second = a.second * b.second;
    res.normalize();
    return res;
}
inline frac operator / (const frac& a, const frac& b) {
    frac res;
    res.first = a.first * b.second;
    res.second = a.second * b.first;
    res.normalize();
    return res;
}
inline frac abs(const frac& a) {
    frac res; res = a; res.normalize(); 
    if (res.first < 0) res.first = res.first * (-1);
    return res;
}
inline const frac& frac::operator += (const frac& x) {*this = *this + x; return *this;}
inline const frac& frac::operator += (long long x) {*this = *this + x; return *this;}
inline const frac& frac::operator -= (const frac& x) {*this = *this - x; return *this;}
inline const frac& frac::operator -= (long long x) {*this = *this + x; return *this;}
inline const frac& frac::operator *= (const frac& x) {*this = *this * x; return *this;}
inline const frac& frac::operator *= (long long x) {*this = *this * x; return *this;}
inline const frac& frac::operator /= (const frac& x) {*this = *this / x; return *this;}
inline const frac& frac::operator /= (long long x) {*this = *this / x; return *this;}



long double cost(const frac &f, long long T) {
    long double df = f.to_d();
    return sqrt(df * T * T / 2 + 0.5 / df);
}

int main() {
    long long T; int N;
    cin >> T >> N;
    const frac center(1, T); // 45 度打ち出しの場合
    using pf = pair<frac,frac>;
    vector<pf> upper, lower, middle;
    for (int i = 0; i < N; ++i) {
        long long x, low, up; cin >> x >> low >> up;
        frac flow(low, x * (T-x));
        frac fup(up, x * (T-x));
        if (fup < center) lower.push_back({flow, fup});
        else if (flow > center) upper.push_back({flow, fup});
        else middle.push_back({flow, fup});
    }
    sort(lower.begin(), lower.end(), [](const pf &a, const pf &b) {
            return a.second < b.second;});
    sort(upper.begin(), upper.end(), [](const pf &a, const pf &b) {
            return a.first > b.first;});

    long double res = 0.0;
    frac left = 0, right = 1000100; // right * 10^12/4 がオーバーフローしないように
    for (auto inter : lower) {
        if (left >= inter.first) continue;
        res += cost(inter.second, T);
        left = inter.second;
    }
    for (auto inter : upper) {
        if (right <= inter.second) continue;
        res += cost(inter.first, T);
        right = inter.first;
    }
    bool remain = false;
    for (auto inter : middle) {
        if (left < inter.first && inter.second < right) remain = true;
    }
    if (remain) res += cost(center, T);
    cout << fixed << setprecision(20) << res << endl;
}

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

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

AOJ 3060 Rings (RUPC 2019 day2-J)

問題概要

制約

考えたこと