にょぐた君と一緒に考察した。面白かった。公式解説は多面体で考えているけど、ここでは LP 双対で解いてみる。

問題概要

$n$ 種類の水溶液がある。水溶液 $i$ は、 $M = 10^{4}$ として、

水溶液全体の重さが $a_{i}$ g
1 g あたりの溶質の重さが $\frac{l_{i}}{M}$ g 以上 $\frac{r_{i}}{M}$ g 以下

であることがわかっている。

これらの水溶液を有理数 g ずつ混ぜていくことで、

重さが $s$ g
1 g あたりの溶質の重さが $\frac{c}{M}$ g

であるようなの水溶液を作りたい。各水溶液の「1 g あたりの溶質の重さ」に幅があるため、作られる水溶液の「1 g あたりの溶質の重さ」にも幅があるが、その幅を加味した上で、 $\frac{c}{M}$ との最悪時の差分が最小となるようなものを求めよ。

最適な水溶液の、水溶液全体の重さと、溶質の重さを出力せよ。

制約

$1 \le n \le 1000$
$1 \le a_{i}, s \le 10^{5}$
$\sum_{i=1}^{n} a_{i} \ge s$
$0 \le c \le M$
$0 \le l_{i} \le r_{i} \le M$

考えたこと

まず、もとの問題を LP で定式化する。混ぜる水溶液 $i$ の重さを $x_{i}$ g、目標との差分を $v$ とすると、次のように定式化できる。

Minimize:

$v$

Subject to:

$0 \le x_{i} \le a_{i}$
$\sum_{i=1}^{n} x_{i} = s$
$\sum_{i=1}^{n} x_{i} l_{i} \ge cs - v$
$\sum_{i=1}^{n} x_{i} r_{i} \le cs + v$

双対をとると、次の LP になる（双対変数を $y_{i}, z, p, q$ とする）。

Maximize:

$-\sum_{i=1}^{n} a_{i}y_{i} + sz + cs(p - q)$

Subject to:

$-y_{i} + z + l_{i} p - r_{i} q \ge 0$
$p + q \le 1$
$y_{i} \ge 0$

まず、 $p + q = 1$ としてよいことがわかる。ある最適解に対して $p, q$ を同じ値だけ増やしても、目的関数の値は変わらず、 $-y_{i} + z + l_{i} p - r_{i} q$ は減少するためだ。そこで、改めて $w = p - q$ とおくと、上の LP は次のように書き直せる。

Maximize:

$-\sum_{i=1}^{n} a_{i}y_{i} + sz + csw$

Subject to:

$-y_{i} + z + \frac{l_{i} + r_{i}}{2}w + \frac{l_{i} - r_{i}}{2} \ge 0$
$y_{i} \ge 0$
$-1 \le w \le 1$

この LP は、 $w$ を固定すると、Greedy に基づく解法で解ける。その部分を考察する。

$w$ を固定したときの最適解

簡単のため、 $b_{i} = \frac{l_{i} + r_{i}}{2}w + \frac{l_{i} - r_{i}}{2}$ とおくと、 $y_{i}$ に関する条件は次のようになる。

$y_{i} \ge z + b_{i}$
$y_{i} \ge 0$

目的関数の式の形から、 $z$ を固定するとき、 $y_{i}$ は極力小さいことが望ましいので、

$y_{i} = \max(z + b_{i}, 0)$

とすればよいことがわかる。つまり、次の関数の最大値を求める問題となる。

$\displaystyle sz - \sum_{i=1}^{n} a_{i} \max(z + b_{i}, 0) + csw$

$\sum_{i=1}^{n} a_{i} \ge s$ であることから、 $z$ はある程度まで小さくした方がよさそうなことがわかる。ただし、 $z$ を小さくしすぎると、 $\max(z + b_{i}, 0) = 0$ となる $i$ が増えていくことに注意しよう。まとめると、次のように整理できる。

