diff --git a/combinatorial_opt/matroid_union.hpp b/combinatorial_opt/matroid_union.hpp index 20ed55eb..aad6dd7a 100644 --- a/combinatorial_opt/matroid_union.hpp +++ b/combinatorial_opt/matroid_union.hpp @@ -1,5 +1,4 @@ #pragma once -#include "../graph/shortest_path.hpp" #include #include #include @@ -16,37 +15,54 @@ template bool augment_union_matroid(M1 &matroid1, M2 &matroid2, State1 &I1, State2 &I2, const std::vector &weights) { const int M = matroid1.size(); - const int gs = M, gt = M + 1; - shortest_path sssp(M + 2); + const int gt = M; + std::vector> rev(M + 1); std::vector color(M, -1); matroid1.set(I1); matroid2.set(I2); for (int e = 0; e < M; e++) { - if (!I1[e] and !I2[e]) sssp.add_edge(gs, e, weights.size() ? weights[e] : 0); if (!I1[e]) { auto c = matroid1.circuit(e); - if (c.empty()) sssp.add_edge(e, gt, 0), color[e] = 0; + if (c.empty()) rev[gt].push_back(e), color[e] = 0; for (int f : c) { - if (f != e) sssp.add_edge(e, f, 1); + if (f != e) rev[f].push_back(e); } } if (!I2[e]) { auto c = matroid2.circuit(e); - if (c.empty()) sssp.add_edge(e, gt, 0), color[e] = 1; + if (c.empty()) rev[gt].push_back(e), color[e] = 1; for (int f : c) { - if (f != e) sssp.add_edge(e, f, 1); + if (f != e) rev[f].push_back(e); } } } - sssp.solve(gs, gt); - auto aug_path = sssp.retrieve_path(gt); - if (aug_path.empty()) return false; - assert(aug_path.size() >= 3); + // Find a shortest exchange path to gt from every augmentable element. + std::vector next(M + 1, -1), q{gt}; + next[gt] = gt; + for (int i = 0; i < int(q.size()); ++i) { + for (int e : rev[q[i]]) { + if (next[e] >= 0) continue; + next[e] = q[i]; + q.push_back(e); + } + } + + // Minimize the added weight independently of the number of exchanges. + int start = -1; + for (int e = 0; e < M; ++e) { + if (I1[e] or I2[e] or next[e] < 0) continue; + if (start < 0 or (!weights.empty() and weights[e] < weights[start])) start = e; + } + if (start < 0) return false; + + std::vector aug_path; + for (int e = start; e != gt; e = next[e]) aug_path.push_back(e); int c0 = -1; - if (I1[aug_path[aug_path.size() - 2]]) c0 = 1; - if (I2[aug_path[aug_path.size() - 2]]) c0 = 0; - if (c0 < 0) c0 = color[aug_path[aug_path.size() - 2]]; - for (int k = int(aug_path.size()) - 2, e = aug_path[k]; k; e = aug_path[--k]) { + if (I1[aug_path.back()]) c0 = 1; + if (I2[aug_path.back()]) c0 = 0; + if (c0 < 0) c0 = color[aug_path.back()]; + for (int k = int(aug_path.size()) - 1; k >= 0; --k) { + int e = aug_path[k]; (c0 ? I2 : I1)[e] = 1, (c0 ? I1 : I2)[e] = 0; c0 ^= 1; } diff --git a/combinatorial_opt/matroid_union.md b/combinatorial_opt/matroid_union.md index d02301ac..523ffe70 100644 --- a/combinatorial_opt/matroid_union.md +++ b/combinatorial_opt/matroid_union.md @@ -3,19 +3,20 @@ title: Matroid union （マトロイドの合併） documentation_of: ./matroid_union.hpp --- -2つのマトロイド $M\_{1} = (E, \mathcal{I}\_{1}), M_{2} = (E, \mathcal{I}\_{2})$, $\mathcal{I}\_{1}$ に関して独立な集合 $I_1$, $\mathcal{I}\_{2}$ に関して独立な集合 $I\_2$ で $I\_1 \cup I\_2 = \emptyset$ を満たすものが与えられたとき，$I'\_1 + I'\_2 = I\_1 + I\_2 + \\{ e \\}$ を満たす新たな排反な独立集合 $I'\_1, I'\_2$ を見つけるアルゴリズム．特に重み最小の $e$ から貪欲に追加を試すことで，「合併したマトロイド」の最小重みサイズ $k$ 独立集合が $k = 1, 2, \dots$ について順次求められる． +2つのマトロイド $M\_{1} = (E, \mathcal{I}\_{1}), M_{2} = (E, \mathcal{I}\_{2})$, $\mathcal{I}\_{1}$ に関して独立な集合 $I_1$, $\mathcal{I}\_{2}$ に関して独立な集合 $I\_2$ で $I\_1 \cap I\_2 = \emptyset$ を満たすものが与えられたとき，$I'\_1 + I'\_2 = I\_1 + I\_2 + \\{ e \\}$ を満たす新たな排反な独立集合 $I'\_1, I'\_2$ を見つけるアルゴリズム．特に重み最小の $e$ から貪欲に追加を試すことで，「合併したマトロイド」の最小重みサイズ $k$ 独立集合が $k = 1, 2, \dots$ について順次求められる． これを応用すると，与えられた重み付き無向グラフについて，辺素な二つの全域木であって重みの総和が最小なものを見つけることも可能． ## アルゴリズムの概要 -$\|E\| = n$ として，$e = 1, \dots, n$ に $s$, $t$ を加えた $n + 2$ 頂点のグラフに次の要領で辺を張る： +$\|E\| = n$ として，各要素 $e \in E$ に終点 $t$ を加えた $n + 1$ 頂点の有向グラフに次の要領で辺を張る．すべての辺は重みなしとする． -- $e \notin I\_1 \cup I\_2$ のとき，$s \rightarrow e$ （重み $0$）を張る． -- $e \notin I\_j$ かつ $I\_j + \\{e \\}$ が $M\_j$ における独立集合のとき， $e \rightarrow t$ （重み $w(e)$）を張る．$(j = 1, 2)$ -- $e \notin I\_j$ かつ $I\_j + \\{e \\}$ が $M\_j$ における従属集合のとき，サーキットに含まれる各 $f$ について $e \rightarrow f$ （重み $1$）を張る．$(j=1, 2)$ +- $e \notin I\_j$ かつ $I\_j + \\{e \\}$ が $M\_j$ における独立集合のとき，$e \rightarrow t$ を張る．$(j = 1, 2)$ +- $e \notin I\_j$ かつ $I\_j + \\{e \\}$ が $M\_j$ における従属集合のとき，サーキットに含まれる $e$ 以外の各 $f$ について $e \rightarrow f$ を張る．$(j = 1, 2)$ -このグラフで $s$ から $t$ への最短路を求め，$s$ の次に通った要素が新たに追加される（$s$ から $t$ に到達不能ならば $I\_1 \cup I\_2$ は既に合併したマトロイド上の最大独立集合である）．それ以降に通った要素は既に $I\_1$ または $I\_2$ に含まれているが，これらを集合間で出し入れすることで $I\_1$ と $I\_2$ の独立性が保たれる． +辺を逆向きにしたグラフで $t$ から BFS を行い，各要素から $t$ への最短路を求める．$e \notin I\_1 \cup I\_2$ で $t$ に到達可能な要素のうち，重み $w(e)$ が最小のものを追加する．このとき重みと交換回数は足し合わせず，まず追加要素の重みを最小化し，その要素からの交換回数を最小化する．該当する要素がなければ $I\_1 \cup I\_2$ は既に合併したマトロイド上の最大独立集合である． + +選んだ $e$ から $t$ への最短路上で，要素を $I\_1$ と $I\_2$ の間で移動させる．最短路を使うことで独立性が保たれ，合併集合には $e$ だけが新たに加わる．重みを省略した場合は，追加可能ないずれかの要素を選ぶ． 本コードには実装されていないが，3 個以上のマトロイドの合併に対しても同様のアルゴリズムが適用可能である． diff --git a/convex_hull_trick/convex_hull_trick.hpp b/convex_hull_trick/convex_hull_trick.hpp index ebdb272d..2cff9540 100644 --- a/convex_hull_trick/convex_hull_trick.hpp +++ b/convex_hull_trick/convex_hull_trick.hpp @@ -30,9 +30,9 @@ struct Line { } bool operator<(const Line &r) const { if (b == T_MIN) { - return r.rp.first == T_MIN ? true : a * r.rp.second < r.rp.first; + return r.rp.first == T_MIN ? true : __int128(a) * r.rp.second < r.rp.first; } else if (r.b == T_MIN) { - return rp.first == T_MIN ? false : !(r.a * rp.second < rp.first); + return rp.first == T_MIN ? false : !(__int128(r.a) * rp.second < rp.first); } else { return a < r.a; } diff --git a/data_structure/wavelet_matrix.hpp b/data_structure/wavelet_matrix.hpp index 3cd30bcf..cee6173d 100644 --- a/data_structure/wavelet_matrix.hpp +++ b/data_structure/wavelet_matrix.hpp @@ -296,14 +296,16 @@ template class wavelet_matrix { if (upper_bound <= 0) return std::nullopt; const int n = index_range_freq(l, r, upper_bound); - return n == 0 ? std::nullopt : index_kth_smallest(l, r, n - 1); + if (n == 0) return std::nullopt; + return index_kth_smallest(l, r, n - 1); } // max y s.t. x in [xl, xr), y < yr std::optional prev_value(Int xl, Int xr, Int yr) const { const int l = to_index_x(xl), r = to_index_x(xr), ub = to_index_y(yr); const auto idx = index_prev_value(l, r, ub); - return idx ? distinct_ys.at(*idx) : std::nullopt; + if (!idx) return std::nullopt; + return distinct_ys.at(*idx); } // min v_i s.t. i in [l, r), v_i >= lower_bound @@ -312,14 +314,16 @@ template class wavelet_matrix { assert(is_built()); if (lower_bound >= (int)distinct_ys.size()) return std::nullopt; const int n = index_range_freq(l, r, lower_bound); - return n >= (r - l) ? std::nullopt : index_kth_smallest(l, r, n); + if (n >= r - l) return std::nullopt; + return index_kth_smallest(l, r, n); } // min y s.t. x in [xl, xr), y >= yl std::optional next_value(Int l, Int r, Int yl) const { const int xl = to_index_x(l), xr = to_index_x(r), yl_idx = to_index_y(yl); const auto idx = index_next_value(xl, xr, yl_idx); - return idx ? distinct_ys.at(*idx) : std::nullopt; + if (!idx) return std::nullopt; + return distinct_ys.at(*idx); } }; /* Sample usage: diff --git a/number/discrete_logarithm.hpp b/number/discrete_logarithm.hpp index f807d1c6..85ceb9d2 100644 --- a/number/discrete_logarithm.hpp +++ b/number/discrete_logarithm.hpp @@ -53,6 +53,9 @@ template long long DiscreteLogarithm(const F &f, const S &s, const S &t, const std::function &mapping, const std::function &composition, long long max_search) { + if (s == t) return 0; + if (max_search <= 0) return -1; + const int giant_stride = ceil(sqrtl(max_search)); F giant = f, tmp = f; for (int n = giant_stride - 1; n; n >>= 1) { diff --git a/other_algorithms/bisect.hpp b/other_algorithms/bisect.hpp index 48c3aa75..08680d39 100644 --- a/other_algorithms/bisect.hpp +++ b/other_algorithms/bisect.hpp @@ -2,16 +2,18 @@ #include #include #include +#include // Calculate next point to check in floating point "binary" search double bisect_mid_fp(double a, double b) { auto encode = [&](double x) -> unsigned long long { auto tmp = std::bit_cast(x); - return x >= 0 ? (tmp ^ (1ULL << 63)) : ~tmp; + // Give -0.0 and +0.0 the same code, as they compare equal in bisect(). + return (tmp >> 63) ? (0ULL - tmp) : (tmp ^ (1ULL << 63)); }; auto decode = [&](unsigned long long x) -> double { - auto tmp = (x >> 63) ? (x ^ (1ULL << 63)) : ~x; + auto tmp = (x >> 63) ? (x ^ (1ULL << 63)) : (0ULL - x); return std::bit_cast(tmp); }; @@ -29,7 +31,12 @@ template auto bisect(T ok, T ng, const std::function &f, T ab }; while (true) { - T mid = std::is_floating_point::value ? bisect_mid_fp(ok, ng) : std::midpoint(ok, ng); + T mid; + if constexpr (std::is_floating_point::value) { + mid = bisect_mid_fp(ok, ng); + } else { + mid = std::midpoint(ok, ng); + } if (mid == ok or mid == ng) break; (f(mid) ? ok : ng) = mid; if (ok - ng <= abs_tol and ng - ok <= abs_tol) break; diff --git a/other_algorithms/doubling.hpp b/other_algorithms/doubling.hpp index bf0c9c2d..06f314a4 100644 --- a/other_algorithms/doubling.hpp +++ b/other_algorithms/doubling.hpp @@ -13,7 +13,7 @@ struct BinaryLifting { std::vector> mat; BinaryLifting() : N(0), lgD(0) {} BinaryLifting(const std::vector &to, int lgd = 0) : N(to.size()), lgD(lgd) { - while ((1LL << lgD) < N) lgD++; + while (lgD == 0 or (1LL << lgD) < N) lgD++; mat.assign(lgD, std::vector(N)); mat[0] = to; diff --git a/string/palindromic_tree.hpp b/string/palindromic_tree.hpp index 31626900..83342d52 100644 --- a/string/palindromic_tree.hpp +++ b/string/palindromic_tree.hpp @@ -24,7 +24,7 @@ template class Node { return (it == children.end()) ? -1 : it->second; } - void set_child(int c, int nxt_idx) { children[c] = nxt_idx; } + void set_child(Key c, int nxt_idx) { children[c] = nxt_idx; } template friend OStream &operator<<(OStream &os, const Node &node) { os << "Node(suffix_link=" << node.suffix_link() << ", length=" << node.length() diff --git a/utilities/integer_segments.hpp b/utilities/integer_segments.hpp index b2d48a58..95d0d7f2 100644 --- a/utilities/integer_segments.hpp +++ b/utilities/integer_segments.hpp @@ -19,7 +19,10 @@ template struct integer_segments { return *itr; } - bool contains(Int x) const { return lower_bound(x) == x; } + bool contains(Int x) const { + auto itr = mp.upper_bound(x); + return itr != mp.begin() and std::prev(itr)->second >= x; + } // Find the min. y in the set that satisfies x <= y Int lower_bound(Int x) const {