diff --git a/combinatorial_opt/convex_sum.hpp b/combinatorial_opt/convex_sum.hpp
index e32bf4dd..553ed36d 100644
--- a/combinatorial_opt/convex_sum.hpp
+++ b/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; }
diff --git a/combinatorial_opt/matroids/binary_matroid.hpp b/combinatorial_opt/matroids/binary_matroid.hpp
index 5a7eed84..2b59a6f0 100644
--- a/combinatorial_opt/matroids/binary_matroid.hpp
+++ b/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;
diff --git a/convolution/fft_arbitrary_mod.hpp b/convolution/fft_arbitrary_mod.hpp
index 36c3200c..0a8eb394 100644
--- a/convolution/fft_arbitrary_mod.hpp
+++ b/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++) {
diff --git a/convolution/fft_double.hpp b/convolution/fft_double.hpp
index 3e0859a8..6cdabd13 100644
--- a/convolution/fft_double.hpp
+++ b/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);
diff --git a/data_structure/fibonacci_heap.hpp b/data_structure/fibonacci_heap.hpp
index 33f8ab52..524a792e 100644
--- a/data_structure/fibonacci_heap.hpp
+++ b/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(); }
diff --git a/data_structure/kd_tree_2d.hpp b/data_structure/kd_tree_2d.hpp
index 417ad58e..54916f2d 100644
--- a/data_structure/kd_tree_2d.hpp
+++ b/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);
             }
         }
diff --git a/data_structure/persistent_queue.hpp b/data_structure/persistent_queue.hpp
index cdcd2927..6b4db591 100644
--- a/data_structure/persistent_queue.hpp
+++ b/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;
     }
 };
diff --git a/data_structure/radix_heap_array.hpp b/data_structure/radix_heap_array.hpp
index 4e0e2a13..ae6f5e42 100644
--- a/data_structure/radix_heap_array.hpp
+++ b/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();
+    }
 };
diff --git a/data_structure/static_range_inversion.hpp b/data_structure/static_range_inversion.hpp
index cb92cea8..e451b07f 100644
--- a/data_structure/static_range_inversion.hpp
+++ b/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) {
diff --git a/formal_power_series/coeff_of_rational_function.hpp b/formal_power_series/coeff_of_rational_function.hpp
index e2c4e281..72c00e8d 100644
--- a/formal_power_series/coeff_of_rational_function.hpp
+++ b/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;
diff --git a/formal_power_series/linear_recurrence.hpp b/formal_power_series/linear_recurrence.hpp
index 2b2f7bbb..5440e696 100644
--- a/formal_power_series/linear_recurrence.hpp
+++ b/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);
diff --git a/formal_power_series/multipoint_evaluation.hpp b/formal_power_series/multipoint_evaluation.hpp
index a727eaee..310b9d6b 100644
--- a/formal_power_series/multipoint_evaluation.hpp
+++ b/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();
diff --git a/geometry/geometry.hpp b/geometry/geometry.hpp
index d22d7cb0..9192d4b0 100644
--- a/geometry/geometry.hpp
+++ b/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};
 }
diff --git a/graph/directed_mst.hpp b/graph/directed_mst.hpp
index 79bc7758..1a806ae3 100644
--- a/graph/directed_mst.hpp
+++ b/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;
diff --git a/graph/grid_graph_template.hpp b/graph/grid_graph_template.hpp
index 300e7053..281d5f42 100644
--- a/graph/grid_graph_template.hpp
+++ b/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);
diff --git a/graph/maximum_independent_set.hpp b/graph/maximum_independent_set.hpp
index 56c293b1..35af1aff 100644
--- a/graph/maximum_independent_set.hpp
+++ b/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;
         }
diff --git a/graph/shortest_cycle01.hpp b/graph/shortest_cycle01.hpp
index be1e1c19..ad358595 100644
--- a/graph/shortest_cycle01.hpp
+++ b/graph/shortest_cycle01.hpp
@@ -1,4 +1,5 @@
 #pragma once
+#include <algorithm>
 #include <cassert>
 #include <deque>
 #include <limits>