$z + b_{i} \gt 0$ であるような $i$ に対する $a_{i}$ の総和が $s$ 以上であるとき、 $z$ は小さくした方がよい
$z + b_{i} \gt 0$ であるような $i$ に対する $a_{i}$ の総和が $s$ 未満であるとき、 $z$ は大きくした方がよい

よって、 $b_{i}$ が小さくなる順に $i$ を並べて、 $i_{1}, i_{2}, \dots, i_{n}$ として、

$a_{i_{1}} + a_{i_{2}} + \dots + a_{i_{K}} \gt \sum_{i=1}^{n} - s$

となる最小の $K$ を見つけて、

$z = -b_{i_{K}}$

とすればよい。以上で、 $w$ を固定したときの LP が解けた！！

元の LP に戻り

我々は、 $w$ を固定したときの、次の LP の最適解を求めることができた。

Maximize:

$-\sum_{i=1}^{n} a_{i}y_{i} + sz + csw$

Subject to:

$-y_{i} + z + \frac{l_{i} + r_{i}}{2}w + \frac{l_{i} - r_{i}}{2} \ge 0$
$y_{i} \ge 0$
$-1 \le w \le 1$

次に、 $w$ を固定したときの、この LP の最適解を $f(w)$ と表して、 $f(w)$ の最適解を求めよう。

関数 $f(w)$ は、高々 $O(n^{2})$ 個の微分不可能点をもつ区分線形関数になる。具体的には、各 $i, j$ に対する、

$\frac{l_{i} + r_{i}}{2}w + \frac{l_{i} - r_{i}}{2} = \frac{l_{j} + r_{j}}{2}w + \frac{l_{j} - r_{j}}{2}$ ...... ①

となるような点である。この左辺と右辺の大小関係が入れ替わるとき、上の議論における $b_{1}, b_{2}, \dots, b_{n}$ の大小関係が入れ替わるためだ。

最初、にょぐた君と一緒に考えたのは、 $O(n^{2})$ 通りの候補を調べればよいのだから、それらの差分を素早く更新できれば、 $O(n^{2} \log n)$ の計算量で解けたのではないか......ということだった。しかしこれは実装が大変だ。

ここで opt さんがスペースにやって来て、ヒントをくれた。関数 $f(w)$ は区分線形凹関数だというのだ。確かに、関数 $f(w)$ の上側の領域は「一次関数の上側の領域を重ね合わせたもの」として表せるのだから、関数 $f(w)$ は凹関数になる！！！

そうであれば、三分探索が使える！！

有理数で答える必要があるから、単純な三分探索ではなく、次のようにする。

① を満たす $w$ を $-1 \le w \le 1$ の範囲で列挙して、 $w = -1, 1$ を含めて、小さい順に $w_{1}, w_{2}, \dots, w_{L}$ とする
新たに自然数関数 $g(x) = f(w_{x})$ を定義すると、これも凹関数である
$1 \le x \le N$ の範囲で三分探索して、 $g(x)$ の最大値を求める

この解法は、候補点のソートがボトルネックであり、計算量は $O(n^{2} \log n)$ となる。

余談

三分探索をするにしても、 $w$ を $-1 \le w \le 1$ を $M$ 等分すれば十分という解法もあるようだ。あと、多面体を用いて幾何的に考察する方法が commentary にある。

コード

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

