F - Sugar Water 2 Editorial by MMNMM

誤差のない値を得る方法

この問題の答えは、既約分数として表したときの分母と分子がたかだか \(2\times10^7\) であるような有理数として表すことができます。 この有理数の分母と分子の値を整数として得る方法について考えます。

Stern-Brocot Tree という木があります。 Stern-Brocot Tree は、すべての正の有理数が \(1\) 度ずつ含まれる無限に大きな二分探索木です。 各ノードには \(4\) つの非負整数 \((a,b,c,d)\) が対応し、そのノードの値は \(\dfrac{a+c}{b+d}\) です。左の子は \((a,b,a+c,b+d)\) 、右の子は \((a+c,b+d,c,d)\) です。
根には \((0,1,1,0)\) が対応します。

「\(K\) 番目の砂糖水の濃度が \(\dfrac rs\) より大きいか?」という判定問題は \(O((N+M)\log(N+M))\) 時間で解くことができるので、これを用いて Stern-Brocot Tree を探索することでこの問題を解くことができます。

しかし、Stern-Brocot Tree における \(\dfrac ab\) の深さは \(\max\lbrace a,b\rbrace\) まで大きくなる場合があるため、単に二分探索木として探索を行うと時間制限に間に合わない場合があります。
そこで、子に潜る際に「同じ向きに何回連続で降りるか」を二分探索することで、\(O(\log(a+b))\) 回判定問題を解くだけでよいことが示せます。

評価の詳細

\(a\gt b\) として一般性を失いません(\(\dfrac ab\) と \(\dfrac ba\) は木の左右対称な位置にあります)。 \(\dfrac ab\) に到達するまでに、右子に \(a _ 1\) 回、左子に \(a _ 2\) 回、⋯、右/左子に \(a _ k\) 回進んだとします。 判定問題を解く回数は \(\displaystyle\sum _ {i=1} ^ k \left(2\log _ 2a _ i+O(1)\right)\) 回です。

\(\dfrac ab=a _ 1+\dfrac 1{a _ 2+\dfrac 1{\ddots+\dfrac 1{a _ k+1}}}\) が成り立っていることに注意します。

連分数の下から、分母 \(+\) 分子 の値を評価することを考えます。

\(S _ 0\coloneqq1,S _ i\coloneqq (a _ {k-i+1}\) 以下の連分数を評価した既約分数の分母と分子の和\(\vphantom{a _ 1})\) と定めます。 すると、

\[S _ {i+1}\geq S _ i\left(1+\dfrac{a _ {k-i}}2\right)\geq S _ i\max\left\lbrace\dfrac32,\dfrac{a _ {k-i}}2\right\rbrace\]

が成り立つことが示せます(\(a _ {k-i}+\dfrac xy\ (x\leq y)\) の分母と分子を計算し、その和を評価すればよいです)。

よって、\(a+b=S _ k\) について \[S _ k\geq\left(\dfrac32\right) ^ k\] および \[S _ k\geq\prod _ {i=1} ^ k\dfrac{a _ i}2\] が成り立ちます。

一つめの式から、\(k=O(\log(a+b))\) が言えます。 二つめの式に \(2 ^ k\) をかけたあと両辺対数をとって、 \[k+\log _ 2(a+b)\geq \sum _ {i=1} ^ k\log _ 2a _ i\] がいえ、\(k=O(\log(a+b))\) より \(\displaystyle\sum _ {i=1} ^ k\log _ 2a _ i=O(\log(a+b))\) がいえました。

よって、判定問題を解く回数 \(\displaystyle\sum _ {i=1} ^ k \left(2\log _ 2a _ i+O(1)\right)\) も \(O(\log(a+b))\) となります。


解 \(\dfrac ab\) に到達しても、それが解であることを知ることは難しい場合があります。 訪れているノードの分母や分子の値が十分大きいとき、その部分木のどのノードの値についても判定関数の結果はすべて等しいです。 ノード \((a,b,c,d)\) の子孫の値の下限は \(\dfrac ab\) で上限は \(\dfrac cd\ \left(\dfrac 10\right.\) は \(+\infty\) とする\(\!\left.\vphantom\dfrac 10\right)\) であることを利用することで、解を求めることができます。

実装例は以下のようになります。 この実装例では、砂糖が \(x\) グラムで水が \(y\) グラムであるような砂糖水の濃度 \(\dfrac {100x}{x+y}\) について探索を行うのではなく、砂糖と水の比 \(\dfrac xy\) によって探索を行なっています。 この \(2\) つのどちらを選択するかで比較結果が変わらないことを示すことができます。 また、この実装例では答えを long double で求めていますが(long double で表すことができる真の値と最も近い値が得られることが示せます)、筆算のような操作を行うことで正確な値を好きな桁数求めることもできます。