@@ -35,43 +36,90 @@ struct ShortestCycle01 {
     std::pair<int, std::pair<int, int>> Solve(int v) {
         assert(0 <= v and v < V);
         dist.assign(V, INF);
-        dist[v] = 0;
         prev.assign(V, -1);
         orig.assign(V, -1);
-        std::deque<std::pair<int, int>> bfsq;
-        std::vector<std::pair<std::pair<int, int>, int>> add_edge;
-        bfsq.emplace_back(v, -1);
+        std::vector<int> root_weight(V, INF);
+        std::deque<int> bfsq;
+        for (auto [nxt, weight] : to[v]) {
+            if (weight >= dist[nxt]) continue;
+            dist[nxt] = weight;
+            root_weight[nxt] = weight;
+            prev[nxt] = v;
+            orig[nxt] = nxt;
+            if (weight)
+                bfsq.emplace_back(nxt);
+            else
+                bfsq.emplace_front(nxt);
+        }
+        int minimum_cycle = INF;
+        int source_a = -1, source_b = -1;
         while (!bfsq.empty()) {
-            int now = bfsq.front().first, prv = bfsq.front().second;
+            int now = bfsq.front();
             bfsq.pop_front();
-            if (prv < 0) {
-                //
-            } else if (prv == v) {
-                orig.at(now) = now;
-            } else {
-                orig.at(now) = orig.at(prv);
+            for (auto [nxt, weight] : to[now]) {
+                if (nxt == v) continue;
+                if (orig[nxt] >= 0 and orig[now] != orig[nxt]) {
+                    int length = dist[now] + dist[nxt] + weight;
+                    if (length < minimum_cycle) {
+                        minimum_cycle = length;
+                        source_a = orig[now], source_b = orig[nxt];
+                    }
+                }
+                int dnext = dist[now] + weight;
+                if (dnext < dist[nxt]) {
+                    dist[nxt] = dnext;
+                    prev[nxt] = now;
+                    orig[nxt] = orig[now];
+                    if (weight)
+                        bfsq.emplace_back(nxt);
+                    else
+                        bfsq.emplace_front(nxt);
+                }
             }
-            for (auto nxt : to[now])
-                if (nxt.first != prv) {
-                    if (dist[nxt.first] == INF) {
-                        dist[nxt.first] = dist[now] + nxt.second;
-                        prev[nxt.first] = now;
-                        if (nxt.second)
-                            bfsq.emplace_back(nxt.first, now);
-                        else
-                            bfsq.emplace_front(nxt.first, now);
-                    } else
-                        add_edge.emplace_back(std::make_pair(now, nxt.first), nxt.second);
+        }
+        if (source_a < 0) return std::make_pair(INF, std::make_pair(-1, -1));
+
+        std::vector<int> path_dist(V, INF), path_prev(V, -1);
+        path_dist[source_a] = 0;
+        bfsq = {source_a};
+        while (!bfsq.empty()) {
+            int now = bfsq.front();
+            bfsq.pop_front();
+            for (auto [nxt, weight] : to[now]) {
+                if (nxt == v) continue;
+                int dnext = path_dist[now] + weight;
+                if (dnext < path_dist[nxt]) {
+                    path_dist[nxt] = dnext;
+                    path_prev[nxt] = now;
+                    if (weight)
+                        bfsq.emplace_back(nxt);
+                    else
+                        bfsq.emplace_front(nxt);
                 }
+            }
         }
-        int minimum_cycle = INF;
-        int s = -1, t = -1;
-        for (auto edge : add_edge) {
-            int a = edge.first.first, b = edge.first.second;
-            if (orig.at(a) == orig.at(b)) continue;
-            int L = dist[a] + dist[b] + edge.second;
-            if (L < minimum_cycle) minimum_cycle = L, s = a, t = b;
+        std::vector<int> path;
+        for (int x = source_b; x >= 0; x = path_prev[x]) {
+            path.push_back(x);
+            if (x == source_a) break;
+        }
+        std::reverse(path.begin(), path.end());
+        int split = (path.size() - 1) / 2;
+        prev[path.front()] = v;
+        orig[path.front()] = path.front();
+        dist[path.front()] = root_weight[path.front()];
+        for (int i = 1; i <= split; i++) {
+            prev[path[i]] = path[i - 1];
+            orig[path[i]] = path.front();
+        }
+        prev[path.back()] = v;
+        orig[path.back()] = path.back();
+        dist[path.back()] = root_weight[path.back()];
+        for (int i = path.size() - 2; i > split; i--) {
+            prev[path[i]] = path[i + 1];
+            orig[path[i]] = path.back();
         }
-        return std::make_pair(minimum_cycle, std::make_pair(s, t));
+        minimum_cycle = root_weight[source_a] + path_dist[source_b] + root_weight[source_b];
+        return std::make_pair(minimum_cycle, std::make_pair(path[split], path[split + 1]));
     }
 };
diff --git a/graph/shortest_path.hpp b/graph/shortest_path.hpp
index 5df21fb0..61e25f44 100644
--- a/graph/shortest_path.hpp
+++ b/graph/shortest_path.hpp
@@ -154,13 +154,14 @@ struct shortest_path {
                 int nxt = nx.first;
                 if (dist[nxt] > dnx) {
                     dist[nxt] = dnx;
+                    prev[nxt] = now;
                     if (!in_queue[nxt]) {
                         if (q.size() and dnx < dist[q.front()]) { // Small label first optimization
                             q.push_front(nxt);
                         } else {
                             q.push_back(nxt);
                         }
-                        prev[nxt] = now, in_queue[nxt] = 1;
+                        in_queue[nxt] = 1;
                     }
                 }
             }
diff --git a/graph/strongly_connected_components.hpp b/graph/strongly_connected_components.hpp
index 827a81f0..578070f4 100644
--- a/graph/strongly_connected_components.hpp
+++ b/graph/strongly_connected_components.hpp
@@ -67,6 +67,9 @@ struct DirectedGraphSCC {
     }
     std::vector<int> DetectCycle() {
         int ns = FindStronglyConnectedComponents();
+        for (int v = 0; v < V; v++) {
+            if (std::find(to[v].begin(), to[v].end(), v) != to[v].end()) return {v};
+        }
         if (ns == V) return {};
         std::vector<int> cnt(ns);
         for (auto x : cmp) cnt[x]++;
diff --git a/graph/strongly_connected_components_bitset.hpp b/graph/strongly_connected_components_bitset.hpp
index 8b25d5ff..fc27da0a 100644
--- a/graph/strongly_connected_components_bitset.hpp
+++ b/graph/strongly_connected_components_bitset.hpp
@@ -11,6 +11,15 @@
 // Complexity: O(V^2/64)
 // Verified: CF1268D <https://codeforces.com/contest/1268/submission/68125495>
 template <int VMAX> struct DirectedGraphSCC64 {
+    static int find_first(const std::bitset<VMAX> &bits) {
+#ifdef __GLIBCXX__
+        return bits._Find_first();
+#else
+        for (int i = 0; i < VMAX; i++)
+            if (bits[i]) return i;
+        return VMAX;
+#endif
+    }
     int V;
     const std::vector<std::bitset<VMAX>> &e, &einv;
     std::vector<int> vs, cmp;
@@ -24,7 +33,7 @@ template <int VMAX> struct DirectedGraphSCC64 {
         while (!_st.empty()) {
             int now = _st.back();
             unvisited.reset(now);
-            int nxt = (unvisited & e[now])._Find_first();
+            int nxt = find_first(unvisited & e[now]);
             if (nxt < V) {
                 unvisited.reset(nxt);
                 _st.push_back(nxt);
@@ -43,7 +52,7 @@ template <int VMAX> struct DirectedGraphSCC64 {
             _st.pop_back();
             cmp[now] = k;
             while (true) {
-                int nxt = (unvisited & einv[now])._Find_first();
+                int nxt = find_first(unvisited & einv[now]);
                 if (nxt >= V) break;
                 _st.push_back(nxt);
                 unvisited.reset(nxt);
diff --git a/linear_algebra_matrix/circular_binary_expansion.hpp b/linear_algebra_matrix/circular_binary_expansion.hpp
index e423f077..e1030e62 100644
--- a/linear_algebra_matrix/circular_binary_expansion.hpp
+++ b/linear_algebra_matrix/circular_binary_expansion.hpp
@@ -31,6 +31,6 @@ std::pair<matrix<MODINT>, matrix<MODINT>> circular_binary_expansion(int lgdim) {
     }
     for (int i = 0; i < D; i++) {
         for (int j = 0; j < D; j++) transinv[i][j] *= invD;
-        return {trans, transinv};
     }
+    return {trans, transinv};
 }
diff --git a/linear_algebra_matrix/matrix.hpp b/linear_algebra_matrix/matrix.hpp
index ae76e796..e7aed541 100644
--- a/linear_algebra_matrix/matrix.hpp
+++ b/linear_algebra_matrix/matrix.hpp
@@ -191,24 +191,21 @@ template <typename T> struct matrix {
         assert(H == W);
         std::vector<std::vector<T>> ret = Identity(H), tmp = *this;
         int rank = 0;
-        for (int i = 0; i < H; i++) {
-            int ti = i;
-            while (ti < H and tmp[ti][i] == T()) ti++;
-            if (ti == H) {
-                continue;
-            } else {
-                rank++;
-            }
-            ret[i].swap(ret[ti]), tmp[i].swap(tmp[ti]);
-            T inv = _T_id<T>() / tmp[i][i];
-            for (int j = 0; j < W; j++) ret[i][j] *= inv;
-            for (int j = i + 1; j < W; j++) tmp[i][j] *= inv;
+        for (int c = 0; c < W; c++) {
+            int ti = rank;
+            while (ti < H and tmp[ti][c] == T()) ti++;
+            if (ti == H) { continue; }
+            ret[rank].swap(ret[ti]), tmp[rank].swap(tmp[ti]);
+            T inv = _T_id<T>() / tmp[rank][c];
+            for (int j = 0; j < W; j++) ret[rank][j] *= inv;
+            for (int j = c + 1; j < W; j++) tmp[rank][j] *= inv;
             for (int h = 0; h < H; h++) {
-                if (i == h) continue;
-                const T c = -tmp[h][i];
-                for (int j = 0; j < W; j++) ret[h][j] += ret[i][j] * c;
-                for (int j = i + 1; j < W; j++) tmp[h][j] += tmp[i][j] * c;
+                if (rank == h) continue;
+                const T coeff = -tmp[h][c];
+                for (int j = 0; j < W; j++) ret[h][j] += ret[rank][j] * coeff;
+                for (int j = c + 1; j < W; j++) tmp[h][j] += tmp[rank][j] * coeff;
             }
+            rank++;
         }
         *this = ret;
         return rank;
diff --git a/linear_algebra_matrix/test/upper_trinaglular_matrix.yuki3530.test.cpp b/linear_algebra_matrix/test/upper_trinaglular_matrix.yuki3530.test.cpp
index b7739dd4..40c9574f 100644
--- a/linear_algebra_matrix/test/upper_trinaglular_matrix.yuki3530.test.cpp
+++ b/linear_algebra_matrix/test/upper_trinaglular_matrix.yuki3530.test.cpp
@@ -13,10 +13,37 @@ using namespace std;
 using S = UpperTriangular3d<ModInt998244353>;
 
 S op(const S &l, const S &r) { return l * r; }
-S e() { return S{1, 0, 0, 1, 0, 1}; }
+S e() {
+    return S{
+        .a00 = 1,
+        .a01 = 0,
+        .a02 = 0,
+        .a11 = 1,
+        .a12 = 0,
+        .a22 = 1,
+    };
+}
 
-S GenR() { return S{ModInt998244353(3) / 4, ModInt998244353(1) / 4, 0, 1, 0, 1}; }
-S GenL() { return S{1, 0, 0, ModInt998244353(3) / 4, ModInt998244353(1) / 4, 1}; }
+S GenR() {
+    return S{
+        .a00 = ModInt998244353(3) / 4,
+        .a01 = ModInt998244353(1) / 4,
+        .a02 = 0,
+        .a11 = 1,
+        .a12 = 0,
+        .a22 = 1,
+    };
+}
+S GenL() {
+    return S{
+        .a00 = 1,
+        .a01 = 0,
+        .a02 = 0,
+        .a11 = ModInt998244353(3) / 4,
+        .a12 = ModInt998244353(1) / 4,
+        .a22 = 1,
+    };
+}
 
 ModInt998244353 Solve(vector<pair<int, int>> ps) {
     vector<tuple<int, int, int>> yxis;
diff --git a/linear_algebra_matrix/upper_triangular_matrix.hpp b/linear_algebra_matrix/upper_triangular_matrix.hpp
index d2d4cdcb..2e03559b 100644
--- a/linear_algebra_matrix/upper_triangular_matrix.hpp
+++ b/linear_algebra_matrix/upper_triangular_matrix.hpp
@@ -1,11 +1,31 @@
 #pragma once
 
+#include <compare>
+
 template <class T> struct UpperTriangular3d {
-    static T explicit_init_required() = delete;
-    T a00 = this->explicit_init_required(), a01 = this->explicit_init_required(),
-      a02 = this->explicit_init_required();
-    T a11 = this->explicit_init_required(), a12 = this->explicit_init_required();
-    T a22 = this->explicit_init_required();
+private:
+    struct DesignatedInitializationOnly {
+    private:
+        constexpr DesignatedInitializationOnly() = default;
+        constexpr DesignatedInitializationOnly(const DesignatedInitializationOnly &) = default;
+        DesignatedInitializationOnly &operator=(const DesignatedInitializationOnly &) = default;
+        friend UpperTriangular3d;
+
+    public:
+        auto operator<=>(const DesignatedInitializationOnly &) const = default;
+    };
+
+    template <class U> static constexpr U explicit_init_required() {
+        static_assert(sizeof(U) == 0, "all matrix entries must be explicitly initialized");
+        return U{};
+    }
+
+public:
+    [[no_unique_address]] DesignatedInitializationOnly designated_initialization_only{};
+    T a00 = explicit_init_required<T>(), a01 = explicit_init_required<T>(),
+      a02 = explicit_init_required<T>();
+    T a11 = explicit_init_required<T>(), a12 = explicit_init_required<T>();
+    T a22 = explicit_init_required<T>();
 
     UpperTriangular3d operator*(const UpperTriangular3d &r) const {
         return UpperTriangular3d{
diff --git a/number/arithmetic_cumsum.hpp b/number/arithmetic_cumsum.hpp
index 63a6f0f7..1d587b7e 100644
--- a/number/arithmetic_cumsum.hpp
+++ b/number/arithmetic_cumsum.hpp
@@ -156,7 +156,7 @@ template <class T> arithmetic_cumsum<T> zeta_shift_cumsum(long long n, int k) {
     auto init_pows = Sieve(k).enumerate_kth_pows<T>(k, k);
     for (int i = 1; i <= ret.K; ++i) {
         ret.a[i - 1] = T(i).pow(k);
-        ret.A[i] = ret.a[i] + (i ? ret.A[i - 1] : 0);
+        ret.A[i - 1] = ret.a[i - 1] + (i > 1 ? ret.A[i - 2] : 0);
     }
     for (int l = 0; l < ret.L; ++l) {
         ret.invA[l] = sum_of_exponential_times_polynomial<T>(1, init_pows, n / (l + 1) + 1);
diff --git a/number/binary_gcd.hpp b/number/binary_gcd.hpp
index e7a0ef54..1ae337ff 100644
--- a/number/binary_gcd.hpp
+++ b/number/binary_gcd.hpp
@@ -1,8 +1,14 @@
 #pragma once
+#include <type_traits>
 
 // CUT begin
 template <typename Int> Int binary_gcd(Int x_, Int y_) {
-    unsigned long long x = x_ < 0 ? -x_ : x_, y = y_ < 0 ? -y_ : y_;
+    using Uint = std::make_unsigned_t<Int>;
+    auto magnitude = [](Int v) -> Uint {
+        Uint u = static_cast<Uint>(v);
+        return v < 0 ? Uint(0) - u : u;
+    };
+    unsigned long long x = magnitude(x_), y = magnitude(y_);
     if (!x or !y) return x + y;
     int n = __builtin_ctzll(x), m = __builtin_ctzll(y);
     x >>= n, y >>= m;
diff --git a/number/modint_runtime.hpp b/number/modint_runtime.hpp
index 937d4679..7d69783c 100644
--- a/number/modint_runtime.hpp
+++ b/number/modint_runtime.hpp
@@ -45,7 +45,9 @@ struct ModIntRuntime {
         get_primitive_root() = 0;
     }
     ModIntRuntime &_setval(lint v) {
-        val_ = (v >= md ? v - md : v);
+        if (v < 0) v += md;
+        if (v >= md) v -= md;
+        val_ = v;
         return *this;
     }
     int val() const noexcept { return val_; }
diff --git a/number/pow_mod.hpp b/number/pow_mod.hpp
index 91f6fe0c..d693d5ae 100644
--- a/number/pow_mod.hpp
+++ b/number/pow_mod.hpp
@@ -10,7 +10,8 @@ template <class Int> Int pow_mod(Int x, long long n, Int md) {
     if (md == 1) return 0;
     if (n == 0) return 1;
 
-    x = (x % md + md) % md;
+    x %= md;
+    if (x < 0) x += md;
     Int ans = 1;
     while (n > 0) {
         if (n & 1) ans = (Long)ans * x % md;
diff --git a/number/sieve.hpp b/number/sieve.hpp
index 49075646..020a9a3e 100644
--- a/number/sieve.hpp
+++ b/number/sieve.hpp
@@ -103,7 +103,9 @@ struct Sieve {
         assert(K >= 0);
         if (K == 0) return std::vector<MODINT>(nmax + 1, 1);
         std::vector<MODINT> ret(nmax + 1);
-        ret[0] = 0, ret[1] = 1;
+        ret[0] = 0;
+        if (nmax == 0) return ret;
+        ret[1] = 1;
         for (int n = 2; n <= nmax; n++) {
             if (min_factor[n] == n) {
                 ret[n] = MODINT(n).pow(K);
diff --git a/number/sqrt_mod.hpp b/number/sqrt_mod.hpp
index 07f36dbb..5aa6f10b 100644
--- a/number/sqrt_mod.hpp
+++ b/number/sqrt_mod.hpp
@@ -17,9 +17,9 @@ template <class Int> Int sqrt_mod(Int a, Int p) {
         }
         return ans;
     };
+    a %= p;
+    if (a < 0) a += p;
     if (a == 0) return 0;
-
-    a = (a % p + p) % p;
     if (p == 2) return a;
     if (pow(a, (p - 1) / 2) != 1) return -1;
 
diff --git a/number/test/binary_gcd.stress.test.cpp b/number/test/binary_gcd.stress.test.cpp
index 17cda2cf..3a3e5537 100644
--- a/number/test/binary_gcd.stress.test.cpp
+++ b/number/test/binary_gcd.stress.test.cpp
@@ -18,6 +18,8 @@ template <typename Int> void test_binary_gcd(Int lo, Int hi) {
 }
 
 int main() {
+    test_binary_gcd<signed char>(-127, 126);
+    test_binary_gcd<short>(-1000, 1000);
     test_binary_gcd<int>(-1000, 1000);
     test_binary_gcd<unsigned int>(0, 2000);
     test_binary_gcd<long long>(-1000, 1000);
diff --git a/other_algorithms/binary_lifting.hpp b/other_algorithms/binary_lifting.hpp
index b080aa9f..b0b37e56 100644
--- a/other_algorithms/binary_lifting.hpp
+++ b/other_algorithms/binary_lifting.hpp
@@ -107,6 +107,7 @@ template <class S, S (*op)(S, S)> class binary_lifting {
             ++d;
         }
         if (d > maxd) return 1LL << maxd;
+        if (d == 0) return 0;
 
         --d;
 
diff --git a/other_algorithms/bounded_knapsack.hpp b/other_algorithms/bounded_knapsack.hpp
index 533572fe..fdebf030 100644
--- a/other_algorithms/bounded_knapsack.hpp
+++ b/other_algorithms/bounded_knapsack.hpp
@@ -102,7 +102,7 @@ std::pair<int, std::vector<bool>> bounded_knapsack(const std::vector<int> &weigh
     auto sol = bounded_knapsack_nonnegative(tmp, capacity);
     for (int i = 0; i < int(weights.size()); ++i) {
         if (weights.at(i) < 0) {
-            capacity += weights.at(i);
+            sol.first += weights.at(i);
             sol.second.at(i).flip();
         }
     }
diff --git a/other_algorithms/mos_algorithm.hpp b/other_algorithms/mos_algorithm.hpp
index 694bc570..d90bcb5c 100644
--- a/other_algorithms/mos_algorithm.hpp
+++ b/other_algorithms/mos_algorithm.hpp
@@ -1,6 +1,8 @@
 #pragma once
+#include <algorithm>
 #include <cassert>
 #include <cmath>
+#include <numeric>
 #include <utility>
 #include <vector>
 
diff --git a/segmenttree/binary_indexed_tree_2d.hpp b/segmenttree/binary_indexed_tree_2d.hpp
index d2e4b4aa..525729d3 100644
--- a/segmenttree/binary_indexed_tree_2d.hpp
+++ b/segmenttree/binary_indexed_tree_2d.hpp
@@ -3,9 +3,9 @@
 
 // 2-dimensional 1-indexed BIT (i : [1, lenX][1, lenY])
 template <typename T, int lenX, int lenY> struct BIT_2D {
-    std::array<T, (lenX + 1) * (lenY + 1)> val;
+    std::array<T, (lenX + 1) * (lenY + 1)> val{};
     constexpr static int M = lenY + 1;
-    BIT_2D() {}
+    BIT_2D() = default;
     void add(int posx, int posy, T v) noexcept {
         for (int x = posx; x <= lenX; x += x & -x) {
             for (int y = posy; y <= lenY; y += y & -y) val[x * M + y] += v;
diff --git a/segmenttree/point-update-range-get_nonrecursive.hpp b/segmenttree/point-update-range-get_nonrecursive.hpp
index a38d32d6..ea99344d 100644
--- a/segmenttree/point-update-range-get_nonrecursive.hpp
+++ b/segmenttree/point-update-range-get_nonrecursive.hpp
@@ -149,7 +149,9 @@ struct CountAndSumLessThan
         return ret;
     }
     TRET data2ret(const TDATA &vec, const TQUERY &q) override {
-        int i = std::lower_bound(vec.begin(), vec.end(), std::make_pair(q, q)) - vec.begin();
+        int i = std::lower_bound(vec.begin(), vec.end(), q,
+                                 [](const auto &p, const TQUERY &v) { return p.first < v; }) -
+                vec.begin();
         if (!i)
             return std::make_pair(0, 0);
         else
diff --git a/segmenttree/range-add-range-min.hpp b/segmenttree/range-add-range-min.hpp
index 5167f0a5..c9bc4b1f 100644
--- a/segmenttree/range-add-range-min.hpp
+++ b/segmenttree/range-add-range-min.hpp
@@ -50,5 +50,7 @@ template <typename Tp, Tp defaultT = std::numeric_limits<Tp>::max() / 2> struct
     }
     // Return f(x_begin, ..., x_{end - 1})
     Tp get(int pos) const noexcept { return prod(pos, pos + 1); }
-    Tp prod(int begin, int end) const noexcept { return _get(begin, end, 1, 0, head); }
+    Tp prod(int begin, int end) const noexcept {
+        return begin == end ? defaultT : _get(begin, end, 1, 0, head);
+    }
 };
diff --git a/string/aho_corasick.hpp b/string/aho_corasick.hpp
index 5dd4fdc8..10773217 100644
--- a/string/aho_corasick.hpp
+++ b/string/aho_corasick.hpp
@@ -18,7 +18,16 @@ template <class T, int (*char2int)(char)> struct AhoCorasick {
     const int D;
     std::vector<T> node;
     AhoCorasick(int D_) : built(false), D(D_), node(1, D) {}
-    AhoCorasick operator=(const AhoCorasick &rhs) { return AhoCorasick(rhs.D); }
+    AhoCorasick(const AhoCorasick &) = default;
+    AhoCorasick &operator=(const AhoCorasick &rhs) {
+        if (this == &rhs) return *this;
+        assert(D == rhs.D);
+        built = rhs.built;
+        node = rhs.node;
+        endpos = rhs.endpos;
+        visorder = rhs.visorder;
+        return *this;
+    }
 
     void enter_child(int n, int nn, int c) { node[n].setch(c, nn); }
 
diff --git a/string/aho_corasick_online.hpp b/string/aho_corasick_online.hpp
index 70196710..a741d598 100644
--- a/string/aho_corasick_online.hpp
+++ b/string/aho_corasick_online.hpp
@@ -1,5 +1,6 @@
 #pragma once
 #include "aho_corasick.hpp"
+#include <bit>
 #include <string>
 #include <vector>
 
@@ -17,7 +18,7 @@ struct OnlineAhoCorasick {
 
     // O(lg(n_keywords) |keyword|) amortized
     void add(const std::string &keyword) {
-        int pos = __builtin_clz(~n_keywords);
+        int pos = std::countr_one(static_cast<unsigned>(n_keywords));
         keywords.push_back(keyword), kwd_ids[pos].push_back(n_keywords);
         automata[pos].add(keyword);
         n_keywords++;
diff --git a/utilities/integer_segments.hpp b/utilities/integer_segments.hpp
index eae49bbf..b2d48a58 100644
--- a/utilities/integer_segments.hpp
+++ b/utilities/integer_segments.hpp
@@ -88,7 +88,7 @@ template <typename Int> struct integer_segments {
     // Count elements strictly less than x, O((# of segments))
     Int order_of_key(Int x) const {
         Int ret = 0;
-        for (auto p : x) {
+        for (auto p : mp) {
             if (p.first >= x) break;
             ret += std::min(x - 1, p.second) - p.first + 1;
         }
diff --git a/utilities/quotients.hpp b/utilities/quotients.hpp
index 79af1519..ede417dd 100644
--- a/utilities/quotients.hpp
+++ b/utilities/quotients.hpp
@@ -8,14 +8,16 @@
 template <class T = long long> std::vector<T> get_quotients(T n) {
     std::vector<T> res;
     for (T x = 1;; ++x) {
-        if (x * x >= n) {
+        const T quotient = n / x;
+        const bool is_square = quotient == x and n % x == 0;
+        if (x > quotient or is_square) {
             const int sz = res.size();
-            if (x * x == n) res.push_back(x);
+            if (is_square) res.push_back(x);
             res.reserve(res.size() + sz);
             for (int i = sz - 1; i >= 0; --i) {
                 T tmp = n / res.at(i);
                 if (tmp < x) continue;
-                if (tmp == x and tmp * tmp == n) continue;
+                if (tmp == x and is_square) continue;
                 res.push_back(tmp);
             }
             return res;
diff --git a/utilities/readline.hpp b/utilities/readline.hpp
index 13aa8ecc..163a7df0 100644
--- a/utilities/readline.hpp
+++ b/utilities/readline.hpp
@@ -12,11 +12,7 @@ std::vector<int> read_ints() {
     if (s.empty()) return {};
     std::stringstream ss(s);
     std::vector<int> ret;
-    while (!ss.eof()) {
-        int t;
-        ss >> t;
-        ret.push_back(t);
-    }
+    for (int t; ss >> t;) ret.push_back(t);
     return ret;
 }