// 有理数
template<class T> struct frac {
    // inner values
    T first, second;

    // constructor
    frac& normalize() {
        if (second < 0) first = -first, second = -second;
        T d = gcd(abs(first), abs(second));
        if (d == 0) first = 0, second = 1;
        else first /= d, second /= d;
        return *this;
    }
    constexpr frac(T f = 0, T s = 1) : first(f), second(s) { 
        normalize(); 
    }
    constexpr frac& operator = (T a) { 
        *this = frac(a, 1); 
        return *this;
    }

    // comparison operators
    constexpr bool operator == (const frac &r) const {
        return this->first == r.first && this->second == r.second;
    }
    constexpr bool operator != (const frac &r) const {
        return this->first != r.first || this->second != r.second;
    }
    constexpr bool operator < (const frac &r) const {
        return this->first * r.second < this->second * r.first;
    }
    constexpr bool operator > (const frac &r) const {
        return this->first * r.second > this->second * r.first;
    }
    constexpr bool operator <= (const frac &r) const {
        return this->first * r.second <= this->second * r.first;
    }
    constexpr bool operator >= (const frac &r) const {
        return this->first * r.second >= this->second * r.first;
    }
    
    // arithmetic operators
    constexpr frac& operator += (const frac &r) {
        this->first = this->first * r.second + this->second * r.first;
        this->second *= r.second;
        this->normalize();
        return *this;
    }
    constexpr frac& operator -= (const frac &r) {
        this->first = this->first * r.second - this->second * r.first;
        this->second *= r.second;
        this->normalize();
        return *this;
    }
    constexpr frac& operator *= (const frac &r) {
        this->first *= r.first;
        this->second *= r.second;
        this->normalize();
        return *this;
    }
    constexpr frac& operator /= (const frac &r) {
        this->first *= r.second;
        this->second *= r.first;
        this->normalize();
        return *this;
    }
    constexpr frac operator + () const { return frac(*this); }
    constexpr frac operator - () const { return frac(0) - frac(*this); }
    constexpr frac operator + (const frac &r) const { return frac(*this) += r; }
    constexpr frac operator - (const frac &r) const { return frac(*this) -= r; }
    constexpr frac operator * (const frac &r) const { return frac(*this) *= r; }
    constexpr frac operator / (const frac &r) const { return frac(*this) /= r; }
    friend constexpr ostream& operator << (ostream &os, const frac<T> &x) {
        os << x.first; 
        if (x.second != 1) os << "/" << x.second;
        return os;
    }
};


void ICPC_Asia_2024_J() {
    using lfrac = frac<long long>;
    const long long M = 10000;
    long long N, S, C, sumA = 0;
    cin >> N >> S >> C;
    vector<long long> A(N), L(N), R(N);
    for (int i = 0; i < N; i++) cin >> A[i] >> L[i] >> R[i], sumA += A[i];

    vector<lfrac> vw({lfrac(-1), lfrac(1)});
    for (int i = 0; i < N; i++) {
        for (int j = i+1; j < N; j++) {
            lfrac f(-L[i]+R[i]+L[j]-R[j], L[i]+R[i]-L[j]-R[j]);
            if (f >= -1 && f <= 1) vw.push_back(f);
        }
    }
    sort(vw.begin(), vw.end());
    vw.erase(unique(vw.begin(), vw.end()), vw.end());

    auto func = [&](lfrac w) -> lfrac { 
        vector<pair<lfrac,int>> v;
        for (int i = 0; i < N; i++) {
            v.emplace_back(w*lfrac(L[i]+R[i])/2 + lfrac(L[i]-R[i])/2, i);
        }
        sort(v.begin(), v.end());

        vector<lfrac> y(N, -1);
        lfrac z, limit = sumA - S;
        for (auto [val, i] : v) {
            if (limit >= 0) {
                limit -= A[i];
                y[i] = 0;
                if (limit < 0) z = -val;
            } else {
                y[i] = z + val;
            }
        }
        lfrac obj = z * S + w * C * S;
        for (int i = 0; i < N; i++) obj -= y[i] * A[i];
        return obj;
    };

    int left = 0, right = vw.size()-1;
    while (right - left > 2) {
        int m1 = (left * 2 + right) / 3, m2 = (left + right * 2) / 3;
        if (func(vw[m1]) > func(vw[m2])) right = m2;
        else left = m1;
    }
    lfrac res = max({func(vw[left]), func(vw[(left+right)/2]), func(vw[right])});
    res /= M;
    cout << res.first << " " << res.second << endl;
}


int main() {
    ICPC_Asia_2024_J();
}

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

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

AOJ 1464 Mixing Solutions (ICPC アジア 2024 J) (6D)

問題概要

制約

考えたこと

$w$ を固定したときの最適解

元の LP に戻り

余談

コード