Dinkelbach algorithm

D: 食塩水 - AtCoder Beginner Contest 034 | AtCoderのD問題「食塩水」で分数計画問題(分数計画問題 - ORWiki)が出ました.

分数計画問題のアルゴリズムについては次のブログに詳しく載っています.
Dinkelbach algorithm - 週刊 spaghetti_source - TopCoder部
また,上のブログの実装例の最小比全域木問題の詳しい説明は,同じ方の 様々な全域木問題 に載っています.

今回は,「D:食塩水」を上のブログで説明しているDinkelbach algorithmで解きました.コンテストの解説(AtCoder Beginner Contest 034 解説)では,上のブログの二分探索による解法について説明を行っています.これは,蟻本のpp.132-134に掲載されている「平均最大化」と同じです.

[D. 食塩水]

問題

日本語なので省略.

方針

まず与えられた問題を分数計画問題の形にします.
 p_i %の食塩水  w_i gに入っている食塩の量 c_i w_i \times p_i / 100 gです.i番目の容器とj番目の容器の食塩水を混ぜると濃度は, (c_i + c_j) / (w_i + w_j) となります.したがって,容器の集合Xを選んだ時の混ぜた後の濃度は,
   \sum_{i \in X}c_i / \sum_{j \in X}w_j
となります.
ここで,
   f(X) = \sum_{i \in X}c_i  g(X) = \sum_{j \in X}w_j
とすると,求める問題は,
  maximize:   f(X) / g(X)
  subjecto to:  X \subseteq 2^{\{1,2, \cdots, N\}}
という分数計画問題となります.ここで, 2^{\{1, 2, \cdots, N\}} は与えられた食塩水が入った容器の部分集合全体とします.

関数fは食塩量,関数gは食塩水の量への離散関数となっており,最小比全域木問題と同じく{0,1}整数線形計画となっています.したがって,実装は上のブログのDinkelbach algorithmの実装例とほとんど同じになっています.
異なる点は最大化(ソート関数の向きに注意)であることと,h(t)の終了判定の誤差を「1e-8」ではなく「1e-3」にしたところです.「1e-8」にすると時間内に処理が終わりませんでした. 100 \times t は求める濃度を表しているのですが,終了判定の誤差は求める値の誤差と関係がなさそうなのでどれぐらい小さくすれば良いのかが難しいです.

ソースコード
#include <bits/stdc++.h>

using namespace std;

typedef long long  ll;

struct Cost {
    double w, p;
    Cost(double w, double p) : w(w), p(p) {}
    Cost() {}
};
int n, k;
vector<Cost> c;

double Check(double t, double &fcost, double &gcost)
{
    fcost = gcost = 0;

    sort(c.begin(), c.end(), [t](Cost a1, Cost a2) {
            return a1.p - t * a1.w > a2.p - t * a2.w; });

    for (int i = 0; i < k; ++i) {
        fcost += c[i].p;
        gcost += c[i].w;
    }

    return fcost - t * gcost;
}

int main()
{
    cin.tie(0);
    ios::sync_with_stdio(false);
    cout << setprecision(8) << setiosflags(ios::fixed); // printf("%.4f", x);

    cin >> n >> k;
    c.resize(n);

    for (int i = 0; i < n; ++i) {
        cin >> c[i].w >> c[i].p;
        c[i].p *= c[i].w * 0.01; // i番目の食塩の量をp[i]とする
    }

    double t = 0.0, fcost, gcost;
    while (true) {
        if (fabs(Check(t, fcost, gcost)) < 1e-3)
            break ;
        t = fcost / gcost;
    }

    cout << 100.0 * t << endl;

    return 0;
}