This documentation is automatically generated by online-judge-tools/verification-helper
View the Project on GitHub maspypy/library
#include "convex/monge.hpp"
#include "convex/larsch.hpp" #include "convex/smawk.hpp" #include "other/fibonacci_search.hpp" // 定義域 [0, N] の範囲で f の monge 性を確認 template <typename T, typename F> bool check_monge(int N, F f) { FOR(l, N + 1) FOR(k, l) FOR(j, k) FOR(i, j) { T lhs = f(i, l) + f(j, k); T rhs = f(i, k) + f(j, l); if (lhs < rhs) { print("monge ng"); print(i, j, k, l, f(i, k), f(i, l), f(j, k), f(j, l), lhs, rhs); return false; } } print("monge ok"); return true; } // newdp[j] = min (dp[i] + f(i,j)) template <typename T, typename F> vc<T> monge_dp_update(int N, vc<T>& dp, F f) { assert(len(dp) == N + 1); auto select = [&](int i, int j, int k) -> int { if (i <= k) return j; return (dp[j] + f(j, i) > dp[k] + f(k, i) ? k : j); }; vc<int> I = SMAWK(N + 1, N + 1, select); vc<T> newdp(N + 1, infty<T>); FOR(j, N + 1) { int i = I[j]; chmin(newdp[j], dp[i] + f(i, j)); } return newdp; } // 遷移回数を問わない場合 template <typename T, typename F> vc<T> monge_shortest_path(int N, F f) { vc<T> dp(N + 1, infty<T>); dp[0] = 0; LARSCH<T> larsch(N, [&](int i, int j) -> T { ++i; if (i <= j) return infty<T>; return dp[j] + f(j, i); }); FOR(r, 1, N + 1) { int l = larsch.get_argmin(); dp[r] = dp[l] + f(l, r); } return dp; } // https://noshi91.github.io/algorithm-encyclopedia/d-edge-shortest-path-monge // |f| の上限 f_lim も渡す // ・larsch が結構重いので、自前で dp できるならその方がよい // ・複数の d で計算するとき:同じ lambda // に対する計算をメモ化しておくと定数倍高速? ・ABC305 template <typename T, typename F> T monge_shortest_path_d_edge(int N, int d, T f_lim, F f) { assert(d <= N); auto calc_L = [&](T lambda) -> T { auto cost = [&](int frm, int to) -> T { return f(frm, to) + lambda; }; vc<T> dp = monge_shortest_path<T>(N, cost); return dp[N] - lambda * d; }; auto [x, fx] = fibonacci_search<T, false>(calc_L, -3 * f_lim, 3 * f_lim + 1); return fx; } // https://topcoder-g-hatena-ne-jp.jag-icpc.org/spaghetti_source/20120915/1347668163.html // Prop 1 // 上三角 monge A, B // C[i][j] = min_k (A[i][k] + B[k][j]) template <typename T, typename F1, typename F2> vvc<T> monge_matrix_product(int N, F1 A, F2 B) { vv(T, C, N + 1, N + 1, infty<T>); vc<int> K(N + 1); FOR(i, N + 1) C[i][i] = A(i, i) + B(i, i), K[i] = i; FOR(s, 1, N + 1) { vc<int> newK(N + 1 - s); FOR(i, N + 1 - s) { int j = i + s; int p = K[i], q = K[i + 1]; FOR(k, p, q + 1) if (chmin(C[i][j], A(i, k) + B(k, j))) newK[i] = k; } swap(K, newK); } return C; }
#line 2 "convex/larsch.hpp" // https://noshi91.github.io/Library/algorithm/larsch.cpp.html template <class T> class LARSCH { struct reduce_row; struct reduce_col; struct reduce_row { int n; std::function<T(int, int)> f; int cur_row; int state; std::unique_ptr<reduce_col> rec; reduce_row(int n_) : n(n_), f(), cur_row(0), state(0), rec() { const int m = n / 2; if (m != 0) { rec = std::make_unique<reduce_col>(m); } } void set_f(std::function<T(int, int)> f_) { f = f_; if (rec) { rec->set_f([&](int i, int j) -> T { return f(2 * i + 1, j); }); } } int get_argmin() { const int cur_row_ = cur_row; cur_row += 1; if (cur_row_ % 2 == 0) { const int prev_argmin = state; const int next_argmin = [&]() { if (cur_row_ + 1 == n) { return n - 1; } else { return rec->get_argmin(); } }(); state = next_argmin; int ret = prev_argmin; for (int j = prev_argmin + 1; j <= next_argmin; j += 1) { if (f(cur_row_, ret) > f(cur_row_, j)) { ret = j; } } return ret; } else { if (f(cur_row_, state) <= f(cur_row_, cur_row_)) { return state; } else { return cur_row_; } } } }; struct reduce_col { int n; std::function<T(int, int)> f; int cur_row; std::vector<int> cols; reduce_row rec; reduce_col(int n_) : n(n_), f(), cur_row(0), cols(), rec(n) {} void set_f(std::function<T(int, int)> f_) { f = f_; rec.set_f([&](int i, int j) -> T { return f(i, cols[j]); }); } int get_argmin() { const int cur_row_ = cur_row; cur_row += 1; const auto cs = [&]() -> std::vector<int> { if (cur_row_ == 0) { return {{0}}; } else { return {{2 * cur_row_ - 1, 2 * cur_row_}}; } }(); for (const int j: cs) { while ([&]() { const int size = cols.size(); return size != cur_row_ && f(size - 1, cols.back()) > f(size - 1, j); }()) { cols.pop_back(); } if (int(cols.size()) != n) { cols.push_back(j); } } return cols[rec.get_argmin()]; } }; std::unique_ptr<reduce_row> base; public: LARSCH(int n, std::function<T(int, int)> f) : base(std::make_unique<reduce_row>(n)) { base->set_f(f); } int get_argmin() { return base->get_argmin(); } }; #line 2 "convex/smawk.hpp" // select(i,j,k) は (i,j) -> (i,k) を行うかどうか // 残念ながら monotone minima より高速な場合が存在しない説がある // https://codeforces.com/contest/1423/problem/M template <typename F> vc<int> smawk(int H, int W, F select) { auto dfs = [&](auto& dfs, vc<int> X, vc<int> Y) -> vc<int> { int N = len(X); if (N == 0) return {}; vc<int> YY; for (auto&& y: Y) { while (len(YY)) { int py = YY.back(), x = X[len(YY) - 1]; if (!select(x, py, y)) break; YY.pop_back(); } if (len(YY) < len(X)) YY.eb(y); } vc<int> XX; FOR(i, 1, len(X), 2) XX.eb(X[i]); vc<int> II = dfs(dfs, XX, YY); vc<int> I(N); FOR(i, len(II)) I[i + i + 1] = II[i]; int p = 0; FOR(i, 0, N, 2) { int LIM = (i + 1 == N ? Y.back() : I[i + 1]); int best = Y[p]; while (Y[p] < LIM) { ++p; if (select(X[i], best, Y[p])) best = Y[p]; } I[i] = best; } return I; }; vc<int> X(H), Y(W); iota(all(X), 0), iota(all(Y), 0); return dfs(dfs, X, Y); } #line 1 "other/fibonacci_search.hpp" // returns: {fx, x} // [L, R) での極小値をひとつ求める、単峰は不要 template <typename T, bool MINIMIZE, typename F> pair<T, ll> fibonacci_search(F f, ll L, ll R) { assert(L < R); --R; ll a = L, b = L + 1, c = L + 2, d = L + 3; int n = 0; while (d < R) { b = c, c = d, d = b + c - a, ++n; } auto get = [&](ll x) -> T { if (R < x) return infty<T>; return (MINIMIZE ? f(x) : -f(x)); }; T ya = get(a), yb = get(b), yc = get(c), yd = get(d); // この中で極小ならば全体でも極小、を維持する FOR(n) { if (yb <= yc) { d = c, c = b, b = a + d - c; yd = yc, yc = yb, yb = get(b); } else { a = b, b = c, c = a + d - b; ya = yb, yb = yc, yc = get(c); } } ll x = a; T y = ya; if (chmin(y, yb)) x = b; if (chmin(y, yc)) x = c; if (chmin(y, yd)) x = d; if (MINIMIZE) return {y, x}; return {-y, x}; } #line 4 "convex/monge.hpp" // 定義域 [0, N] の範囲で f の monge 性を確認 template <typename T, typename F> bool check_monge(int N, F f) { FOR(l, N + 1) FOR(k, l) FOR(j, k) FOR(i, j) { T lhs = f(i, l) + f(j, k); T rhs = f(i, k) + f(j, l); if (lhs < rhs) { print("monge ng"); print(i, j, k, l, f(i, k), f(i, l), f(j, k), f(j, l), lhs, rhs); return false; } } print("monge ok"); return true; } // newdp[j] = min (dp[i] + f(i,j)) template <typename T, typename F> vc<T> monge_dp_update(int N, vc<T>& dp, F f) { assert(len(dp) == N + 1); auto select = [&](int i, int j, int k) -> int { if (i <= k) return j; return (dp[j] + f(j, i) > dp[k] + f(k, i) ? k : j); }; vc<int> I = SMAWK(N + 1, N + 1, select); vc<T> newdp(N + 1, infty<T>); FOR(j, N + 1) { int i = I[j]; chmin(newdp[j], dp[i] + f(i, j)); } return newdp; } // 遷移回数を問わない場合 template <typename T, typename F> vc<T> monge_shortest_path(int N, F f) { vc<T> dp(N + 1, infty<T>); dp[0] = 0; LARSCH<T> larsch(N, [&](int i, int j) -> T { ++i; if (i <= j) return infty<T>; return dp[j] + f(j, i); }); FOR(r, 1, N + 1) { int l = larsch.get_argmin(); dp[r] = dp[l] + f(l, r); } return dp; } // https://noshi91.github.io/algorithm-encyclopedia/d-edge-shortest-path-monge // |f| の上限 f_lim も渡す // ・larsch が結構重いので、自前で dp できるならその方がよい // ・複数の d で計算するとき:同じ lambda // に対する計算をメモ化しておくと定数倍高速? ・ABC305 template <typename T, typename F> T monge_shortest_path_d_edge(int N, int d, T f_lim, F f) { assert(d <= N); auto calc_L = [&](T lambda) -> T { auto cost = [&](int frm, int to) -> T { return f(frm, to) + lambda; }; vc<T> dp = monge_shortest_path<T>(N, cost); return dp[N] - lambda * d; }; auto [x, fx] = fibonacci_search<T, false>(calc_L, -3 * f_lim, 3 * f_lim + 1); return fx; } // https://topcoder-g-hatena-ne-jp.jag-icpc.org/spaghetti_source/20120915/1347668163.html // Prop 1 // 上三角 monge A, B // C[i][j] = min_k (A[i][k] + B[k][j]) template <typename T, typename F1, typename F2> vvc<T> monge_matrix_product(int N, F1 A, F2 B) { vv(T, C, N + 1, N + 1, infty<T>); vc<int> K(N + 1); FOR(i, N + 1) C[i][i] = A(i, i) + B(i, i), K[i] = i; FOR(s, 1, N + 1) { vc<int> newK(N + 1 - s); FOR(i, N + 1 - s) { int j = i + s; int p = K[i], q = K[i + 1]; FOR(k, p, q + 1) if (chmin(C[i][j], A(i, k) + B(k, j))) newK[i] = k; } swap(K, newK); } return C; }