Fix library bugs 2

combinatorial_opt/matroid_union.hpp

@@ -1,5 +1,4 @@
 #pragma once
-#include "../graph/shortest_path.hpp"
 #include <cassert>
 #include <utility>
 #include <vector>
@@ -16,37 +15,54 @@ template <class M1, class M2, class State1, class State2, class T = int>
 bool augment_union_matroid(M1 &matroid1, M2 &matroid2, State1 &I1, State2 &I2,
                            const std::vector<T> &weights) {
     const int M = matroid1.size();
-    const int gs = M, gt = M + 1;
-    shortest_path<T> sssp(M + 2);
+    const int gt = M;
+    std::vector<std::vector<int>> rev(M + 1);
     std::vector<int> 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<int> 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<int> 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;
     }

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 個以上のマトロイドの合併に対しても同様のアルゴリズムが適用可能である．
 

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;
         }

data_structure/wavelet_matrix.hpp

@@ -296,14 +296,16 @@ template <class Int> 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<Int> 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 Int> 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<Int> 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:

number/discrete_logarithm.hpp

@@ -53,6 +53,9 @@ template <class S, class F, class Container>
 long long
 DiscreteLogarithm(const F &f, const S &s, const S &t, const std::function<S(F, S)> &mapping,
                   const std::function<F(F, F)> &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) {

other_algorithms/bisect.hpp

@@ -2,16 +2,18 @@
 #include <bit>
 #include <functional>
 #include <numeric>
+#include <type_traits>
 
 // 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<unsigned long long>(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<double>(tmp);
     };
 
@@ -29,7 +31,12 @@ template <class T> auto bisect(T ok, T ng, const std::function<bool(T)> &f, T ab
     };
 
     while (true) {
-        T mid = std::is_floating_point<T>::value ? bisect_mid_fp(ok, ng) : std::midpoint(ok, ng);
+        T mid;
+        if constexpr (std::is_floating_point<T>::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;

other_algorithms/doubling.hpp

@@ -13,7 +13,7 @@ struct BinaryLifting {
     std::vector<std::vector<int>> mat;
     BinaryLifting() : N(0), lgD(0) {}
     BinaryLifting(const std::vector<int> &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<int>(N));
         mat[0] = to;
 

string/palindromic_tree.hpp

@@ -24,7 +24,7 @@ template <class Key> 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 <class OStream> friend OStream &operator<<(OStream &os, const Node &node) {
         os << "Node(suffix_link=" << node.suffix_link() << ", length=" << node.length()

utilities/integer_segments.hpp

@@ -19,7 +19,10 @@ template <typename Int> 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 {