Fix library bugs

combinatorial_opt/convex_sum.hpp

@@ -10,7 +10,9 @@ struct Quadratic {
     Quadratic(Int a, Int b, Int c, Int lb, Int ub) : a(a), b(b), c(c), lb(lb), ub(ub) {}
     Int slope(Int s) const noexcept {
         if (a == 0) return b <= s ? ub : lb;
-        auto ret = (s + a - b) / (a * 2);
+        const Int num = s + a - b, den = a * 2;
+        auto ret = num / den;
+        if (num < 0 and num % den) --ret;
         return ret > ub ? ub : ret < lb ? lb : ret;
     }
     Int eval(Int x) const noexcept { return (a * x + b) * x + c; }

combinatorial_opt/matroids/binary_matroid.hpp

@@ -9,8 +9,17 @@
 // Verified: CF102156D 2019 Petrozavodsk Winter Camp, Yandex Cup D. Pick Your Own Nim
 template <int VDIM> class BinaryMatroid {
     using Element = int;
+    static int find_first(const std::bitset<VDIM> &bits) {
+#ifdef __GLIBCXX__
+        return bits._Find_first();
+#else
+        for (int i = 0; i < VDIM; i++)
+            if (bits[i]) return i;
+        return VDIM;
+#endif
+    }
     static void chxormin(std::bitset<VDIM> &l, const std::bitset<VDIM> &r) {
-        int i = r._Find_first();
+        int i = find_first(r);
         if (i < VDIM and l[i]) l ^= r;
     }
     std::vector<std::bitset<VDIM>> mat;

convolution/fft_arbitrary_mod.hpp

@@ -63,6 +63,7 @@ void fft(int n, std::vector<cmplx> &a) {
 // retval[i] = \sum_j a[j] b[i - j]
 template <typename MODINT>
 std::vector<MODINT> convolution_mod(std::vector<MODINT> a, std::vector<MODINT> b) {
+    if (a.empty() or b.empty()) return {};
     int need = int(a.size() + b.size()) - 1;
     int nbase = 0;
     while ((1 << nbase) < need) nbase++;
@@ -100,7 +101,7 @@ std::vector<MODINT> convolution_mod(std::vector<MODINT> a, std::vector<MODINT> b
     }
     fft(sz, fa);
     fft(sz, fb);
-    std::vector<MODINT> ret(sz);
+    std::vector<MODINT> ret(need);
     long long bp = MODINT(2).pow(D_FFT).val();
     long long cp = MODINT(2).pow(D_FFT * 2).val();
     for (int i = 0; i < need; i++) {

convolution/fft_double.hpp

@@ -52,6 +52,7 @@ std::vector<cmplx> conv_cmplx(const std::vector<T> &a, const std::vector<T> &b)
 // Requirement: length * max(a) * max(b) < 10^15
 template <typename T>
 std::vector<long long int> fftconv(const std::vector<T> &a, const std::vector<T> &b) {
+    if (a.empty() or b.empty()) return {};
     std::vector<cmplx> ans = conv_cmplx(a, b);
     std::vector<long long int> ret(ans.size());
     for (int i = 0; i < (int)ans.size(); i++) ret[i] = floor(ans[i].real() + 0.5);

data_structure/fibonacci_heap.hpp

@@ -1,8 +1,9 @@
 #pragma once
-#include <array>
 #include <cassert>
 #include <iostream>
 #include <list>
+#include <utility>
+#include <vector>
 
 // CUT begin
 // Fibonacci heap
@@ -52,9 +53,10 @@ template <typename Tp> struct fibonacci_heap {
     bool empty() const noexcept { return sz == 0; }
     int size() const noexcept { return sz; }
 
-    std::array<Node *, 30> _arr;
+    std::vector<Node *> _arr;
     void _fmerge(Node *ptr) {
         int d = ptr->deg;
+        if (d >= int(_arr.size())) _arr.resize(d + 1, nullptr);
         if (_arr[d] == nullptr)
             _arr[d] = ptr;
         else {
@@ -74,7 +76,7 @@ template <typename Tp> struct fibonacci_heap {
         }
     }
     void _consolidate() {
-        _arr.fill(nullptr);
+        _arr.assign(1, nullptr);
         for (auto ptr : roots)
             if (ptr != nullptr) {
                 if (ptr->deg < 0)
@@ -125,12 +127,14 @@ template <typename Tp> struct fibonacci_heap {
     }
 
     void _deldfs(Node *now) {
-        while (now != nullptr) {
+        if (now == nullptr) return;
+        Node *start = now;
+        do {
             if (now->child != nullptr) _deldfs(now->child);
             Node *nxt = now->right;
             delete now;
             now = nxt;
-        }
+        } while (now != start);
     }
     void clear() {
         for (auto root : roots) _deldfs(root);
@@ -198,8 +202,6 @@ template <typename Tp> struct fibonacci_heap {
     }
 };
 
-#include <utility>
-#include <vector>
 template <typename Tp> struct heap {
     using P = std::pair<Tp, int>;
     fibonacci_heap<P> _heap;
@@ -229,6 +231,7 @@ template <typename Tp> struct heap {
     P pop() {
         P ret = _heap.top();
         _heap.pop();
+        vp[ret.second] = nullptr;
         return ret;
     }
     int size() { return _heap.size(); }

data_structure/kd_tree_2d.hpp

@@ -59,8 +59,8 @@ template <class T> struct kd_tree {
                 // split y
                 std::nth_element(
                     _tmp.begin() + l, _tmp.begin() + c, _tmp.begin() + r,
-                    [&](const Tpl &l, const Tpl &r) { return std::get<2>(l) < std : get<2>(r); });
-                _nodes[_node_id].lch = _build(l, c, 0, nsplity + 1);:
+                    [&](const Tpl &l, const Tpl &r) { return std::get<2>(l) < std::get<2>(r); });
+                _nodes[_node_id].lch = _build(l, c, 0, nsplity + 1);
                 _nodes[_node_id].rch = _build(c, r, 0, nsplity + 1);
             }
         }

data_structure/persistent_queue.hpp

@@ -9,10 +9,10 @@
 template <typename T, int D> struct persistent_queue {
     int now;
 
-    std::vector<T> data;                 // Elements on each node of tree
-    std::vector<std::array<int, D>> par; // binary-lifted parents
-    std::vector<int> back_id;            // back_id[t] = leaf id of the tree at time t
-    std::vector<int> size;               // size[t] = size of the queue at time t
+    std::vector<T> data;                     // Elements on each node of tree
+    std::vector<std::array<int, D + 1>> par; // binary-lifted parents
+    std::vector<int> back_id;                // back_id[t] = leaf id of the tree at time t
+    std::vector<int> size;                   // size[t] = size of the queue at time t
 
     persistent_queue() : now(0), data(1), par(1), back_id(1, 0), size(1, 0) {}
 
@@ -23,7 +23,7 @@ template <typename T, int D> struct persistent_queue {
         assert(now < 1 << (D + 1));
         int r = back_id[t], len = size[t] - 1;
         back_id.emplace_back(r), size.emplace_back(len);
-        for (int d = 0; d < D; d++)
+        for (int d = 0; d <= D; d++)
             if ((len >> d) & 1) r = par[r][d];
         return std::make_pair(now, data[r]);
     }
@@ -37,7 +37,7 @@ template <typename T, int D> struct persistent_queue {
         data.emplace_back(dat);
         par.push_back({}), par.back()[0] = back_id[t];
         back_id.emplace_back(newid), size.emplace_back(size[t] + 1);
-        for (int d = 1; d < D; d++) par[newid][d] = par[par[newid][d - 1]][d - 1];
+        for (int d = 1; d <= D; d++) par[newid][d] = par[par[newid][d - 1]][d - 1];
         return now;
     }
 };

data_structure/radix_heap_array.hpp

@@ -73,5 +73,8 @@ template <class Uint> class radix_heap_array {
         chmin(vnew, i);
     }
 
-    void clear() noexcept { sz = 0, last = 0, i2bj.clear(); }
+    void clear() noexcept {
+        sz = 0, last = 0, i2bj.clear();
+        for (auto &bucket : v) bucket.clear();
+    }
 };

data_structure/static_range_inversion.hpp

@@ -71,6 +71,7 @@ template <typename T> struct StaticRangeInversion {
     }
     long long get(int l, int r) const {
         assert(l >= 0 and l <= N and r >= 0 and r <= N and l <= r);
+        if (l == r) return 0;
         const int lb = (l + bs - 1) / bs, rb = (r == N ? nb_bc : r / bs) - 1;
         long long ret = 0;
         if (l / bs == (r - 1) / bs) {

formal_power_series/coeff_of_rational_function.hpp

@@ -12,6 +12,7 @@ Tp coefficient_of_rational_function(long long N, std::vector<Tp> num, std::vecto
     assert(N >= 0);
     while (den.size() and den.back() == 0) den.pop_back();
     assert(den.size());
+    if (num.empty()) return Tp(0);
     int h = 0;
     while (den[h] == 0) h++;
     N += h;

formal_power_series/linear_recurrence.hpp

@@ -67,9 +67,10 @@ std::vector<Tfield> monomial_mod_polynomial(long long N, const std::vector<Tfiel
     assert(!f_reversed.empty() and f_reversed[0] == 1);
     int K = f_reversed.size() - 1;
     if (!K) return {};
-    int D = 64 - __builtin_clzll(N);
     std::vector<Tfield> ret(K, 0);
     ret[0] = 1;
+    if (N == 0) return ret;
+    int D = 64 - __builtin_clzll(N);
     auto self_conv = [](std::vector<Tfield> x) -> std::vector<Tfield> {
         int d = x.size();
         std::vector<Tfield> ret(d * 2 - 1);

formal_power_series/multipoint_evaluation.hpp

@@ -13,7 +13,7 @@ template <typename Tfield> struct MultipointEvaluation {
     using polynomial = FormalPowerSeries<Tfield>;
     std::vector<polynomial> segtree;
     MultipointEvaluation(const std::vector<Tfield> &xs) : nx(xs.size()) {
-        segtree.resize(nx * 2 - 1);
+        segtree.resize(nx ? nx * 2 - 1 : 0);
         for (int i = 0; i < nx; i++) { segtree[nx - 1 + i] = {-xs[i], 1}; }
         for (int i = nx - 2; i >= 0; i--) { segtree[i] = segtree[2 * i + 1] * segtree[2 * i + 2]; }
     }
@@ -29,6 +29,7 @@ template <typename Tfield> struct MultipointEvaluation {
         _eval_rec(f, 2 * now + 2);
     }
     std::vector<Tfield> evaluate_polynomial(const polynomial &f) {
+        if (!nx) return {};
         ret.resize(nx);
         _eval_rec(f, 0);
         return ret;
@@ -47,6 +48,7 @@ template <typename Tfield> struct MultipointEvaluation {
     }
     std::vector<Tfield> polynomial_interpolation(std::vector<Tfield> ys) {
         assert(nx == int(ys.size()));
+        if (!nx) return {};
         if (_interpolate_coeffs.empty()) {
             _interpolate_coeffs = evaluate_polynomial(segtree[0].differential());
             for (auto &x : _interpolate_coeffs) x = x.inv();

geometry/geometry.hpp

@@ -79,7 +79,7 @@ int ccw(const Point2d<T_P> &a, const Point2d<T_P> &b, const Point2d<T_P> &c) {
     if (v1.det(v2) > Point2d<T_P>::EPS) return 1;   // 左折
     if (v1.det(v2) < -Point2d<T_P>::EPS) return -1; // 右折
     if (v1.dot(v2) < -Point2d<T_P>::EPS) return 2;  // c-a-b
-    if (v1.norm() < v2.norm()) return -2;           // a-b-c
+    if (v1.norm2() < v2.norm2()) return -2;         // a-b-c
     return 0;                                       // a-c-b
 }
 
@@ -184,6 +184,7 @@ IntersectTwoCircles(const Point2d<T_P> &Ca, T_P Ra, const Point2d<T_P> &Cb, T_P
     static_assert(std::is_floating_point<T_P>::value == true);
     T_P d = (Ca - Cb).norm();
     if (Ra + Rb < d) return {};
+    if (d == 0) return {};
     T_P rc = (d * d + Ra * Ra - Rb * Rb) / (2 * d);
     T_P rs2 = Ra * Ra - rc * rc;
     if (rs2 < 0) return {};
@@ -199,7 +200,9 @@ std::vector<Float> IntersectCircleLine(const PointNd &x0, const PointNd &v, Floa
     Float b = Float(x0.dot(v)) / v.norm2();
     Float c = Float(x0.norm2() - Float(R) * R) / v.norm2();
     if (b * b - c < 0) return {};
-    Float ret1 = -b + sqrtl(b * b - c) * (b > 0 ? -1 : 1);
+    Float discriminant_root = sqrtl(b * b - c);
+    if (discriminant_root == 0) return {-b, -b};
+    Float ret1 = -b + discriminant_root * (b > 0 ? -1 : 1);
     Float ret2 = c / ret1;
     return ret1 < ret2 ? std::vector<Float>{ret1, ret2} : std::vector<Float>{ret2, ret1};
 }

graph/directed_mst.hpp

@@ -92,6 +92,7 @@ template <typename Weight = long long> struct MinimumSpanningArborescence {
         void pop() {
             data[root].push();
             root = _meld(data[root].r, data[root].l);
+            sz--;
         }
         int size() const { return sz; }
         bool empty() const { return sz == 0; }
@@ -162,7 +163,7 @@ template <typename Weight = long long> struct MinimumSpanningArborescence {
         }
     }
 };
-template <>
-std::vector<MinimumSpanningArborescence<long long>::skew_heap::node>
-    MinimumSpanningArborescence<long long>::skew_heap::data = {};
+template <typename Weight>
+std::vector<typename MinimumSpanningArborescence<Weight>::skew_heap::node>
+    MinimumSpanningArborescence<Weight>::skew_heap::data = {};
 template <typename T> unsigned MinimumSpanningArborescence<T>::skew_heap::len = 0;

graph/grid_graph_template.hpp

@@ -60,11 +60,12 @@ template <typename T_E, T_E INF, T_E (*edge_cost)(int, int, int, int)> struct Gr
             for (unsigned d = 0; d < dx.size(); d++) {
                 int xn = x + dx[d], yn = y + dy[d];
                 if (xn < 0 or yn < 0 or xn >= H or yn >= W) continue;
-                auto dnxt = dnow + edge_cost(x, y, xn, yn);
+                auto weight = edge_cost(x, y, xn, yn);
+                auto dnxt = dnow + weight;
                 if (dnxt < dist[xn][yn]) {
                     dist[xn][yn] = dnxt;
                     prv[xn][yn] = std::make_pair(x, y);
-                    if (dnxt)
+                    if (weight)
                         deq.emplace_back(xn, yn);
                     else
                         deq.emplace_front(xn, yn);

graph/maximum_independent_set.hpp

@@ -10,6 +10,24 @@
 // Verified: https://judge.yosupo.jp/submission/1864 / https://yukicoder.me/problems/no/382
 // Reference: https://www.slideshare.net/wata_orz/ss-12131479
 template <int BS = 64> struct maximum_independent_set {
+    static int find_first(const std::bitset<BS> &bits) {
+#ifdef __GLIBCXX__
+        return bits._Find_first();
+#else
+        for (int i = 0; i < BS; i++)
+            if (bits[i]) return i;
+        return BS;
+#endif
+    }
+    static int find_next(const std::bitset<BS> &bits, int pos) {
+#ifdef __GLIBCXX__
+        return bits._Find_next(pos);
+#else
+        for (int i = pos + 1; i < BS; i++)
+            if (bits[i]) return i;
+        return BS;
+#endif
+    }
     std::vector<std::bitset<BS>> conn;
     int V;               // # of vertices
     int nret;            // Largest possible size of independent set
@@ -22,13 +40,13 @@ template <int BS = 64> struct maximum_independent_set {
         std::stack<int> st;
         while (retry) {
             retry = false;
-            for (int i = _avail._Find_first(); i < V; i = _avail._Find_next(i)) {
+            for (int i = find_first(_avail); i < V; i = find_next(_avail, i)) {
                 int nb = (_avail & conn[i]).count();
                 if (nb <= 1) {
                     st.emplace(i), _avail.reset(i), _tmp_state.set(i);
                     retry = true;
                     if (nb == 1) {
-                        int j = (_avail & conn[i])._Find_first();
+                        int j = find_first(_avail & conn[i]);
                         st.emplace(j), _avail.reset(j);
                     }
                 }
@@ -39,7 +57,7 @@ template <int BS = 64> struct maximum_independent_set {
         if (t > nret) nret = t, ret = _tmp_state;
 
         int d = -1, n = -1;
-        for (int i = _avail._Find_first(); i < V; i = _avail._Find_next(i)) {
+        for (int i = find_first(_avail); i < V; i = find_next(_avail, i)) {
             int c = (_avail & conn[i]).count();
             if (c > d) d = c, n = i;
         }