#include <utility>

// 準備: pair の和、スカラー倍
template<typename T, typename S>
std::pair<T, S> operator+(const std::pair<T, S> &lhs, const std::pair<T, S> &rhs) {
    return {lhs.first + rhs.first, lhs.second + rhs.second};
}

template<typename T>
std::pair<T, T> operator*(const T &lhs, const std::pair<T, T> &rhs) {
    return {lhs * rhs.first, lhs * rhs.second};
}

// Stern-Brocot Tree の探索
//   judge_function(a, b) は非負有理数から {true, false} への写像であって、
//   ある分子・分母ともに max_value 以下であるような有理数 x/y によって a/b と x/y の大小関係で true, false が定まるもの。
//   x/y を O(log max_value) 時間で求める。
template<typename Integer, typename Judge>
std::pair<Integer, Integer> stern_brocot_tree_search(Judge &&judge_function, const Integer &max_value) {
    std::pair<Integer, Integer> lower{0, 1}, upper{1, 0}, now{1, 1}; // 下界 0/1 上界 1/0
    while (true) {
        now = lower + upper;
        const auto now_judge{judge_function(now)};
        auto &from{now_judge ? lower : upper}, &to{now_judge ? upper : lower}; // from から to へ向かって潜っていく

        // 指数探索
        // 1. 上限を探索
        unsigned long L{1}, R{2};
        while (judge_function(from + R * to) == now_judge) {
            L *= 2;
            R *= 2;
            // max_value より下まで潜ったら、それ以降無限に潜る -> 答えは to (= from + ∞ * to)
            if ((from + L * to).first > max_value || (from + L * to).second > max_value)return to;
        }
        // 2. 二分探索
        while (L + 1 < R) {
            const auto M{(L + R) / 2};
            (judge_function(from + M * to) == now_judge ? L : R) = M;
        }

        from = from + L * to;
    }
}

#include <iostream>
#include <vector>
#include <algorithm>
#include <numeric>
#include <iomanip>

int main() {
    using namespace std;

    unsigned long N, M, K;
    cin >> N >> M >> K;

    vector<pair<unsigned long, unsigned long>> takahashi_sugar_water(N), aoki_sugar_water(M);
    for (auto&&[A, B] : takahashi_sugar_water)
        cin >> A >> B;
    for (auto&&[C, D] : aoki_sugar_water)
        cin >> C >> D;

    constexpr unsigned long INF{40000000000};

    // check(a, b) := 混ぜた砂糖水のうち、砂糖/水 <= a/b であるようなものが K 個以下か
    auto check{[N, M, K,
                &takahashi_sugar_water,
                &aoki_sugar_water,
                value{vector<unsigned long>(N + M)},
                index{[N, M] {
                    vector<unsigned long> ret(N + M);
                    iota(begin(ret), end(ret), 0UL);
                    return ret;
                }()}
               ](auto &&fraction) mutable {
        const auto&[x, y]{fraction};
        for (unsigned long i{}; i < N; ++i) {
            const auto&[A, B]{takahashi_sugar_water[i]};
            value[i] = B * x - A * y + INF;
        }
        for (unsigned long i{}; i < M; ++i) {
            const auto&[C, D]{aoki_sugar_water[i]};
            value[i + N] = C * y - D * x + INF;
        }
        sort(begin(index), end(index), [&value](auto i, auto j) { return make_pair(value[i], i) > make_pair(value[j], j); });

        unsigned long takahashi_count{}, pair_count{};
        for (const auto i : index)(i < N ? takahashi_count : pair_count) += (i < N ? 1 : takahashi_count);
        return pair_count < N * M - K + 1;
    }};

    // 分母・分子はたかだか max(max(A) + max(C), max(B) + max(D))
    const auto max_value{[chmax{[](auto &&x, const auto &y) { if (x < y)x = y; }}](auto &&a, auto &&b) {
        unsigned long ans_x{}, ans_y{}, ans_z{}, ans_w{};
        for (const auto&[x, y] : a) {
            chmax(ans_x, x);
            chmax(ans_y, y);
        }
        for (const auto&[z, w] : b) {
            chmax(ans_z, z);
            chmax(ans_w, w);
        }
        return max(ans_x + ans_z, ans_y + ans_w);
    }(takahashi_sugar_water, aoki_sugar_water)};

    const auto[x, y]{stern_brocot_tree_search(check, max_value)};
    // 答えは 100 x / (x + y)
    cout << fixed << setprecision(20) << static_cast<long double>(100 * x) / (x + y) << endl;

    return 0;
}

posted:
last update: