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

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

Codeforces #189 (Div. 1) C. Kalila and Dimna in the Logging Industry (R2400)

sky さんの Monotone Minimma の例題!!!
練習として解いてみた。

問題へのリンク

問題概要 (意訳)

 N 個の値の組  a_{0}, \dots, a_{N-1},  b_{0}, \dots, b_{N-1} が与えられる。 a_{0} = 1 b_{N-1} = 0 であり、 a は狭義単調増加、 b は狭義単調減少である。

 0 = k_{0} \lt k_{1} \lt \dots \lt k_{K} = N-1 を適切に定めたときのスコアが、

  •  \sum_{i = 0}^{K-1} a_{k_{i+1}} \times b_{k_{i}}

で与えられる。スコアの最小値を求めよ。

制約

  •  1 \le N \le 10^{5}

あらすじ

文字通り、以下のような DP が立つ。

  • dp[  i ] :=  i 番目の要素までについてのスコア

として、

  • dp[  i ] =  \min_{0 \le j \lt i} (dp[  j ] +  a_{i} \times b_{j}

このままだと  O(N^{2}) であるが、色々な高速化手法が考えられる。

解法 1: Convex Hull Trick

Convex Hull Trick

  • insert (a, b): 直線  y = ax + b を追加する
  • query (x):  {\rm min}_{i} (a_{i}*x + b_{i}) を答える

を高速にできるデータ構造で、一般にともに  O(\log{n}) でできる。しかし嬉しい仮定として

  • insert される順番に  a が単調減少 (直線の傾きが単調減少)
  • query される順番に  x が単調増加 (クエリが単調増加)

の場合はならし  O(1) でできる。今回の DP は、

  • dp[  i ] の値を更新したら、直線  y = b_{i} + dp[  i ] を追加する
  • dp[  i ] の値を求めるためには、 x = a_{i} として、 \min_{0 \le j \lt i} (b_{j} x + dp[  j )] を答える

という風にすればよい。まさに直線の傾きもクエリも単調なので、 O(1) な更新を行うことができる。計算量は  O(N)

なお、この問題は普通に CHT をやると long long 型でもオーバーフローしかねないことに注意。

  • 係数  a 10^{9} オーダー
  • dp の値は  10^{14} くらいになりうる
  • ということで、不用意に  ax + b を挿入すると、 b の値が  10^{14} とかのため、線分の傾き比較とかで分母を払ったりしたときに long long 型のままやるとオーバーフローする (double 型なら通った)
#include <iostream>
#include <vector>
#include <deque>
#include <algorithm>
using namespace std;

// Convex Hull Trick
/*
    追加クエリの直線の傾きが単調減少の場合
 
    - insert (a, b): add y = ax + b
    - query (x): min_i{a[i]*x + b}
*/
template<class T> struct CHT {
    using Line = pair<T, T>;
    deque<Line> lines;
    
    inline bool check(const Line &a, const Line &b, const Line &c) {
        return (double)(b.first - a.first) * (c.second - b.second)
            >= (double)(b.second - a.second) * (c.first - b.first);
    }
    
    inline double get(int k, T x) {
        return (double)lines[k].first * x + lines[k].second;
    }
    
    void insert(T a, T b) {
        Line l(a, b);
        while (lines.size() >= 2
               && check(lines[(int)lines.size()-2], lines[(int)lines.size()-1], l))
            lines.pop_back();
        lines.push_back(l);
    }
    
    T query(T x) {
        int low = -1, high = (int)lines.size();
        while (high - low > 1) {
            int mid = (low + high) / 2;
            if (get(mid, x) >= get(mid + 1, x)) low = mid;
            else high = mid;
        }
        return get(high, x);
    }
    
    // クエリの単調性も成り立つ場合 (x が単調増加)
    T query_monotone(T x) {
        while (lines.size() >= 2 && get(0, x) >= get(1, x)) lines.pop_front();
        return lines[0].first * x + lines[0].second;
    }
};

int main() {
    int N;
    cin >> N;
    vector<long long> a(N), b(N);
    for (int i = 0; i < N; ++i) cin >> a[i];
    for (int i = 0; i < N; ++i) cin >> b[i];

    CHT<long long> cht;
    vector<long long> dp(N+1, 1LL<<61);
    dp[0] = 0;
    cht.insert(b[0], dp[0]);
    for (int i = 1; i < N; ++i) {
        dp[i] = cht.query_monotone(a[i]);
        cht.insert(b[i], dp[i]);
    }
    cout << dp[N-1] << endl;
}

解法 2: Monotone Minima を利用した分割統治法

大前提として、


二変数関数  f(i, j) が monotone ( g(j) = {\rm argmin}_{j} f(i, j) が広義単調増加) であるとき、分割統治法を用いて  g(j) O(N\log{N}) で求められる


これを活用して、

  • dp[  i ] =  \min_{0 \le j \lt i} (dp[  j ] +  w(i, j))

という形をした更新が  O(N(\log{N})^{2}) でできることがある。具体的には再び分割統治法をすることにする。

  • dp[ left : mid ] の区間で完結する遷移の分は先に緩和してしまう
  • dp[ left : mid ] の区間から dp[ mid : right ] の区間への遷移部分に Monotone Minima を活用する
  • dp[ mid : right ] の区間で完結する遷移の分を緩和する

という二重の分割統治法になる。ここで、真ん中の Monotone Minima は、具体的には

  •  i = {\rm mid}, {\rm mid} + 1, \dots, {\rm right} - 1 に対して、 f(i, j) = dp[  j ] +  w(i, j) ( j = {\rm left}, {\rm left} + 1, \dots, {\rm mid - 1}) で定まる  f(i, j) が monotone であるとき、Monotone Minima を行う

ということになる。ここで  f が monotone であるための重要な十分条件として、遷移が交差しない (ここ参照) というのがあるようだ。

今回は  a が単調増加であることから直接言える。

#include <iostream>
#include <vector>
using namespace std;
template<class T> inline bool chmin(T& a, T b) { if (a > b) { a = b; return 1; } return 0; }

// 入力
int N;
vector<long long> a, b;

// DP
vector<long long> dp;

// monotone minima function
long long calc(int i, int j) {
    return dp[j] + a[i] * b[j];
}

// i:[mid:rght), j:[left:mid)
void rec2(int ileft, int iright, int jleft, int jright) {
    if (iright - ileft < 1) return;
    int imid = (ileft + iright) / 2;
    long long mi = calc(imid, jleft);
    int pj = jleft;
    for (int j = jleft; j < jright; ++j) if (chmin(mi, calc(imid, j))) pj = j;
    chmin(dp[imid], mi);
    rec2(ileft, imid, jleft, pj+1);
    rec2(imid+1, iright, pj, jright);
}

void rec1(int left, int right) {
    if (right - left < 2) return;
    int mid = (left + right) / 2;
    rec1(left, mid);
    rec2(mid, right, left, mid);
    rec1(mid, right);
}

int main() {
    cin >> N;
    a.resize(N); b.resize(N);
    for (int i = 0; i < N; ++i) cin >> a[i];
    for (int i = 0; i < N; ++i) cin >> b[i];
    dp.assign(N+1, 1LL<<60);
    dp[0] = 0;
    rec1(0, N);
    cout << dp[N-1] << endl;
}