From 4e0235a598bd44db05145b09d43414241c95ff42 Mon Sep 17 00:00:00 2001 From: Not_Leonian <75620009+NotLeonian@users.noreply.github.com> Date: Sat, 11 Jul 2026 16:06:24 +0900 Subject: [PATCH 1/4] =?UTF-8?q?=E3=83=AA=E3=83=95=E3=82=A1=E3=82=AF?= =?UTF-8?q?=E3=82=BF=E3=83=AA=E3=83=B3=E3=82=B0?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../line-convex-polygon-intersection.md | 4 +- docs/graph/others/graph-isomorphism.md | 8 +- docs/graph/tree/hl-rec-dp.md | 4 +- ...determinant-of-linear-matrix-polynomial.md | 6 +- docs/math/matrix/dynamic-matrix-rank.md | 28 +- .../prime-counting-modulo.md | 4 +- .../generalized-floor-sum-degree-le-2.md | 6 +- docs/math/number-theory/generalized-garner.md | 2 +- ...od-range-increment-decrement-operations.md | 6 +- docs/other/rectangle-add-max-get.md | 8 +- docs/structure/others/dynamic-median.md | 8 +- geometry/line-convex-polygon-intersection.hpp | 181 ++++---- graph/others/graph-isomorphism.hpp | 182 ++++---- graph/tree/01-on-tree.hpp | 42 +- graph/tree/hl-rec-dp.hpp | 68 ++- internal/int128.hpp | 103 +++-- math/combinatorics/online-binomial-sum.hpp | 45 +- ...eterminant-of-linear-matrix-polynomial.hpp | 83 ++-- math/matrix/dynamic-matrix-rank.hpp | 404 +++++++++++------- .../prime-counting-modulo.hpp | 84 ++-- math/number-theory/floor-sum.hpp | 85 ++-- .../generalized-floor-sum-degree-le-2.hpp | 131 +++--- math/number-theory/generalized-garner.hpp | 112 +++-- ...rate-maximum-independent-set-path-sums.hpp | 61 +-- ...d-range-increment-decrement-operations.hpp | 15 +- other/rectangle-add-max-get.hpp | 94 ++-- structure/others/dynamic-median.hpp | 40 +- verify/aoj-cgl-4-c.test.cpp | 130 +----- verify/standalone-01-on-tree.test.cpp | 125 ++++++ ...inant-of-linear-matrix-polynomial.test.cpp | 216 ++++++++++ .../standalone-dynamic-matrix-rank.test.cpp | 2 +- verify/standalone-dynamic-median.test.cpp | 49 +++ verify/standalone-floor-sum.test.cpp | 57 +++ verify/standalone-generalized-garner.test.cpp | 30 +- verify/standalone-graph-isomorphism.test.cpp | 14 + verify/standalone-int128.test.cpp | 17 + ...-line-convex-polygon-intersection.test.cpp | 29 ++ ...ge-increment-decrement-operations.test.cpp | 4 + .../standalone-prime-counting-modulo.test.cpp | 11 + ...cal-order-with-minimum-inversions.test.cpp | 39 +- verify/yosupo-sum-of-floor-of-linear.test.cpp | 66 --- 41 files changed, 1638 insertions(+), 965 deletions(-) create mode 100644 verify/standalone-01-on-tree.test.cpp create mode 100644 verify/standalone-determinant-of-linear-matrix-polynomial.test.cpp create mode 100644 verify/standalone-floor-sum.test.cpp diff --git a/docs/geometry/line-convex-polygon-intersection.md b/docs/geometry/line-convex-polygon-intersection.md index d7eaa1d..edb38e8 100644 --- a/docs/geometry/line-convex-polygon-intersection.md +++ b/docs/geometry/line-convex-polygon-intersection.md @@ -68,5 +68,5 @@ documentation_of: geometry/line-convex-polygon-intersection.hpp `hull` のサイズ を $N$ とおく。 -- 時間: $O(\log N)$ -- 空間: $O(1)$ +- 時間 $O(\log N)$ +- 空間 $O(1)$ diff --git a/docs/graph/others/graph-isomorphism.md b/docs/graph/others/graph-isomorphism.md index c5f926e..7dd19b7 100644 --- a/docs/graph/others/graph-isomorphism.md +++ b/docs/graph/others/graph-isomorphism.md @@ -32,8 +32,8 @@ documentation_of: graph/others/graph-isomorphism.hpp また、 $i=1,2$ について、`edges_i` を $E_i$ 、 $M_i=\lvert E_i\rvert$ とする。 - $M_i$ は多重辺を重複込みで数え、自己ループも入力 $1$ 個を $1$ 本として数える。 -- $U$ を $2$ つのグラフの内部辺要素の個数の合計とする。内部辺要素は、同じ頂点対に入った辺を $1$ つにまとめ、本数を `count` として持つ。自己ループの頂点対も $1$ つと数える。 -- 多重辺は重複を消して同一視するのではなく、内部辺要素の `count` まで比較される。 +- $U$ を、各グラフで同じ頂点対に入った辺を $1$ つと数えたときの、 $2$ つのグラフでの個数の合計とする。自己ループの頂点対も $1$ つと数える。 +- 多重辺は重複を消して同一視するのではなく、各頂点対に入った辺の本数まで比較される。 - $D$ を内部隣接リストの全長の合計とする。 $D\le 2U$ である。 - $L$ を色分割中に作られる signature の最大長とする。 $L\le 2n+2$ である。 - $S$ を DFS で訪問した状態数、 $R$ を全状態を通した色分割更新回数の合計とする。 $R\le S(n+1)$ である。 @@ -41,10 +41,8 @@ documentation_of: graph/others/graph-isomorphism.hpp - $C_{\mathrm{build}}=M_1\log M_1+M_2\log M_2+n\log n+n+M_1+M_2$ とする。 - $C_{\mathrm{ref}}=D\log n+Ln\log n+Ln+D+n$ とする。 - $C_{\mathrm{memo}}=n\log S+n$ とする。 -- $C_{\mathrm{check}}=U\log U+U+n$ とする。 - 構築は時間 $O(C_{\mathrm{build}})$ である。 - 色分割更新 $1$ 回は時間 $O(C_{\mathrm{ref}})$ である。 - 探索メモの検索または挿入は時間 $O(C_{\mathrm{memo}})$ である。 -- 色がすべて単独になった状態での辺集合比較は時間 $O(C_{\mathrm{check}})$ である。 -- 全体は時間 $O(C_{\mathrm{build}}+RC_{\mathrm{ref}}+S(C_{\mathrm{memo}}+C_{\mathrm{check}}))$ である。 +- 全体は時間 $O(C_{\mathrm{build}}+RC_{\mathrm{ref}}+SC_{\mathrm{memo}})$ である。 - 最悪の場合、 $S$ は指数的に大きくなる。 diff --git a/docs/graph/tree/hl-rec-dp.md b/docs/graph/tree/hl-rec-dp.md index cbaa2a6..d4b9a86 100644 --- a/docs/graph/tree/hl-rec-dp.md +++ b/docs/graph/tree/hl-rec-dp.md @@ -133,8 +133,8 @@ int main() { ## 計算量 -- 木の根付き化と重い子の選択: $O(n)$ 。 -- 軽い辺だけで再帰が深くなるため、再帰段数: $O(\log n)$ 。 +- 木の根付き化と重い子の選択は時間 $O(n)$ 。 +- 軽い辺だけで再帰が深くなるため、再帰段数は $O(\log n)$ 。 - $K$ を `Spec::K` の値とする。部分木 1 回分の値計算の呼び出し回数は $O(n^{\log_2(K+1)})$ 。 - 全頂点の `before_vertex` と `after_vertex` を回収する実行については、同様に $K\ge 2$ ならば $O(n^{\log_2(K+1)})$ 回、 $K=1$ ならば $O(n\log n)$ 回である。 - 各呼び出しで、`Spec` 側の処理時間が掛かる。典型的な場合として $K=2$ 、各処理が $O(X)$ ならば $O(n^{\log_2 3}X)$ である。 diff --git a/docs/math/matrix/determinant-of-linear-matrix-polynomial.md b/docs/math/matrix/determinant-of-linear-matrix-polynomial.md index 9d5b5ee..cc006fa 100644 --- a/docs/math/matrix/determinant-of-linear-matrix-polynomial.md +++ b/docs/math/matrix/determinant-of-linear-matrix-polynomial.md @@ -30,6 +30,6 @@ documentation_of: math/matrix/determinant-of-linear-matrix-polynomial.hpp ## 計算量 -- `hessenberg_reduction`: $O(N^3)$ -- `characteristic_polynomial`: $O(N^3)$ -- `determinant_of_linear_matrix_polynomial`: $O(N^3)$ +- `hessenberg_reduction`: 時間 $O(N^3)$ +- `characteristic_polynomial`: 時間 $O(N^3)$ +- `determinant_of_linear_matrix_polynomial`: 時間 $O(N^3)$ diff --git a/docs/math/matrix/dynamic-matrix-rank.md b/docs/math/matrix/dynamic-matrix-rank.md index afabc1d..f46f137 100644 --- a/docs/math/matrix/dynamic-matrix-rank.md +++ b/docs/math/matrix/dynamic-matrix-rank.md @@ -21,9 +21,10 @@ documentation_of: math/matrix/dynamic-matrix-rank.hpp - 空に構築する。後で `build(matrix)` を呼ぶ。 - `DynamicMatrixRank(const std::vector>& matrix)` - `matrix` で構築する。 - - 前提: `matrix` は長方形であり、`T` は除算ができる。 + - 前提: `matrix` は長方形であり、`T` は体を成す。 - `void build(const std::vector>& matrix)` - `matrix` を現在の行列として前処理し直す。 + - 前提: `matrix` は長方形であり、`T` は体を成す。 - `void build()` - 現在保持している行列から前処理し直す。 - `int rank() const` @@ -50,23 +51,26 @@ documentation_of: math/matrix/dynamic-matrix-rank.hpp - 備考: 内部状態は変更しない。 - `int apply_rank_one_update(const std::vector& column_vector, const std::vector& row_vector)` - $A+uv^{\top}$ に内部状態を更新し、その階数を返す。 + - 前提: `column_vector` の長さは行数、`row_vector` の長さは列数。 - `int apply_row_replacement(int row_index, const std::vector& new_row)` - `row_index` 行目を `new_row` に差し替え、変更後の階数を返す。 + - 前提: $0\le i& new_column)` - `column_index` 列目を `new_column` に差し替え、変更後の階数を返す。 + - 前提: $0\le j0$ 。 - 前提: 内部計算が `Internal` の範囲に収まる。 - 備考: `T` が符号付き整数型の場合、返り値が `T` の範囲に収まる必要がある。 @@ -39,5 +39,5 @@ documentation_of: math/number-theory/generalized-floor-sum-degree-le-2.hpp ## 計算量 -- 時間: $O(\log m)$ -- 空間: $O(\log m)$ （再帰） +- 時間 $O(\log m)$ +- 空間 $O(\log m)$ （再帰） diff --git a/docs/math/number-theory/generalized-garner.md b/docs/math/number-theory/generalized-garner.md index 8dde6bd..d1dce7a 100644 --- a/docs/math/number-theory/generalized-garner.md +++ b/docs/math/number-theory/generalized-garner.md @@ -29,7 +29,7 @@ documentation_of: math/number-theory/generalized-garner.hpp ## 計算量 -- $N=\lvert a\rvert,\;V=\max_i M_i$ とする。 +- 入力する式の個数を $N$ 、全ての法と $2$ の最大値を $V$ とする。 - 標準の 64 bit 以下の整数型では時間 $O(N\log V)$ 。 - それ以外の整数型で通常乗算を使えない場合は時間 $O(N\log^2 V)$ 。 - 空間 $O(1)$ 。 diff --git a/docs/other/minimum-mod-range-increment-decrement-operations.md b/docs/other/minimum-mod-range-increment-decrement-operations.md index aa3b90f..f50844f 100644 --- a/docs/other/minimum-mod-range-increment-decrement-operations.md +++ b/docs/other/minimum-mod-range-increment-decrement-operations.md @@ -12,7 +12,7 @@ documentation_of: other/minimum-mod-range-increment-decrement-operations.hpp ## 使い方 -- `T minimum_mod_range_increment_decrement_operations(std::vector a, std::vector b, T m)` +- `T minimum_mod_range_increment_decrement_operations(const std::vector &a, const std::vector &b, T m)` - $a$ を $b$ に一致させるための最小操作回数を返す。 - 前提: `a[i]` を $a_i$ 、`b[i]` を $b_i$ として、 $\lvert a\rvert=\lvert b\rvert,\,m>0,\,0\le a_i,b_i Calc to_calc_y(const Point &point) { return static_cast(get_y(point)); } -template -Calc diff_x(const Point &lhs, const Point &rhs) { - return to_calc_x(lhs) - to_calc_x(rhs); -} +template struct LineParameters { + Calc base_x; + Calc base_y; + Calc direction_x; + Calc direction_y; +}; template -Calc diff_y(const Point &lhs, const Point &rhs) { - return to_calc_y(lhs) - to_calc_y(rhs); +LineParameters make_line_parameters(const Point &line_a, + const Point &line_b) { + const Calc base_x = to_calc_x(line_a); + const Calc base_y = to_calc_y(line_a); + return {base_x, base_y, to_calc_x(line_b) - base_x, + to_calc_y(line_b) - base_y}; } template -Calc cross_diff(const Point &lhs_a, const Point &lhs_b, const Point &rhs_a, - const Point &rhs_b) { - return diff_x(lhs_a, lhs_b) * - diff_y(rhs_a, rhs_b) - - diff_y(lhs_a, lhs_b) * - diff_x(rhs_a, rhs_b); +Calc line_height(const Point &point, const LineParameters &line) { + return (to_calc_x(point) - line.base_x) * line.direction_y - + (to_calc_y(point) - line.base_y) * line.direction_x; } template @@ -216,12 +219,14 @@ template int last_non_negative(int length, Func &&func) { return left; } -inline int positive_mod(int index, int mod) { - int ret = index % mod; - if (ret < 0) { - ret += mod; +inline int wrap_nearby_index(int index, int size) { + if (index < 0) { + return index + size; } - return ret; + if (index >= size) { + return index - size; + } + return index; } inline int distance_forward(int from, int to, int n) { @@ -237,7 +242,6 @@ template int compare_value(const T &lhs, const T &rhs) { template std::pair find_extreme_vertices(int n, Func &&height) { - assert(n >= 3); using Value = std::remove_cvref_t; int minimum_index = -1; @@ -288,13 +292,13 @@ std::pair find_extreme_vertices(int n, Func &&height) { const int offset = first_non_negative(n, [&](int index) -> Value { return height(base + index) - height(base + index + 1); }); - maximum_index = positive_mod(minimum_index + offset, n); + maximum_index = wrap_nearby_index(minimum_index + offset, n); } else { const int base = maximum_index - n; const int offset = first_non_negative(n, [&](int index) -> Value { return height(base + index + 1) - height(base + index); }); - minimum_index = positive_mod(maximum_index + offset, n); + minimum_index = wrap_nearby_index(maximum_index + offset, n); } return {minimum_index, maximum_index}; @@ -347,8 +351,16 @@ template struct LineConvexHullIntersectionPoint { return Point{static_cast(x_numerator), static_cast(y_numerator)}; } else { - return Point{static_cast(x_as()), - static_cast(y_as())}; + assert(denominator > 0); + const long double real_denominator = + lpi_internal::number_as_real(denominator); + return Point{ + static_cast( + lpi_internal::number_as_real(x_numerator) / + real_denominator), + static_cast( + lpi_internal::number_as_real(y_numerator) / + real_denominator)}; } } @@ -401,21 +413,20 @@ LineConvexHullIntersectionPoint make_integral_point(const Point &point) { template LineConvexHullIntersectionPoint -line_edge_intersection_integral(const Point &line_a, const Point &line_b, +line_edge_intersection_integral(const LineParameters &line, const Point &segment_a, const Point &segment_b) { - const Calc edge_x = diff_x(segment_b, segment_a); - const Calc edge_y = diff_y(segment_b, segment_a); + const Calc segment_a_x = to_calc_x(segment_a); + const Calc segment_a_y = to_calc_y(segment_a); + const Calc edge_x = to_calc_x(segment_b) - segment_a_x; + const Calc edge_y = to_calc_y(segment_b) - segment_a_y; - Calc denominator = - cross_diff(segment_b, segment_a, line_b, line_a); - assert(denominator != 0); - Calc numerator = cross_diff(line_a, segment_a, line_b, line_a); + Calc denominator = edge_x * line.direction_y - edge_y * line.direction_x; + Calc numerator = (line.base_x - segment_a_x) * line.direction_y - + (line.base_y - segment_a_y) * line.direction_x; - Calc x_numerator = - to_calc_x(segment_a) * denominator + edge_x * numerator; - Calc y_numerator = - to_calc_y(segment_a) * denominator + edge_y * numerator; + Calc x_numerator = segment_a_x * denominator + edge_x * numerator; + Calc y_numerator = segment_a_y * denominator + edge_y * numerator; if (denominator < 0) { denominator = -denominator; @@ -424,32 +435,30 @@ line_edge_intersection_integral(const Point &line_a, const Point &line_b, } const Calc g = gcd3_value(x_numerator, y_numerator, denominator); - if (g != 0) { - x_numerator /= g; - y_numerator /= g; - denominator /= g; - } + x_numerator /= g; + y_numerator /= g; + denominator /= g; return {x_numerator, y_numerator, denominator}; } template -Point line_edge_intersection_floating(const Point &line_a, const Point &line_b, +Point line_edge_intersection_floating(const LineParameters &line, const Point &segment_a, const Point &segment_b) { - const Calc edge_x = diff_x(segment_b, segment_a); - const Calc edge_y = diff_y(segment_b, segment_a); + const Calc segment_a_x = to_calc_x(segment_a); + const Calc segment_a_y = to_calc_y(segment_a); + const Calc edge_x = to_calc_x(segment_b) - segment_a_x; + const Calc edge_y = to_calc_y(segment_b) - segment_a_y; const Calc denominator = - cross_diff(segment_b, segment_a, line_b, line_a); - assert(sign_value(denominator) != 0); - - const Calc numerator = - cross_diff(line_a, segment_a, line_b, line_a); + edge_x * line.direction_y - edge_y * line.direction_x; + const Calc numerator = (line.base_x - segment_a_x) * line.direction_y - + (line.base_y - segment_a_y) * line.direction_x; const Calc t = numerator / denominator; - const Calc x = to_calc_x(segment_a) + edge_x * t; - const Calc y = to_calc_y(segment_a) + edge_y * t; + const Calc x = segment_a_x + edge_x * t; + const Calc y = segment_a_y + edge_y * t; return make_point(x, y); } @@ -464,15 +473,15 @@ result_value_t make_vertex_result(const Point &point) { } template -result_value_t -make_edge_result(const Point &line_a, const Point &line_b, - const Point &segment_a, const Point &segment_b) { +result_value_t make_edge_result(const LineParameters &line, + const Point &segment_a, + const Point &segment_b) { if constexpr (is_integer_v>) { - return line_edge_intersection_integral( - line_a, line_b, segment_a, segment_b); + return line_edge_intersection_integral(line, segment_a, + segment_b); } else { - return line_edge_intersection_floating( - line_a, line_b, segment_a, segment_b); + return line_edge_intersection_floating(line, segment_a, + segment_b); } } @@ -484,16 +493,10 @@ bool equals_point(const Point &lhs, const Point &rhs) { to_calc_y(rhs)); } -template -Calc line_height(const Point &point, const Point &line_a, const Point &line_b) { - return cross_diff(point, line_a, line_b, line_a); -} - template std::vector> line_degenerate_convex_hull_intersection(const std::vector &hull, - const Point &line_a, - const Point &line_b) { + const LineParameters &line) { std::vector> result; if (hull.empty()) { @@ -501,23 +504,21 @@ line_degenerate_convex_hull_intersection(const std::vector &hull, } if (hull.size() == 1) { - if (sign_value(line_height(hull[0], line_a, line_b)) == - 0) { + if (sign_value(line_height(hull[0], line)) == 0) { result.push_back(make_vertex_result(hull[0])); } return result; } - assert(hull.size() == 2); assert((!equals_point(hull[0], hull[1]))); - const Calc first_height = line_height(hull[0], line_a, line_b); - const Calc second_height = - line_height(hull[1], line_a, line_b); + const Calc first_height = line_height(hull[0], line); + const Calc second_height = line_height(hull[1], line); const int first_sign = sign_value(first_height); const int second_sign = sign_value(second_height); if (first_sign == 0 && second_sign == 0) { + result.reserve(2); result.push_back(make_vertex_result(hull[0])); result.push_back(make_vertex_result(hull[1])); return result; @@ -535,8 +536,7 @@ line_degenerate_convex_hull_intersection(const std::vector &hull, if ((first_sign < 0 && second_sign > 0) || (first_sign > 0 && second_sign < 0)) { - result.push_back( - make_edge_result(line_a, line_b, hull[0], hull[1])); + result.push_back(make_edge_result(line, hull[0], hull[1])); } return result; @@ -545,24 +545,19 @@ line_degenerate_convex_hull_intersection(const std::vector &hull, template std::vector> line_strict_convex_polygon_intersection(const std::vector &polygon, - const Point &line_a, - const Point &line_b) { + const LineParameters &line) { const int n = static_cast(polygon.size()); - assert(n >= 3); auto vertex = [&](int index) -> const Point & { - return polygon[positive_mod(index, n)]; + return polygon[wrap_nearby_index(index, n)]; }; auto height = [&](int index) -> Calc { - return line_height(vertex(index), line_a, line_b); + return line_height(vertex(index), line); }; auto chain_index = [&](int start, int step, int offset) -> int { - assert(0 <= start && start < n); - assert(step == 1 || step == -1); - assert(0 <= offset && offset <= n); - return line_convex_polygon_intersection_internal::positive_mod( + return line_convex_polygon_intersection_internal::wrap_nearby_index( start + step * offset, n); }; @@ -570,18 +565,21 @@ line_strict_convex_polygon_intersection(const std::vector &polygon, find_extreme_vertices(n, height); const Calc minimum_value = height(minimum_index); const Calc maximum_value = height(maximum_index); + const int minimum_sign = sign_value(minimum_value); + const int maximum_sign = sign_value(maximum_value); LinePolygonIntersectionResult result; - if (sign_value(minimum_value) > 0 || sign_value(maximum_value) < 0) { + if (minimum_sign > 0 || maximum_sign < 0) { return result; } + result.reserve(2); const int forward_length = distance_forward(minimum_index, maximum_index, n); const int backward_length = n - forward_length; - if (sign_value(minimum_value) == 0) { + if (minimum_sign == 0) { const int forward_zero = last_non_positive(forward_length, [&](int offset) { return height(minimum_index + offset); @@ -604,7 +602,7 @@ line_strict_convex_polygon_intersection(const std::vector &polygon, return result; } - if (sign_value(maximum_value) == 0) { + if (maximum_sign == 0) { const int forward_length_from_max = distance_forward(maximum_index, minimum_index, n); const int backward_length_from_max = n - forward_length_from_max; @@ -635,14 +633,15 @@ line_strict_convex_polygon_intersection(const std::vector &polygon, return height(minimum_index + offset); }); const int first_cross_index = chain_index(minimum_index, +1, first_cross); + const Calc first_cross_height = height(first_cross_index); - if (sign_value(height(first_cross_index)) == 0) { + if (sign_value(first_cross_height) == 0) { result.push_back( make_vertex_result(vertex(first_cross_index))); } else { const int prev_index = chain_index(minimum_index, +1, first_cross - 1); result.push_back(make_edge_result( - line_a, line_b, vertex(prev_index), vertex(first_cross_index))); + line, vertex(prev_index), vertex(first_cross_index))); } const int second_cross = @@ -650,14 +649,15 @@ line_strict_convex_polygon_intersection(const std::vector &polygon, return height(minimum_index - offset); }); const int second_cross_index = chain_index(minimum_index, -1, second_cross); + const Calc second_cross_height = height(second_cross_index); - if (sign_value(height(second_cross_index)) == 0) { + if (sign_value(second_cross_height) == 0) { result.push_back( make_vertex_result(vertex(second_cross_index))); } else { const int prev_index = chain_index(minimum_index, -1, second_cross - 1); result.push_back(make_edge_result( - line_a, line_b, vertex(prev_index), vertex(second_cross_index))); + line, vertex(prev_index), vertex(second_cross_index))); } return result; @@ -682,16 +682,19 @@ line_convex_hull_intersection(const std::vector &hull, lpi_internal::is_signed_v, "integer calculation type must be signed"); - assert((!lpi_internal::equals_point(line_a, line_b))); + const auto line = + lpi_internal::make_line_parameters(line_a, line_b); + assert(lpi_internal::sign_value(line.direction_x) != 0 || + lpi_internal::sign_value(line.direction_y) != 0); if (hull.size() <= 2) { return lpi_internal::line_degenerate_convex_hull_intersection( - hull, line_a, line_b); + hull, line); } return lpi_internal::line_strict_convex_polygon_intersection( - hull, line_a, line_b); + hull, line); } template diff --git a/graph/others/graph-isomorphism.hpp b/graph/others/graph-isomorphism.hpp index b8d17cf..7954d74 100644 --- a/graph/others/graph-isomorphism.hpp +++ b/graph/others/graph-isomorphism.hpp @@ -22,7 +22,6 @@ struct GraphIsomorphism { std::vector>> adjacency; std::vector loop_count; std::vector degree; - std::vector> edges; int edge_count; Graph(int n_, const std::vector> &input_edges) @@ -46,7 +45,6 @@ struct GraphIsomorphism { } std::sort(sorted_edges.begin(), sorted_edges.end()); - edges.reserve(sorted_edges.size()); for (int i = 0; i < static_cast(sorted_edges.size());) { int j = i + 1; while (j < static_cast(sorted_edges.size()) && @@ -56,7 +54,6 @@ struct GraphIsomorphism { const int u = sorted_edges[i][0]; const int v = sorted_edges[i][1]; const int count = j - i; - edges.push_back({u, v, count}); if (u == v) { adjacency[u].push_back({v, count}); loop_count[u] += count; @@ -90,9 +87,7 @@ struct GraphIsomorphism { GraphIsomorphism(int n_, const std::vector> &edges_1, const std::vector> &edges_2) - : n(n_), graph{Graph(n_, edges_1), Graph(n_, edges_2)} { - assert(n >= 0); - } + : n(n_), graph{Graph(n_, edges_1), Graph(n_, edges_2)} {} bool run() { if (graph[0].edge_count != graph[1].edge_count) { @@ -102,26 +97,25 @@ struct GraphIsomorphism { return true; } - std::vector> keys; + std::vector> keys; keys.reserve(2 * n); for (int t = 0; t < 2; ++t) { for (int v = 0; v < n; ++v) { - keys.push_back({graph[t].loop_count[v], graph[t].degree[v]}); + keys.push_back( + {graph[t].loop_count[v], graph[t].degree[v], t, v}); } } std::sort(keys.begin(), keys.end()); - keys.erase(std::unique(keys.begin(), keys.end()), keys.end()); std::array, 2> color{std::vector(n), std::vector(n)}; - for (int t = 0; t < 2; ++t) { - for (int v = 0; v < n; ++v) { - const std::array key{graph[t].loop_count[v], - graph[t].degree[v]}; - color[t][v] = static_cast( - std::lower_bound(keys.begin(), keys.end(), key) - - keys.begin()); + int color_count = 0; + for (int i = 0; i < static_cast(keys.size()); ++i) { + if (i == 0 || keys[i][0] != keys[i - 1][0] || + keys[i][1] != keys[i - 1][1]) { + ++color_count; } + color[keys[i][2]][keys[i][3]] = color_count - 1; } if (!same_color_count(color)) { return false; @@ -132,63 +126,97 @@ struct GraphIsomorphism { } private: - bool refine(std::array, 2> &color) const { + struct RefineWorkspace { + std::vector> signatures; + std::vector multiplicity; + std::vector touched_colors; + std::vector order; + std::vector> count; + std::array, 2> next_color; + }; + + RefineWorkspace refine_workspace; + + bool refine(std::array, 2> &color) { + auto &signatures = refine_workspace.signatures; + auto &multiplicity = refine_workspace.multiplicity; + auto &touched_colors = refine_workspace.touched_colors; + auto &order = refine_workspace.order; + auto &count = refine_workspace.count; + auto &next_color = refine_workspace.next_color; + signatures.resize(2 * n); + order.resize(2 * n); + next_color[0].resize(n); + next_color[1].resize(n); while (true) { - std::vector> signatures(2 * n); + int color_count = 0; + for (int v = 0; v < n; ++v) { + color_count = std::max(color_count, color[0][v] + 1); + } + + multiplicity.assign(color_count, 0); + touched_colors.clear(); + touched_colors.reserve(color_count); for (int t = 0; t < 2; ++t) { for (int v = 0; v < n; ++v) { - std::vector> neighbor_colors; - neighbor_colors.reserve(graph[t].adjacency[v].size()); for (const auto &[to, count] : graph[t].adjacency[v]) { if (to != v) { - neighbor_colors.push_back({color[t][to], count}); + const int neighbor_color = color[t][to]; + if (multiplicity[neighbor_color] == 0) { + touched_colors.push_back(neighbor_color); + } + multiplicity[neighbor_color] += count; } } - std::sort(neighbor_colors.begin(), neighbor_colors.end()); + std::sort(touched_colors.begin(), touched_colors.end()); - std::vector signature; - signature.reserve(2 + 2 * neighbor_colors.size()); + std::vector &signature = signatures[t * n + v]; + signature.clear(); + signature.reserve(2 + 2 * touched_colors.size()); signature.push_back(color[t][v]); signature.push_back(graph[t].loop_count[v]); - for (int i = 0; - i < static_cast(neighbor_colors.size());) { - int j = i + 1; - int count_sum = neighbor_colors[i].second; - while (j < static_cast(neighbor_colors.size()) && - neighbor_colors[i].first == - neighbor_colors[j].first) { - count_sum += neighbor_colors[j].second; - ++j; - } - signature.push_back(neighbor_colors[i].first); - signature.push_back(count_sum); - i = j; + for (int neighbor_color : touched_colors) { + signature.push_back(neighbor_color); + signature.push_back(multiplicity[neighbor_color]); + multiplicity[neighbor_color] = 0; } - signatures[t * n + v] = std::move(signature); + touched_colors.clear(); } } - std::vector> keys = signatures; - std::sort(keys.begin(), keys.end()); - keys.erase(std::unique(keys.begin(), keys.end()), keys.end()); + for (int i = 0; i < 2 * n; ++i) { + order[i] = i; + } + std::sort(order.begin(), order.end(), [&](int lhs, int rhs) { + return signatures[lhs] < signatures[rhs]; + }); + + int next_color_count = 0; + for (int i = 0; i < 2 * n; ++i) { + if (i == 0 || + signatures[order[i - 1]] != signatures[order[i]]) { + ++next_color_count; + } + const int t = order[i] / n; + const int v = order[i] % n; + next_color[t][v] = next_color_count - 1; + } - std::array, 2> next_color{std::vector(n), - std::vector(n)}; + count.assign(next_color_count, std::array{}); for (int t = 0; t < 2; ++t) { for (int v = 0; v < n; ++v) { - next_color[t][v] = static_cast( - std::lower_bound(keys.begin(), keys.end(), - signatures[t * n + v]) - - keys.begin()); + ++count[next_color[t][v]][t]; } } - if (!same_color_count(next_color)) { - return false; + for (const auto &color_count_pair : count) { + if (color_count_pair[0] != color_count_pair[1]) { + return false; + } } - if (next_color[0] == color[0] && next_color[1] == color[1]) { + if (next_color_count == color_count) { return true; } - color = std::move(next_color); + color.swap(next_color); } } @@ -197,11 +225,6 @@ struct GraphIsomorphism { return false; } - StateKey key = make_state_key(color); - if (dead_states.find(key) != dead_states.end()) { - return false; - } - int color_count = 0; for (int v = 0; v < n; ++v) { color_count = std::max(color_count, color[0][v] + 1); @@ -222,7 +245,12 @@ struct GraphIsomorphism { } if (branch_color == -1) { - return check_mapping(color); + return true; + } + + StateKey key = make_state_key(color); + if (dead_states.contains(key)) { + return false; } int u = -1; @@ -232,7 +260,6 @@ struct GraphIsomorphism { break; } } - assert(u != -1); const int new_color = color_count; for (int v = 0; v < n; ++v) { @@ -256,42 +283,17 @@ struct GraphIsomorphism { color_count = std::max(color_count, color[t][v] + 1); } } - std::vector count0(color_count, 0), count1(color_count, 0); - for (int v = 0; v < n; ++v) { - ++count0[color[0][v]]; - ++count1[color[1][v]]; - } - return count0 == count1; - } - - bool check_mapping(const std::array, 2> &color) const { - int color_count = 0; - for (int v = 0; v < n; ++v) { - color_count = std::max(color_count, color[0][v] + 1); - } - std::vector position(color_count, -1); + std::vector> count(color_count); for (int v = 0; v < n; ++v) { - position[color[1][v]] = v; + ++count[color[0][v]][0]; + ++count[color[1][v]][1]; } - - std::vector permutation(n, -1); - for (int v = 0; v < n; ++v) { - permutation[v] = position[color[0][v]]; - assert(permutation[v] != -1); - } - - std::vector> mapped_edges; - mapped_edges.reserve(graph[0].edges.size()); - for (const auto &edge : graph[0].edges) { - int u = permutation[edge[0]]; - int v = permutation[edge[1]]; - if (v < u) { - std::swap(u, v); + for (const auto &color_count_pair : count) { + if (color_count_pair[0] != color_count_pair[1]) { + return false; } - mapped_edges.push_back({u, v, edge[2]}); } - std::sort(mapped_edges.begin(), mapped_edges.end()); - return mapped_edges == graph[1].edges; + return true; } static std::uint64_t mix(std::uint64_t x) { diff --git a/graph/tree/01-on-tree.hpp b/graph/tree/01-on-tree.hpp index 79b5363..5477583 100644 --- a/graph/tree/01-on-tree.hpp +++ b/graph/tree/01-on-tree.hpp @@ -1,5 +1,5 @@ -#ifndef GRAPH_TREE_01_ON_TREE_HPP -#define GRAPH_TREE_01_ON_TREE_HPP +#ifndef GRAPH_TREE_ZERO_ONE_ON_TREE_HPP +#define GRAPH_TREE_ZERO_ONE_ON_TREE_HPP // 根付き木上で、親が子より左に出る順序の 01 列の転倒数最小値を求める。 // 頂点 v には c0[v] 個の 0 の後に c1[v] 個の 1 を置いた列が書かれている。 @@ -86,8 +86,10 @@ solve_01_on_tree(int n, const std::vector> &edges, std::vector parent(n, -2); parent[root] = -1; - int visited = 0; - std::vector stack{root}; + [[maybe_unused]] int visited = 0; + std::vector stack; + stack.reserve(n); + stack.push_back(root); while (!stack.empty()) { const int v = stack.back(); stack.pop_back(); @@ -112,6 +114,9 @@ solve_01_on_tree(int n, const std::vector> &edges, struct QueueCompare { bool operator()(const QueueNode &a, const QueueNode &b) const { + if (a.zero == b.zero && a.one == b.one) { + return a.vertex < b.vertex; + } if (a.one == Count{} && b.one == Count{}) { return a.vertex < b.vertex; } @@ -121,6 +126,15 @@ solve_01_on_tree(int n, const std::vector> &edges, if (b.one == Count{}) { return true; } + if (a.zero == Count{} && b.zero == Count{}) { + return a.vertex < b.vertex; + } + if (a.zero == Count{}) { + return true; + } + if (b.zero == Count{}) { + return false; + } const int result = zero_one_on_tree_impl::compare_fraction( a.zero, a.one, b.zero, b.one); if (result < 0) { @@ -133,7 +147,7 @@ solve_01_on_tree(int n, const std::vector> &edges, } }; - std::vector leader(n), up = parent, version(n, 0); + std::vector leader(n), up = std::move(parent), version(n, 0); for (int v = 0; v < n; ++v) { leader[v] = v; } @@ -146,7 +160,10 @@ solve_01_on_tree(int n, const std::vector> &edges, return v; }; - std::priority_queue, QueueCompare> que; + std::vector queue_storage; + queue_storage.reserve(n); + std::priority_queue, QueueCompare> que( + QueueCompare{}, std::move(queue_storage)); auto push = [&](int v) { if (up[v] != -1) { que.push(QueueNode{v, version[v], zero[v], one[v]}); @@ -159,20 +176,19 @@ solve_01_on_tree(int n, const std::vector> &edges, Answer answer{}; for (int merge_count = 0; merge_count < n - 1; ++merge_count) { - QueueNode current{}; + int v = -1; while (true) { - assert(!que.empty()); - current = que.top(); + const QueueNode ¤t = que.top(); + const int current_vertex = current.vertex; + const int current_version = current.version; que.pop(); - const int v = find(current.vertex); - if (v == current.vertex && current.version == version[v]) { + v = find(current_vertex); + if (v == current_vertex && current_version == version[v]) { break; } } - const int v = current.vertex; const int p = find(up[v]); - assert(p != v); answer += static_cast(one[p]) * static_cast(zero[v]); zero[p] += zero[v]; diff --git a/graph/tree/hl-rec-dp.hpp b/graph/tree/hl-rec-dp.hpp index e9c22a3..806aaf3 100644 --- a/graph/tree/hl-rec-dp.hpp +++ b/graph/tree/hl-rec-dp.hpp @@ -41,31 +41,6 @@ inline void put_heavy_child_first(std::vector> &child, } } -inline std::vector -validate_rooted_tree(const std::vector> &child, int root) { - const int n = static_cast(child.size()); - assert(n >= 1); - assert(0 <= root && root < n); - std::vector parent(n, -2); - std::vector order; - order.reserve(n); - std::vector stack{root}; - parent[root] = -1; - while (!stack.empty()) { - const int v = stack.back(); - stack.pop_back(); - order.push_back(v); - for (int to : child[v]) { - assert(0 <= to && to < n); - assert(parent[to] == -2); - parent[to] = v; - stack.push_back(to); - } - } - assert(static_cast(order.size()) == n); - return order; -} - inline std::vector> make_rooted_tree(int n, const std::vector> &edges, int root) { @@ -81,29 +56,33 @@ make_rooted_tree(int n, const std::vector> &edges, graph[v].push_back(u); } - std::vector> child(n); std::vector parent(n, -2); std::vector order; order.reserve(n); - std::vector stack{root}; + std::vector stack; + stack.reserve(n); + stack.push_back(root); parent[root] = -1; while (!stack.empty()) { const int v = stack.back(); stack.pop_back(); order.push_back(v); + int child_count = 0; for (int to : graph[v]) { if (to == parent[v]) { continue; } assert(parent[to] == -2); parent[to] = v; - child[v].push_back(to); + graph[v][child_count] = to; + ++child_count; stack.push_back(to); } + graph[v].resize(child_count); } assert(static_cast(order.size()) == n); - put_heavy_child_first(child, order); - return child; + put_heavy_child_first(graph, order); + return graph; } template struct hl_rec_dp_runner { @@ -113,17 +92,20 @@ template struct hl_rec_dp_runner { using Pack = std::array; std::vector> child; + std::vector path; Spec &spec; hl_rec_dp_runner(std::vector> child, Spec &spec) - : child(std::move(child)), spec(spec) {} + : child(std::move(child)), spec(spec) { + path.reserve(this->child.size()); + } Pack run(int root, const State &initial_state) { - return dfs_heavy_path(root, initial_state, true); + return dfs_heavy_path(root, initial_state); } - Pack dfs_heavy_path(int start, const State &in, bool collect) { - std::vector path; + template Pack dfs_heavy_path(int start, const State &in) { + const int path_begin = static_cast(path.size()); for (int v = start;; v = child[v][0]) { path.push_back(v); if (child[v].empty()) { @@ -131,31 +113,33 @@ template struct hl_rec_dp_runner { } } - Pack cur = spec.make_pack(path.back(), in); - for (int i = static_cast(path.size()) - 1; i >= 0; --i) { + const int path_end = static_cast(path.size()); + Pack cur = spec.make_pack(path[path_end - 1], in); + for (int i = path_end - 1; i >= path_begin; --i) { const int v = path[i]; - if (i + 1 < static_cast(path.size())) { + if (i + 1 < path_end) { cur = spec.take_heavy(v, path[i + 1], std::move(cur)); } for (int idx = 1; idx < static_cast(child[v].size()); ++idx) { const int to = child[v][idx]; for (int lane = 0; lane < K; ++lane) { - Pack got = dfs_heavy_path(to, cur[lane], false); - cur[lane] = spec.take_light(v, to, lane, std::move(got)); + cur[lane] = spec.take_light( + v, to, lane, dfs_heavy_path(to, cur[lane])); } } - if (collect) { + if constexpr (Collect) { spec.before_vertex(v, cur); } spec.add_vertex(v, cur); - if (collect) { + if constexpr (Collect) { spec.after_vertex(v, cur); for (int idx = 1; idx < static_cast(child[v].size()); ++idx) { - dfs_heavy_path(child[v][idx], in, true); + dfs_heavy_path(child[v][idx], in); } } } + path.resize(path_begin); return cur; } }; diff --git a/internal/int128.hpp b/internal/int128.hpp index f42950c..c90e898 100644 --- a/internal/int128.hpp +++ b/internal/int128.hpp @@ -1,5 +1,5 @@ -#ifndef INTERNAL_INT128_HPP -#define INTERNAL_INT128_HPP +#ifndef INTERNAL_INT_ONE_TWO_EIGHT_HPP +#define INTERNAL_INT_ONE_TWO_EIGHT_HPP // 128 bit 符号なし整数型 UInt128 と 128 bit 符号付き整数型 Int128 を提供する。 // GCC 拡張には依存しない。 @@ -7,8 +7,9 @@ #include #include #include -#include +#include #include +#include #include #include @@ -48,6 +49,9 @@ class UInt128 { explicit constexpr operator bool() const { return high_ != 0 || low_ != 0; } explicit constexpr operator long double() const { + if (high_ == 0) { + return static_cast(low_); + } return static_cast(high_) * 0x1p64L + static_cast(low_); } @@ -131,26 +135,7 @@ class UInt128 { static constexpr void div_mod(UInt128 lhs, UInt128 rhs, UInt128 "ient, UInt128 &remainder) { assert(rhs != UInt128{}); - - UInt128 q, r; - if (rhs.high_ == 0) { - if (lhs.high_ == 0) { - q = UInt128(lhs.low_ / rhs.low_); - r = UInt128(lhs.low_ % rhs.low_); - } else if (rhs.low_ == 1) { - q = lhs; - r = UInt128{}; - } else { - div_mod_32bit_words(lhs, rhs, q, r); - } - } else if (lhs < rhs) { - q = UInt128{}; - r = lhs; - } else { - div_mod_32bit_words(lhs, rhs, q, r); - } - quotient = q; - remainder = r; + div_mod_unchecked(lhs, rhs, quotient, remainder); } friend constexpr UInt128 operator/(UInt128 lhs, UInt128 rhs) { @@ -168,6 +153,8 @@ class UInt128 { } private: + friend class Int128; + static constexpr std::uint64_t word_base() { return std::uint64_t{1} << 32; } @@ -195,11 +182,10 @@ class UInt128 { return length; } + // divisor != 0 が呼び出し元で保証されていることを仮定する。 static constexpr UInt128 div_mod_uint32(UInt128 value, std::uint32_t divisor, std::uint32_t &remainder) { - assert(divisor != 0); - const std::uint32_t words[4] = { static_cast(value.low_ & word_mask()), static_cast(value.low_ >> 32), @@ -351,7 +337,7 @@ class UInt128 { } } - // 呼び出し元で lhs >= rhs が保証されていることを仮定する。 + // 呼び出し元で lhs >= rhs > 0 が保証されていることを仮定する。 static constexpr void div_mod_32bit_words(UInt128 lhs, UInt128 rhs, UInt128 "ient, UInt128 &remainder) { @@ -362,7 +348,6 @@ class UInt128 { const int lhs_length = word_length(lhs_words); const int rhs_length = word_length(rhs_words); - assert(rhs_length > 0); if (rhs_length == 1) { std::uint32_t rem = 0; @@ -380,6 +365,36 @@ class UInt128 { remainder = from_32bit_words(remainder_words); } + // rhs != 0 が呼び出し元で保証されていることを仮定する。 + static constexpr void div_mod_unchecked(UInt128 lhs, UInt128 rhs, + UInt128 "ient, + UInt128 &remainder) { + UInt128 q, r; + if (rhs.high_ == 0) { + if (lhs.high_ == 0) { + q = UInt128(lhs.low_ / rhs.low_); + r = UInt128(lhs.low_ % rhs.low_); + } else if (rhs.low_ == 1) { + q = lhs; + r = UInt128{}; + } else if (rhs.low_ <= std::numeric_limits::max()) { + std::uint32_t rem = 0; + q = div_mod_uint32(lhs, static_cast(rhs.low_), + rem); + r = UInt128(rem); + } else { + div_mod_32bit_words(lhs, rhs, q, r); + } + } else if (lhs < rhs) { + q = UInt128{}; + r = lhs; + } else { + div_mod_32bit_words(lhs, rhs, q, r); + } + quotient = q; + remainder = r; + } + static constexpr UInt128 multiply_u64(std::uint64_t lhs, std::uint64_t rhs) { constexpr std::uint64_t mask = (std::uint64_t{1} << 32) - 1; @@ -428,9 +443,10 @@ class Int128 { constexpr bool is_negative() const { return (value_.high() >> 63) != 0; } explicit constexpr operator long double() const { + const bool negative = is_negative(); const long double magnitude = - static_cast(abs_unsigned(*this)); - return is_negative() ? -magnitude : magnitude; + static_cast(negative ? -value_ : value_); + return negative ? -magnitude : magnitude; } friend constexpr bool operator==(Int128 lhs, Int128 rhs) { @@ -515,16 +531,18 @@ class Int128 { assert(rhs != Int128{}); assert(!(lhs == min_value() && rhs == Int128(-1))); - const bool quotient_negative = lhs.is_negative() != rhs.is_negative(); - const bool remainder_negative = lhs.is_negative(); + const bool lhs_negative = lhs.is_negative(); + const bool rhs_negative = rhs.is_negative(); + const bool quotient_negative = lhs_negative != rhs_negative; UInt128 quotient_abs; UInt128 remainder_abs; - UInt128::div_mod(abs_unsigned(lhs), abs_unsigned(rhs), quotient_abs, - remainder_abs); + UInt128::div_mod_unchecked(abs_unsigned(lhs), abs_unsigned(rhs), + quotient_abs, remainder_abs); - const Int128 q = from_unsigned(quotient_abs, quotient_negative); - const Int128 r = from_unsigned(remainder_abs, remainder_negative); + const Int128 q = + from_unsigned_unchecked(quotient_abs, quotient_negative); + const Int128 r = from_unsigned_unchecked(remainder_abs, lhs_negative); quotient = q; remainder = r; @@ -565,10 +583,16 @@ class Int128 { static constexpr Int128 from_unsigned(UInt128 value, bool negative) { if (!negative) { assert((value.high() >> 63) == 0); - return from_twos_complement(value); + } else { + assert(value <= UInt128::from_words(std::uint64_t{1} << 63, 0)); } - assert(value <= UInt128::from_words(std::uint64_t{1} << 63, 0)); - return from_twos_complement(-value); + return from_unsigned_unchecked(value, negative); + } + + // 符号を付けた value が Int128 の範囲に収まることを仮定する。 + static constexpr Int128 from_unsigned_unchecked(UInt128 value, + bool negative) { + return from_twos_complement(negative ? -value : value); } friend std::istream &operator>>(std::istream &input, Int128 &value); @@ -578,10 +602,9 @@ class Int128 { }; namespace int128_internal { +// first < text.size() が呼び出し元で保証されていることを仮定する。 inline UInt128 read_uint128_decimal(const std::string &text, std::size_t first) { - assert(first < text.size()); - UInt128 value = 0; std::size_t i = first; diff --git a/math/combinatorics/online-binomial-sum.hpp b/math/combinatorics/online-binomial-sum.hpp index 0792ea7..c9a0483 100644 --- a/math/combinatorics/online-binomial-sum.hpp +++ b/math/combinatorics/online-binomial-sum.hpp @@ -119,6 +119,21 @@ template struct OnlineBinomialSum { assert(m >= 0); assert(m <= max_m); + return binom_prefix_sum_unchecked(n, m); + } + + T binom_sum(int l, int u, int m) const { + assert(l >= 0); + assert(l <= u); + assert(m >= 0); + assert(m <= max_m); + + return binom_prefix_sum_unchecked(u, m) - + binom_prefix_sum_unchecked(l, m); + } + + private: + T binom_prefix_sum_unchecked(int n, int m) const { if (n == 0) { return T(); } @@ -132,6 +147,9 @@ template struct OnlineBinomialSum { if (m == 0) { return T(1); } + if (n > m) { + return T(); + } T ans = binomial(m - 1, n - 1); if ((n - 1) % 2 == 1) { @@ -149,18 +167,24 @@ template struct OnlineBinomialSum { sample_sum_table[sample_m_index * sample_n_count + sample_n_index]; while (current_n < n) { - sum += power_r[current_n] * binomial(current_m, current_n); + if (current_n <= current_m) { + sum += power_r[current_n] * binomial(current_m, current_n); + } ++current_n; } while (current_n > n) { --current_n; - sum -= power_r[current_n] * binomial(current_m, current_n); + if (current_n <= current_m) { + sum -= power_r[current_n] * binomial(current_m, current_n); + } } while (current_m < m) { sum *= r_plus_one; - sum -= power_r[current_n] * binomial(current_m, current_n - 1); + if (current_n - 1 <= current_m) { + sum -= power_r[current_n] * binomial(current_m, current_n - 1); + } ++current_m; } @@ -173,19 +197,7 @@ template struct OnlineBinomialSum { return sum; } - T binom_sum(int l, int u, int m) const { - assert(l >= 0); - assert(l <= u); - - return binom_prefix_sum(u, m) - binom_prefix_sum(l, m); - } - - private: T binomial(int n, int k) const { - if (k < 0 || k > n) { - return T(); - } - return factorial[n] * inverse_factorial[k] * inverse_factorial[n - k]; } @@ -207,9 +219,6 @@ template struct OnlineBinomialSum { int value) const { int index = value / bucket_size; const int sample_count = static_cast(sample_list.size()); - if (index >= sample_count) { - index = sample_count - 1; - } if (index + 1 >= sample_count) { return index; } diff --git a/math/matrix/determinant-of-linear-matrix-polynomial.hpp b/math/matrix/determinant-of-linear-matrix-polynomial.hpp index 3b30266..55fc2ed 100644 --- a/math/matrix/determinant-of-linear-matrix-polynomial.hpp +++ b/math/matrix/determinant-of-linear-matrix-polynomial.hpp @@ -12,14 +12,21 @@ #include #include +namespace determinant_of_linear_matrix_polynomial_internal { template -void hessenberg_reduction(std::vector> &matrix) { +bool is_square_matrix(const std::vector> &matrix) { const int n = static_cast(matrix.size()); - assert(n == 0 || static_cast(matrix[0].size()) == n); - for (int i = 1; i < n; ++i) { - assert(static_cast(matrix[i].size()) == n); + for (const std::vector &row : matrix) { + if (static_cast(row.size()) != n) { + return false; + } } + return true; +} +template +void hessenberg_reduction(std::vector> &matrix) { + const int n = static_cast(matrix.size()); for (int r = 0; r < n - 2; ++r) { int piv = -1; for (int h = r + 1; h < n; ++h) { @@ -45,7 +52,8 @@ void hessenberg_reduction(std::vector> &matrix) { if (coef == T()) { continue; } - for (int j = 0; j < n; ++j) { + matrix[i][r] = T(); + for (int j = r + 1; j < n; ++j) { matrix[i][j] -= matrix[r + 1][j] * coef; } for (int j = 0; j < n; ++j) { @@ -58,12 +66,8 @@ void hessenberg_reduction(std::vector> &matrix) { template std::vector characteristic_polynomial(std::vector> matrix) { const int n = static_cast(matrix.size()); - assert(n == 0 || static_cast(matrix[0].size()) == n); - for (int i = 1; i < n; ++i) { - assert(static_cast(matrix[i].size()) == n); - } - - hessenberg_reduction(matrix); + determinant_of_linear_matrix_polynomial_internal::hessenberg_reduction( + matrix); // p[i] = det(x I_i - matrix[0..i-1][0..i-1])（係数は昇順） std::vector> p(n + 1); @@ -71,16 +75,21 @@ std::vector characteristic_polynomial(std::vector> matrix) { for (int i = 0; i < n; ++i) { p[i + 1].assign(i + 2, T()); - for (int j = 0; j <= i; ++j) { - p[i + 1][j + 1] += p[i][j]; - } - for (int j = 0; j <= i; ++j) { - p[i + 1][j] -= p[i][j] * matrix[i][i]; + p[i + 1][0] = T() - p[i][0] * matrix[i][i]; + for (int j = 1; j <= i; ++j) { + p[i + 1][j] = p[i][j - 1] - p[i][j] * matrix[i][i]; } + p[i + 1][i + 1] = p[i][i]; T betas = T(1); for (int j = i - 1; j >= 0; --j) { betas *= matrix[j + 1][j]; + if (betas == T()) { + break; + } + if (matrix[j][i] == T()) { + continue; + } const T hb = (T() - matrix[j][i]) * betas; for (int k = 0; k <= j; ++k) { p[i + 1][k] += hb * p[j][k]; @@ -89,6 +98,25 @@ std::vector characteristic_polynomial(std::vector> matrix) { } return p[n]; } +} // namespace determinant_of_linear_matrix_polynomial_internal + +template +void hessenberg_reduction(std::vector> &matrix) { + assert(determinant_of_linear_matrix_polynomial_internal::is_square_matrix( + matrix)); + + determinant_of_linear_matrix_polynomial_internal::hessenberg_reduction( + matrix); +} + +template +std::vector characteristic_polynomial(std::vector> matrix) { + assert(determinant_of_linear_matrix_polynomial_internal::is_square_matrix( + matrix)); + + return determinant_of_linear_matrix_polynomial_internal:: + characteristic_polynomial(std::move(matrix)); +} template std::vector @@ -96,13 +124,13 @@ determinant_of_linear_matrix_polynomial(std::vector> M0, std::vector> M1) { const int n = static_cast(M0.size()); assert(static_cast(M1.size()) == n); + assert( + determinant_of_linear_matrix_polynomial_internal::is_square_matrix(M0)); + assert( + determinant_of_linear_matrix_polynomial_internal::is_square_matrix(M1)); if (n == 0) { return {T(1)}; } - for (int i = 0; i < n; ++i) { - assert(static_cast(M0[i].size()) == n); - assert(static_cast(M1[i].size()) == n); - } int multiply_by_x = 0; // 特定の列に x を掛ける操作の回数 T det_inv = T(1); // 1 / (det A det B) @@ -134,9 +162,10 @@ determinant_of_linear_matrix_polynomial(std::vector> M0, } } - // (M0 + x M1) の p 列に x を掛ける（M1 の p 列が 0 のとき swap で実現できる）。 + // M1 の p 列が 0 なので、M0 の p 列を M1 に移して列に x を掛ける。 for (int i = 0; i < n; ++i) { - std::swap(M0[i][p], M1[i][p]); + M1[i][p] = std::move(M0[i][p]); + M0[i][p] = T(); } --p; // 同じ列をやり直す（高々 n 回） @@ -150,11 +179,13 @@ determinant_of_linear_matrix_polynomial(std::vector> M0, } const T v = M1[p][p]; - assert(v != T()); det_inv *= v; const T vinv = T(1) / v; for (int col = 0; col < n; ++col) { M0[p][col] *= vinv; + } + M1[p][p] = T(1); + for (int col = p + 1; col < n; ++col) { M1[p][col] *= vinv; } @@ -168,6 +199,9 @@ determinant_of_linear_matrix_polynomial(std::vector> M0, } for (int col = 0; col < n; ++col) { M0[row][col] -= M0[p][col] * coef; + } + M1[row][p] = T(); + for (int col = p + 1; col < n; ++col) { M1[row][col] -= M1[p][col] * coef; } } @@ -179,7 +213,8 @@ determinant_of_linear_matrix_polynomial(std::vector> M0, M0[i][j] = T() - M0[i][j]; } } - std::vector poly = characteristic_polynomial(M0); + std::vector poly = determinant_of_linear_matrix_polynomial_internal:: + characteristic_polynomial(std::move(M0)); for (T &c : poly) { c *= det_inv; } diff --git a/math/matrix/dynamic-matrix-rank.hpp b/math/matrix/dynamic-matrix-rank.hpp index 48939d6..7e84d22 100644 --- a/math/matrix/dynamic-matrix-rank.hpp +++ b/math/matrix/dynamic-matrix-rank.hpp @@ -5,7 +5,8 @@ // さらに、内部状態を O((k + 1)(r + c)) で更新しつつ変更後の階数を返せる。 // 現在の行列は左右の階数分解と片側逆元で保持する。 // 1 行差し替え、1 列差し替えはそれぞれ e_i (b-a_i)^T, (b-a_j) e_j^T に帰着する。 -// 前処理は O(rc min(r, c) + k^2(r + c))、更新や判定は O((k + 1)(r + c)) である。 +// 前処理は O(r + rc min(r, c) + k^2(r + c))、更新や判定は +// O((k + 1)(r + c)) である。 // T は体を成し、零判定と四則演算ができることを仮定する。 #include @@ -28,17 +29,20 @@ template struct DynamicMatrixRank { } void build(const std::vector> &matrix) { + assert(is_rectangular(matrix)); + + build_unchecked(matrix); + } + + private: + void build_unchecked(const std::vector> &matrix) { row_size = static_cast(matrix.size()); column_size = row_size == 0 ? 0 : static_cast(matrix[0].size()); - for (int i = 1; i < row_size; ++i) { - assert(static_cast(matrix[i].size()) == column_size); - } - - const std::vector basis_columns = find_independent_columns(matrix); - const std::vector basis_rows = - find_independent_columns(transpose_matrix(matrix)); + const IndependentSubmatrix independent_submatrix = + find_independent_submatrix(matrix); + const std::vector &basis_rows = independent_submatrix.rows; + const std::vector &basis_columns = independent_submatrix.columns; matrix_rank = static_cast(basis_columns.size()); - assert(static_cast(basis_rows.size()) == matrix_rank); std::vector> intersection_matrix( matrix_rank, std::vector(matrix_rank, T())); @@ -87,7 +91,8 @@ template struct DynamicMatrixRank { } } - void build() { build(materialize_matrix()); } + public: + void build() { build_unchecked(materialize_matrix()); } int rank() const { return matrix_rank; } @@ -109,10 +114,13 @@ template struct DynamicMatrixRank { std::vector get_column(int column_index) const { assert(0 <= column_index && column_index < column_size); std::vector column(row_size, T()); - for (int i = 0; i < row_size; ++i) { - for (int j = 0; j < matrix_rank; ++j) { - column[i] += - column_space_basis[i][j] * row_space_basis[j][column_index]; + for (int j = 0; j < matrix_rank; ++j) { + const T value = row_space_basis[j][column_index]; + if (value == T()) { + continue; + } + for (int i = 0; i < row_size; ++i) { + column[i] += column_space_basis[i][j] * value; } } return column; @@ -137,39 +145,45 @@ template struct DynamicMatrixRank { int rank_after_rank_one_update(const std::vector &column_vector, const std::vector &row_vector) const { - return analyze_rank_one_update(column_vector, row_vector).next_rank; + assert(static_cast(column_vector.size()) == row_size); + assert(static_cast(row_vector.size()) == column_size); + + return rank_after_rank_one_update_unchecked(column_vector, row_vector); } int rank_after_row_replacement(int row_index, const std::vector &new_row) const { assert(0 <= row_index && row_index < row_size); assert(static_cast(new_row.size()) == column_size); - std::vector difference = new_row; - const std::vector current_row = get_row(row_index); - for (int j = 0; j < column_size; ++j) { - difference[j] -= current_row[j]; - } + const std::vector difference = + row_replacement_difference(row_index, new_row); std::vector column_vector(row_size, T()); column_vector[row_index] = T(1); - return rank_after_rank_one_update(column_vector, difference); + return rank_after_rank_one_update_unchecked(column_vector, difference); } int rank_after_column_replacement(int column_index, const std::vector &new_column) const { assert(0 <= column_index && column_index < column_size); assert(static_cast(new_column.size()) == row_size); - std::vector difference = new_column; - const std::vector current_column = get_column(column_index); - for (int i = 0; i < row_size; ++i) { - difference[i] -= current_column[i]; - } + const std::vector difference = + column_replacement_difference(column_index, new_column); std::vector row_vector(column_size, T()); row_vector[column_index] = T(1); - return rank_after_rank_one_update(difference, row_vector); + return rank_after_rank_one_update_unchecked(difference, row_vector); } int apply_rank_one_update(const std::vector &column_vector, const std::vector &row_vector) { + assert(static_cast(column_vector.size()) == row_size); + assert(static_cast(row_vector.size()) == column_size); + + return apply_rank_one_update_unchecked(column_vector, row_vector); + } + + private: + int apply_rank_one_update_unchecked(const std::vector &column_vector, + const std::vector &row_vector) { RankOneUpdateInfo info = analyze_rank_one_update(column_vector, row_vector); @@ -181,34 +195,31 @@ template struct DynamicMatrixRank { const std::vector rho = normalized_right_annihilator( pivot_column, info.row_residual[pivot_column]); - const std::vector> old_row_space_right_inverse = - row_space_right_inverse; - append_column(column_space_basis, info.column_residual); column_space_left_inverse.push_back(lambda); row_space_basis.push_back(row_vector); for (int i = 0; i < matrix_rank; ++i) { + if (info.alpha[i] == T()) { + continue; + } for (int j = 0; j < column_size; ++j) { row_space_basis[i][j] += info.alpha[i] * row_vector[j]; } } - std::vector right_alpha(column_size, T()); - for (int i = 0; i < column_size; ++i) { - for (int j = 0; j < matrix_rank; ++j) { - right_alpha[i] += - old_row_space_right_inverse[i][j] * info.alpha[j]; - } - } + const std::vector right_alpha = + multiply_row_space_right_inverse(info.alpha); - row_space_right_inverse.assign( - column_size, std::vector(matrix_rank + 1, T())); for (int i = 0; i < column_size; ++i) { + row_space_right_inverse[i].resize(matrix_rank + 1, T()); + if (rho[i] == T()) { + row_space_right_inverse[i][matrix_rank] = + T() - right_alpha[i]; + continue; + } for (int j = 0; j < matrix_rank; ++j) { - row_space_right_inverse[i][j] = - old_row_space_right_inverse[i][j] - - rho[i] * info.beta[j]; + row_space_right_inverse[i][j] -= rho[i] * info.beta[j]; } row_space_right_inverse[i][matrix_rank] = T() - right_alpha[i] + rho[i] * info.schur; @@ -222,18 +233,22 @@ template struct DynamicMatrixRank { const std::vector lambda = normalized_left_annihilator( pivot_row, info.column_residual[pivot_row]); for (int i = 0; i < row_size; ++i) { + if (column_vector[i] == T()) { + continue; + } for (int j = 0; j < matrix_rank; ++j) { column_space_basis[i][j] += column_vector[i] * info.beta[j]; } } - std::vector> new_left_inverse = - column_space_left_inverse; for (int i = 0; i < matrix_rank; ++i) { + if (info.alpha[i] == T()) { + continue; + } for (int j = 0; j < row_size; ++j) { - new_left_inverse[i][j] -= info.alpha[i] * lambda[j]; + column_space_left_inverse[i][j] -= + info.alpha[i] * lambda[j]; } } - column_space_left_inverse.swap(new_left_inverse); return matrix_rank; } @@ -242,38 +257,46 @@ template struct DynamicMatrixRank { const std::vector rho = normalized_right_annihilator( pivot_column, info.row_residual[pivot_column]); for (int i = 0; i < matrix_rank; ++i) { + if (info.alpha[i] == T()) { + continue; + } for (int j = 0; j < column_size; ++j) { row_space_basis[i][j] += info.alpha[i] * row_vector[j]; } } - std::vector> new_right_inverse = - row_space_right_inverse; for (int i = 0; i < column_size; ++i) { + if (rho[i] == T()) { + continue; + } for (int j = 0; j < matrix_rank; ++j) { - new_right_inverse[i][j] -= rho[i] * info.beta[j]; + row_space_right_inverse[i][j] -= rho[i] * info.beta[j]; } } - row_space_right_inverse.swap(new_right_inverse); return matrix_rank; } if (info.schur != T()) { - std::vector right_alpha(column_size, T()); - for (int i = 0; i < column_size; ++i) { - for (int j = 0; j < matrix_rank; ++j) { - right_alpha[i] += - row_space_right_inverse[i][j] * info.alpha[j]; - } + const T schur_inverse = T(1) / info.schur; + for (T &value : info.beta) { + value *= schur_inverse; } + const std::vector right_alpha = + multiply_row_space_right_inverse(info.alpha); for (int i = 0; i < matrix_rank; ++i) { + if (info.alpha[i] == T()) { + continue; + } for (int j = 0; j < column_size; ++j) { row_space_basis[i][j] += info.alpha[i] * row_vector[j]; } } for (int i = 0; i < column_size; ++i) { + if (right_alpha[i] == T()) { + continue; + } for (int j = 0; j < matrix_rank; ++j) { row_space_right_inverse[i][j] -= - right_alpha[i] * (info.beta[j] / info.schur); + right_alpha[i] * info.beta[j]; } } return matrix_rank; @@ -282,12 +305,8 @@ template struct DynamicMatrixRank { const int removed = first_nonzero(info.alpha); const T removed_inverse = T(1) / info.alpha[removed]; - std::vector right_alpha(column_size, T()); - for (int i = 0; i < column_size; ++i) { - for (int j = 0; j < matrix_rank; ++j) { - right_alpha[i] += row_space_right_inverse[i][j] * info.alpha[j]; - } - } + const std::vector right_alpha = + multiply_row_space_right_inverse(info.alpha); std::vector> new_column_space_basis( row_size, std::vector(matrix_rank - 1, T())); @@ -335,34 +354,34 @@ template struct DynamicMatrixRank { return matrix_rank; } + public: int apply_row_replacement(int row_index, const std::vector &new_row) { assert(0 <= row_index && row_index < row_size); assert(static_cast(new_row.size()) == column_size); - std::vector difference = new_row; - const std::vector current_row = get_row(row_index); - for (int j = 0; j < column_size; ++j) { - difference[j] -= current_row[j]; - } + const std::vector difference = + row_replacement_difference(row_index, new_row); std::vector column_vector(row_size, T()); column_vector[row_index] = T(1); - return apply_rank_one_update(column_vector, difference); + return apply_rank_one_update_unchecked(column_vector, difference); } int apply_column_replacement(int column_index, const std::vector &new_column) { assert(0 <= column_index && column_index < column_size); assert(static_cast(new_column.size()) == row_size); - std::vector difference = new_column; - const std::vector current_column = get_column(column_index); - for (int i = 0; i < row_size; ++i) { - difference[i] -= current_column[i]; - } + const std::vector difference = + column_replacement_difference(column_index, new_column); std::vector row_vector(column_size, T()); row_vector[column_index] = T(1); - return apply_rank_one_update(difference, row_vector); + return apply_rank_one_update_unchecked(difference, row_vector); } private: + struct IndependentSubmatrix { + std::vector rows; + std::vector columns; + }; + struct RankOneUpdateInfo { std::vector alpha; std::vector beta; @@ -374,37 +393,67 @@ template struct DynamicMatrixRank { int next_rank = 0; }; + static bool is_rectangular(const std::vector> &matrix) { + if (matrix.empty()) { + return true; + } + const int column_count = static_cast(matrix[0].size()); + for (const std::vector &row : matrix) { + if (static_cast(row.size()) != column_count) { + return false; + } + } + return true; + } + RankOneUpdateInfo analyze_rank_one_update(const std::vector &column_vector, const std::vector &row_vector) const { - assert(static_cast(column_vector.size()) == row_size); - assert(static_cast(row_vector.size()) == column_size); - RankOneUpdateInfo info; info.alpha = multiply_left_inverse(column_vector); info.beta = multiply_right_inverse(row_vector); - const std::vector projected_column = - reconstruct_column_vector(info.alpha); - const std::vector projected_row = reconstruct_row_vector(info.beta); + info.column_residual = column_vector; + for (int j = 0; j < matrix_rank; ++j) { + const T value = info.alpha[j]; + if (value == T()) { + continue; + } + for (int i = 0; i < row_size; ++i) { + info.column_residual[i] -= column_space_basis[i][j] * value; + } + } + info.row_residual = row_vector; + for (int i = 0; i < matrix_rank; ++i) { + const T value = info.beta[i]; + if (value == T()) { + continue; + } + for (int j = 0; j < column_size; ++j) { + info.row_residual[j] -= value * row_space_basis[i][j]; + } + } - info.column_residual.resize(row_size, T()); - info.row_residual.resize(column_size, T()); info.column_inside = true; info.row_inside = true; for (int i = 0; i < row_size; ++i) { - info.column_residual[i] = column_vector[i] - projected_column[i]; if (info.column_residual[i] != T()) { info.column_inside = false; + break; } } for (int j = 0; j < column_size; ++j) { - info.row_residual[j] = row_vector[j] - projected_row[j]; if (info.row_residual[j] != T()) { info.row_inside = false; + break; } } + if (info.column_inside != info.row_inside) { + info.next_rank = matrix_rank; + return info; + } + info.schur = T(1); for (int i = 0; i < matrix_rank; ++i) { info.schur += info.beta[i] * info.alpha[i]; @@ -412,61 +461,128 @@ template struct DynamicMatrixRank { if (!info.column_inside && !info.row_inside) { info.next_rank = matrix_rank + 1; - } else if (info.column_inside != info.row_inside) { - info.next_rank = matrix_rank; } else { info.next_rank = info.schur == T() ? matrix_rank - 1 : matrix_rank; } return info; } + int rank_after_rank_one_update_unchecked( + const std::vector &column_vector, + const std::vector &row_vector) const { + const std::vector alpha = multiply_left_inverse(column_vector); + const std::vector beta = multiply_right_inverse(row_vector); + + bool column_inside = true; + for (int i = 0; i < row_size; ++i) { + T residual = column_vector[i]; + for (int j = 0; j < matrix_rank; ++j) { + residual -= column_space_basis[i][j] * alpha[j]; + } + if (residual != T()) { + column_inside = false; + break; + } + } + + bool row_inside = true; + for (int j = 0; j < column_size; ++j) { + T residual = row_vector[j]; + for (int i = 0; i < matrix_rank; ++i) { + residual -= beta[i] * row_space_basis[i][j]; + } + if (residual != T()) { + row_inside = false; + break; + } + } + + if (!column_inside && !row_inside) { + return matrix_rank + 1; + } + if (column_inside != row_inside) { + return matrix_rank; + } + + T schur = T(1); + for (int i = 0; i < matrix_rank; ++i) { + schur += beta[i] * alpha[i]; + } + return schur == T() ? matrix_rank - 1 : matrix_rank; + } + std::vector - multiply_left_inverse(const std::vector &column_vector) const { - std::vector result(matrix_rank, T()); + row_replacement_difference(int row_index, + const std::vector &new_row) const { + std::vector difference = new_row; for (int i = 0; i < matrix_rank; ++i) { - for (int j = 0; j < row_size; ++j) { - result[i] += column_space_left_inverse[i][j] * column_vector[j]; + const T value = column_space_basis[row_index][i]; + if (value == T()) { + continue; + } + for (int j = 0; j < column_size; ++j) { + difference[j] -= value * row_space_basis[i][j]; } } - return result; + return difference; } std::vector - multiply_right_inverse(const std::vector &row_vector) const { + column_replacement_difference(int column_index, + const std::vector &new_column) const { + std::vector difference = new_column; + for (int j = 0; j < matrix_rank; ++j) { + const T value = row_space_basis[j][column_index]; + if (value == T()) { + continue; + } + for (int i = 0; i < row_size; ++i) { + difference[i] -= column_space_basis[i][j] * value; + } + } + return difference; + } + + std::vector + multiply_left_inverse(const std::vector &column_vector) const { std::vector result(matrix_rank, T()); - for (int j = 0; j < column_size; ++j) { - const T value = row_vector[j]; + for (int j = 0; j < row_size; ++j) { + const T value = column_vector[j]; if (value == T()) { continue; } for (int i = 0; i < matrix_rank; ++i) { - result[i] += value * row_space_right_inverse[j][i]; + result[i] += column_space_left_inverse[i][j] * value; } } return result; } std::vector - reconstruct_column_vector(const std::vector &coefficients) const { - std::vector result(row_size, T()); - for (int i = 0; i < row_size; ++i) { - for (int j = 0; j < matrix_rank; ++j) { - result[i] += column_space_basis[i][j] * coefficients[j]; + multiply_right_inverse(const std::vector &row_vector) const { + std::vector result(matrix_rank, T()); + for (int j = 0; j < column_size; ++j) { + const T value = row_vector[j]; + if (value == T()) { + continue; + } + for (int i = 0; i < matrix_rank; ++i) { + result[i] += value * row_space_right_inverse[j][i]; } } return result; } std::vector - reconstruct_row_vector(const std::vector &coefficients) const { + multiply_row_space_right_inverse(const std::vector &coefficients) const { std::vector result(column_size, T()); - for (int i = 0; i < matrix_rank; ++i) { - const T value = coefficients[i]; + for (int j = 0; j < matrix_rank; ++j) { + const T value = coefficients[j]; if (value == T()) { continue; } - for (int j = 0; j < column_size; ++j) { - result[j] += value * row_space_basis[i][j]; + for (int i = 0; i < column_size; ++i) { + result[i] += row_space_right_inverse[i][j] * value; } } return result; @@ -485,8 +601,9 @@ template struct DynamicMatrixRank { result[i] -= value * column_space_left_inverse[j][i]; } } + const T pivot_inverse = T(1) / pivot_value; for (int i = 0; i < row_size; ++i) { - result[i] /= pivot_value; + result[i] *= pivot_inverse; } return result; } @@ -495,22 +612,24 @@ template struct DynamicMatrixRank { const T &pivot_value) const { std::vector result(column_size, T()); result[pivot_column] = T(1); - for (int i = 0; i < column_size; ++i) { - for (int j = 0; j < matrix_rank; ++j) { - result[i] -= row_space_right_inverse[i][j] * - row_space_basis[j][pivot_column]; + for (int j = 0; j < matrix_rank; ++j) { + const T value = row_space_basis[j][pivot_column]; + if (value == T()) { + continue; + } + for (int i = 0; i < column_size; ++i) { + result[i] -= row_space_right_inverse[i][j] * value; } } + const T pivot_inverse = T(1) / pivot_value; for (int i = 0; i < column_size; ++i) { - result[i] /= pivot_value; + result[i] *= pivot_inverse; } return result; } static void append_column(std::vector> &matrix, const std::vector &column) { - assert(static_cast(matrix.size()) == - static_cast(column.size())); for (int i = 0; i < static_cast(matrix.size()); ++i) { matrix[i].push_back(column[i]); } @@ -522,37 +641,27 @@ template struct DynamicMatrixRank { return i; } } - assert(false); return -1; } - static std::vector> - transpose_matrix(const std::vector> &matrix) { + static IndependentSubmatrix + find_independent_submatrix(const std::vector> &matrix) { const int h = static_cast(matrix.size()); if (h == 0) { return {}; } const int w = static_cast(matrix[0].size()); - std::vector> transposed(w, std::vector(h, T())); + std::vector> b = matrix; + std::vector original_rows(h); for (int i = 0; i < h; ++i) { - assert(static_cast(matrix[i].size()) == w); - for (int j = 0; j < w; ++j) { - transposed[j][i] = matrix[i][j]; - } + original_rows[i] = i; } - return transposed; - } - static std::vector - find_independent_columns(const std::vector> &matrix) { - const int h = static_cast(matrix.size()); - if (h == 0) { - return {}; - } - const int w = static_cast(matrix[0].size()); - std::vector> b = matrix; int rank = 0; - std::vector pivot_columns; + IndependentSubmatrix result; + const int maximum_rank = h < w ? h : w; + result.rows.reserve(maximum_rank); + result.columns.reserve(maximum_rank); for (int column = 0; column < w; ++column) { int pivot = -1; for (int row = rank; row < h; ++row) { @@ -566,10 +675,14 @@ template struct DynamicMatrixRank { } if (pivot != rank) { b[pivot].swap(b[rank]); + const int original_row = original_rows[pivot]; + original_rows[pivot] = original_rows[rank]; + original_rows[rank] = original_row; } const T inverse = T(1) / b[rank][column]; - for (int j = column; j < w; ++j) { + b[rank][column] = T(1); + for (int j = column + 1; j < w; ++j) { b[rank][j] *= inverse; } for (int row = rank + 1; row < h; ++row) { @@ -577,26 +690,24 @@ template struct DynamicMatrixRank { if (factor == T()) { continue; } - for (int j = column; j < w; ++j) { + b[row][column] = T(); + for (int j = column + 1; j < w; ++j) { b[row][j] -= b[rank][j] * factor; } } - pivot_columns.push_back(column); + result.rows.push_back(original_rows[rank]); + result.columns.push_back(column); ++rank; if (rank == h) { break; } } - return pivot_columns; + return result; } static std::vector> inverse_matrix(std::vector> matrix) { const int n = static_cast(matrix.size()); - for (int i = 0; i < n; ++i) { - assert(static_cast(matrix[i].size()) == n); - } - std::vector> inverse(n, std::vector(n, T())); for (int i = 0; i < n; ++i) { inverse[i][i] = T(1); @@ -610,15 +721,17 @@ template struct DynamicMatrixRank { break; } } - assert(pivot >= 0); if (pivot != column) { matrix[pivot].swap(matrix[column]); inverse[pivot].swap(inverse[column]); } const T diagonal_inverse = T(1) / matrix[column][column]; - for (int j = 0; j < n; ++j) { + matrix[column][column] = T(1); + for (int j = column + 1; j < n; ++j) { matrix[column][j] *= diagonal_inverse; + } + for (int j = 0; j < n; ++j) { inverse[column][j] *= diagonal_inverse; } for (int row = 0; row < n; ++row) { @@ -629,8 +742,11 @@ template struct DynamicMatrixRank { if (factor == T()) { continue; } - for (int j = 0; j < n; ++j) { + matrix[row][column] = T(); + for (int j = column + 1; j < n; ++j) { matrix[row][j] -= matrix[column][j] * factor; + } + for (int j = 0; j < n; ++j) { inverse[row][j] -= inverse[column][j] * factor; } } diff --git a/math/multiplicative-function/prime-counting-modulo.hpp b/math/multiplicative-function/prime-counting-modulo.hpp index 92d6498..654896d 100644 --- a/math/multiplicative-function/prime-counting-modulo.hpp +++ b/math/multiplicative-function/prime-counting-modulo.hpp @@ -8,12 +8,12 @@ // 計算量 O(m N^{3/4} / log N)、空間 O(m sqrt(N))。 #include +#include #include #include namespace prime_counting_modulo_internal { inline long long integer_sqrt(long long n) { - assert(n >= 0); long long ok = 0, ng = 1; while (ng <= n / ng) { ng <<= 1; @@ -29,10 +29,8 @@ inline long long integer_sqrt(long long n) { return ok; } -inline long long count_residue_2_to_n(long long n, long long m, long long r) { - assert(n >= 0); - assert(m > 0); - assert(0 <= r && r < m); +inline long long count_residue_2_to_n(long long n, long long m, long long r, + long long one_residue) { long long res = 0; if (r <= n) { res = (n - r) / m + 1; @@ -40,30 +38,22 @@ inline long long count_residue_2_to_n(long long n, long long m, long long r) { if (r == 0) { --res; } - if (n >= 1 && r == 1 % m) { + if (n >= 1 && r == one_residue) { --res; } return res; } -inline void add_mod(long long &x, long long a, long long m) { - assert(0 <= x && x < m); - assert(0 <= a && a < m); - if (a == 0) { - return; - } +inline void add_mod_nonzero(long long &x, long long a, long long m) { if (x >= m - a) { x -= m - a; } else { x += a; } } -} // namespace prime_counting_modulo_internal inline std::pair, std::vector>> -prime_counting_modulo_table(long long N, long long m) { - assert(N >= 0); - assert(m > 0); +make_table(long long N, long long m) { using i64 = long long; std::vector ns{0}; for (i64 i = N; i > 0;) { @@ -77,10 +67,11 @@ prime_counting_modulo_table(long long N, long long m) { const i64 sq = prime_counting_modulo_internal::integer_sqrt(N); const i64 nsz = static_cast(ns.size()); std::vector> h(m, std::vector(nsz)); + const i64 one_residue = m == 1 ? 0 : 1; for (i64 r = 0; r < m; ++r) { for (i64 i = 0; i < nsz; ++i) { h[r][i] = prime_counting_modulo_internal::count_residue_2_to_n( - ns[i], m, r); + ns[i], m, r, one_residue); } } for (i64 x = 2; x <= sq; ++x) { @@ -91,18 +82,41 @@ prime_counting_modulo_table(long long N, long long m) { continue; } const i64 x2 = x * x; + const i64 direct_index_limit = sq / x; for (i64 i = 1; i < nsz && ns[i] >= x2; ++i) { const i64 n = ns[i]; const i64 q = n / x; - const i64 q_idx = i <= sq / x ? i * x : nsz - q; - i64 to = 0; - for (i64 r = 0; r < m; ++r) { - h[to][i] -= h[r][q_idx] - h[r][prev_idx]; - prime_counting_modulo_internal::add_mod(to, x_mod, m); + const i64 q_idx = i <= direct_index_limit ? i * x : nsz - q; + if (x_mod == 0) { + i64 removed = 0; + for (i64 r = 0; r < m; ++r) { + removed += h[r][q_idx] - h[r][prev_idx]; + } + h[0][i] -= removed; + } else if (x_mod == 1) { + for (i64 r = 0; r < m; ++r) { + h[r][i] -= h[r][q_idx] - h[r][prev_idx]; + } + } else { + i64 to = 0; + for (i64 r = 0; r < m; ++r) { + h[to][i] -= h[r][q_idx] - h[r][prev_idx]; + prime_counting_modulo_internal::add_mod_nonzero(to, x_mod, + m); + } } } } - return {ns, h}; + return {std::move(ns), std::move(h)}; +} +} // namespace prime_counting_modulo_internal + +inline std::pair, std::vector>> +prime_counting_modulo_table(long long N, long long m) { + assert(N >= 0); + assert(m > 0); + + return prime_counting_modulo_internal::make_table(N, m); } inline std::vector prime_counting_modulo(long long N, long long m) { @@ -112,7 +126,7 @@ inline std::vector prime_counting_modulo(long long N, long long m) { if (N == 0) { return res; } - const auto table = prime_counting_modulo_table(N, m).second; + const auto table = prime_counting_modulo_internal::make_table(N, m).second; for (long long r = 0; r < m; ++r) { res[r] = table[r][1]; } @@ -124,19 +138,23 @@ std::vector> prime_counting_modulo_mf_prefix_sum_table(long long N, long long m) { assert(N >= 0); assert(m > 0); - std::vector> res(m); if (N == 0) { - return res; + return std::vector>(m); } - const auto table = prime_counting_modulo_table(N, m).second; - for (long long r = 0; r < m; ++r) { - res[r].resize(table[r].size()); - for (long long i = 0; i < static_cast(table[r].size()); - ++i) { - res[r][i] = static_cast(table[r][i]); + auto table = prime_counting_modulo_internal::make_table(N, m).second; + if constexpr (std::is_same_v) { + return table; + } else { + std::vector> res(m); + for (long long r = 0; r < m; ++r) { + res[r].resize(table[r].size()); + for (long long i = 0; i < static_cast(table[r].size()); + ++i) { + res[r][i] = static_cast(table[r][i]); + } } + return res; } - return res; } #endif diff --git a/math/number-theory/floor-sum.hpp b/math/number-theory/floor-sum.hpp index ed7a278..ebd29c3 100644 --- a/math/number-theory/floor-sum.hpp +++ b/math/number-theory/floor-sum.hpp @@ -13,42 +13,34 @@ #include namespace floor_sum_internal { -template T floor_div(T x, T y) { - assert(y > 0); +template std::pair floor_div_mod(T x, T y) { if constexpr (std::numeric_limits::is_signed) { + if (x >= 0 && x < y) { + return {0, x}; + } T q = x / y; - T r = x % y; + T r = x - q * y; if (r < 0) { q -= 1; - } - return q; - } else { - return x / y; - } -} - -template T floor_mod(T x, T y) { - assert(y > 0); - if constexpr (std::numeric_limits::is_signed) { - T r = x % y; - if (r < 0) { r += y; } - return r; + return {q, r}; } else { - return x % y; + if (x < y) { + return {0, x}; + } + const T q = x / y; + return {q, x - q * y}; } } // Σ_{i=0}^{n-1} i = n(n-1)/2 template T sum_0_to_n_minus_1(T n) { - if (n == 0) { - return 0; + const T half = n / 2; + if (n == half * 2) { + return half * (n - 1); } - if (n % 2 == 0) { - return (n / 2) * (n - 1); - } - return n * ((n - 1) / 2); + return n * half; } } // namespace floor_sum_internal @@ -61,40 +53,49 @@ template T floor_sum(T n, T m, T a, T b) { if (n == 0) { return 0; } + if (n == 1) { + return floor_sum_internal::floor_div_mod(b, m).first; + } T ans = 0; { - const T q = floor_sum_internal::floor_div(a, m); - a = floor_sum_internal::floor_mod(a, m); - ans += q * floor_sum_internal::sum_0_to_n_minus_1(n); - } - { - const T q = floor_sum_internal::floor_div(b, m); - b = floor_sum_internal::floor_mod(b, m); - ans += q * n; - } - - while (true) { - if (a >= m) { - const T q = a / m; + const auto [q, r] = floor_sum_internal::floor_div_mod(a, m); + a = r; + if (q != 0) { ans += q * floor_sum_internal::sum_0_to_n_minus_1(n); - a %= m; } - if (b >= m) { - const T q = b / m; + } + { + const auto [q, r] = floor_sum_internal::floor_div_mod(b, m); + b = r; + if (q != 0) { ans += q * n; - b %= m; } + } + if (a == 0) { + return ans; + } + while (true) { const T y_max = a * n + b; if (y_max < m) { break; } - n = y_max / m; - b = y_max % m; + const T new_n = y_max / m; + b = y_max - new_n * m; + n = new_n; std::swap(m, a); + + const T qa = a / m; + ans += qa * floor_sum_internal::sum_0_to_n_minus_1(n); + a -= qa * m; + if (b >= m) { + const T qb = b / m; + ans += qb * n; + b -= qb * m; + } } return ans; } diff --git a/math/number-theory/generalized-floor-sum-degree-le-2.hpp b/math/number-theory/generalized-floor-sum-degree-le-2.hpp index 1abcba9..2aa5773 100644 --- a/math/number-theory/generalized-floor-sum-degree-le-2.hpp +++ b/math/number-theory/generalized-floor-sum-degree-le-2.hpp @@ -1,9 +1,11 @@ -#ifndef MATH_NUMBER_THEORY_GENERALIZED_FLOOR_SUM_DEGREE_LE_2_HPP -#define MATH_NUMBER_THEORY_GENERALIZED_FLOOR_SUM_DEGREE_LE_2_HPP - -// 一般化 floor sum（次数 2 以下）のうち、(0,1),(1,1),(0,2) をまとめて求める。 -// ans_01 = Σ floor((a i + b)/m)、ans_11 = Σ i*floor((a i + b)/m)、ans_02 = Σ floor((a i + b)/m)^2。 -// n >= 0、m > 0 を仮定する。a, b は負でもよい。 +#ifndef MATH_NUMBER_THEORY_GENERALIZED_FLOOR_SUM_DEGREE_LE_TWO_HPP +#define MATH_NUMBER_THEORY_GENERALIZED_FLOOR_SUM_DEGREE_LE_TWO_HPP + +// 一般化 floor sum（次数 2 以下）のうち、（0, 1）、（1, 1）、（0, 2）をまとめて求める。 +// ans_01 = Σ floor((a i + b) / m)、ans_11 = Σ i floor((a i + b) / m)、 +// ans_02 = Σ floor((a i + b) / m)^2。 +// n >= 0、m > 0 を仮定する。T が符号付きの場合、a, b は負でもよい。 +// T および明示的に指定する Internal は整数型であり、内部計算が内部型の範囲に収まることを仮定する。 // 計算量 O(log m)。 #include @@ -39,65 +41,68 @@ using default_internal_t = std::conditional_t, NicheLibrary::Int128, T>; template std::pair floor_div_mod(T x, T y) { - assert(y > 0); if constexpr (is_signed_v) { + if (x >= 0 && x < y) { + return {0, x}; + } T q = x / y; - T r = x % y; + T r = x - q * y; if (r < 0) { q -= T(1); r += y; } return {q, r}; } else { - return {x / y, x % y}; + if (x < y) { + return {0, x}; + } + const T q = x / y; + return {q, x - q * y}; } } // Σ_{i=0}^{n-1} i = n(n-1)/2 template Value sum_0_to_n_minus_1(State n) { - if (n == 0) { - return 0; - } - if (n % 2 == 0) { - return static_cast(n / 2) * static_cast(n - 1); + const State half = n / 2; + if (n == half * 2) { + return static_cast(half) * static_cast(n - 1); } - return static_cast(n) * static_cast((n - 1) / 2); + return static_cast(n) * static_cast(half); } // Σ_{i=0}^{n-1} i^2 = (n-1)n(2n-1)/6 template Value sum_0_to_n_minus_1_sq(State n) { - if (n == 0) { - return 0; - } Value a = static_cast(n - 1); Value b = static_cast(n); Value c = static_cast(2) * static_cast(n) - Value(1); - if (a % Value(2) == 0) { - a /= Value(2); - } else if (b % Value(2) == 0) { + const State n_mod_6 = n % 6; + if (n_mod_6 == 0) { + b /= Value(6); + } else if (n_mod_6 == 1) { + a /= Value(6); + } else if (n_mod_6 == 2) { b /= Value(2); - } else { - c /= Value(2); - } - if (a % Value(3) == 0) { - a /= Value(3); - } else if (b % Value(3) == 0) { + c /= Value(3); + } else if (n_mod_6 == 3) { + a /= Value(2); b /= Value(3); + } else if (n_mod_6 == 4) { + a /= Value(3); + b /= Value(2); } else { + a /= Value(2); c /= Value(3); } return a * b * c; } -template Value sum_range(State l, State r) { - // Σ_{i=l}^{r} i - if (l > r) { - return 0; - } - const State count = r - l + 1; - const Value sum = static_cast(l) + static_cast(r); - if (count % 2 == 0) { - return static_cast(count / 2) * sum; +// Σ_{i=l}^{n-1} i +template Value sum_l_to_n_minus_1(State l, State n) { + const State count = n - l; + const Value sum = static_cast(l) + static_cast(n) - Value(1); + const State half_count = count / 2; + if (count == half_count * 2) { + return static_cast(half_count) * sum; } return static_cast(count) * (sum / Value(2)); } @@ -110,26 +115,11 @@ template struct Result { template Result solve(State n, State m, State a, State b) { - if constexpr (is_signed_v) { - assert(n >= 0); - } - assert(m > 0); - if (n == 0) { - return {0, 0, 0}; - } - const auto [qa_state, a_mod] = floor_div_mod(a, m); const auto [qb_state, b_mod] = floor_div_mod(b, m); a = a_mod; b = b_mod; - if constexpr (is_signed_v) { - assert(a >= 0); - assert(b >= 0); - } - assert(a < m); - assert(b < m); - const Value qa = static_cast(qa_state); const Value qb = static_cast(qb_state); @@ -147,13 +137,15 @@ Result solve(State n, State m, State a, State b) { const Value x_max = y_max_value * static_cast(m) - static_cast(b); - const State t = - static_cast((x_max + static_cast(a) - Value(1)) / - static_cast(a)); - - const State b2 = static_cast( - (static_cast(a) - (x_max % static_cast(a))) % - static_cast(a)); + const Value a_value = static_cast(a); + const Value x_max_div_a = x_max / a_value; + const Value x_max_mod_a = x_max - x_max_div_a * a_value; + const bool has_remainder = x_max_mod_a != Value(0); + const State t = static_cast(x_max_div_a) + + (has_remainder ? State(1) : State(0)); + const State b2 = has_remainder + ? static_cast(a_value - x_max_mod_a) + : State(0); const auto rec = solve(y_max, a, m, b2); @@ -173,7 +165,7 @@ Result solve(State n, State m, State a, State b) { const Value y = static_cast(y_max); const Value tail_01 = tail_len * y; - const Value tail_11 = y * sum_range(t, n - 1); + const Value tail_11 = y * sum_l_to_n_minus_1(t, n); const Value tail_02 = tail_len * y * y; base = {head_01 + tail_01, head_11 + tail_11, head_02 + tail_02}; @@ -181,8 +173,21 @@ Result solve(State n, State m, State a, State b) { } const Value n_value = static_cast(n); + if (qa == Value(0)) { + if (qb == Value(0)) { + return base; + } + const Value si = sum_0_to_n_minus_1(n); + return {qb * n_value + base.ans_01, qb * si + base.ans_11, + qb * qb * n_value + Value(2) * qb * base.ans_01 + base.ans_02}; + } + const Value si = sum_0_to_n_minus_1(n); const Value si2 = sum_0_to_n_minus_1_sq(n); + if (qb == Value(0)) { + return {qa * si + base.ans_01, qa * si2 + base.ans_11, + qa * qa * si2 + Value(2) * qa * base.ans_11 + base.ans_02}; + } Result ans; ans.ans_01 = qa * si + qb * n_value + base.ans_01; @@ -220,6 +225,14 @@ generalized_floor_sum_degree_le_2(T n, T m, T a, T b) { assert(n >= 0); } assert(m > 0); + if (n == 0) { + return {T(0), T(0), T(0)}; + } + if (n == 1) { + const T q = gfs_internal::floor_div_mod(b, m).first; + const Value q_value = static_cast(q); + return {q, T(0), static_cast(q_value * q_value)}; + } const auto res = gfs_internal::solve(n, m, a, b); return {static_cast(res.ans_01), static_cast(res.ans_11), diff --git a/math/number-theory/generalized-garner.hpp b/math/number-theory/generalized-garner.hpp index a8dbdf4..40bd233 100644 --- a/math/number-theory/generalized-garner.hpp +++ b/math/number-theory/generalized-garner.hpp @@ -7,9 +7,11 @@ // 係数や法は互いに素でなくてよい。解がなければ (0, 0) を返す。 // 標準の 64 bit 以下の整数型では剰余乗算に 128 bit 整数型を用いる。 // それ以外の整数型では加算と 2 倍による剰余乗算を用いる。 -// 計算量 O(N log max M_i)。ただし、通常乗算を使えない型では O(N log^2 max M_i)。 +// 式の個数を N、2 と全ての法の最大値を V とすると、計算量 O(N log V)。 +// ただし、通常乗算を使えない型では O(N log^2 V)。 #include +#include #include #include #include @@ -37,7 +39,13 @@ constexpr bool use_uint128_mul_v = std::numeric_limits>::digits <= 64; template T safe_mod(T x, T m) { - assert(m > 0); + if constexpr (is_signed_v) { + if (x >= 0 && x < m) { + return x; + } + } else if (x < m) { + return x; + } x %= m; if constexpr (is_signed_v) { if (x < 0) { @@ -48,8 +56,6 @@ template T safe_mod(T x, T m) { } template T gcd(T a, T b) { - assert(a >= 0); - assert(b >= 0); while (b != 0) { const T c = a % b; a = b; @@ -59,28 +65,37 @@ template T gcd(T a, T b) { } template T add_mod(T a, T b, T m) { - assert(0 <= a && a < m); - assert(0 <= b && b < m); - if (a >= m - b) { - return a - (m - b); + const T gap = m - b; + if (a >= gap) { + return a - gap; } return a + b; } -template T mul_mod(T a, T b, T m) { - assert(m > 0); - a = safe_mod(a, m); - b = safe_mod(b, m); +template T mul_mod_normalized(T a, T b, T m) { + if (a == 0 || b == 0) { + return 0; + } + if (a == 1) { + return b; + } + if (b == 1) { + return a; + } if constexpr (use_uint128_mul_v) { using U = NicheLibrary::UInt128; return static_cast((U(a) * U(b)) % U(m)); } else { + if (a < b) { + std::swap(a, b); + } T res = 0; while (b != 0) { - if (b % 2 != 0) { + const T half_b = b / 2; + if (b != half_b * 2) { res = add_mod(res, a, m); } - b /= 2; + b = half_b; if (b != 0) { a = add_mod(a, a, m); } @@ -90,56 +105,83 @@ template T mul_mod(T a, T b, T m) { } template T inv_mod(T a, T m) { - assert(m > 1); - a = safe_mod(a, m); - assert(gcd(a, m) == 1); - T b = m; - T x = 1; - T y = 0; - while (b != 0) { + if (a == 1) { + return 1; + } + T b = a; + a = m; + T x = 0; + T y = 1; + while (b != 1) { const T q = a / b; const T c = a - q * b; a = b; b = c; - const T z = safe_mod(x - mul_mod(q, y, m), m); + // 初期値と呼び出し元が保証する互いに素という条件から 0 < q < m。 + const T qy = mul_mod_normalized(q, y, m); + const T z = x >= qy ? x - qy : x + (m - qy); x = y; y = z; } - assert(a == 1); - return x; + return y; } template std::pair merge_congruence(T r0, T m0, T r1, T m1) { - assert(0 <= r0 && r0 < m0); - assert(0 <= r1 && r1 < m1); + if (m0 == 1) { + return {r1, m1}; + } + if (m1 == 1) { + return {r0, m0}; + } const T g = gcd(m0, m1); const T diff = r1 - r0; - if (diff % g != 0) { + if (diff == 0) { + return {r0, m0 * (m1 / g)}; + } + const T diff_div_g = diff / g; + if (diff_div_g * g != diff) { return {0, 0}; } const T u1 = m1 / g; if (u1 == 1) { return {r0, m0}; } - const T t = mul_mod(diff / g, inv_mod(m0 / g, u1), u1); - r0 += t * m0; + const T factor = safe_mod(diff_div_g, u1); + if (factor != 0) { + const T t = + mul_mod_normalized(factor, inv_mod(safe_mod(m0 / g, u1), u1), u1); + r0 += t * m0; + } m0 *= u1; return {r0, m0}; } template std::pair solve_linear_congruence(T a, T b, T m) { - assert(m > 0); + if (m == 1) { + return {0, 1}; + } a = safe_mod(a, m); b = safe_mod(b, m); - const T g = gcd(a, m); - if (b % g != 0) { + if (a == 0) { + return b == 0 ? std::pair{0, 1} : std::pair{0, 0}; + } + if (a == 1) { + return {b, m}; + } + // 0 < a < m なので、引数順をこの向きにすると最初の剰余計算を省ける。 + const T g = gcd(m, a); + const T b_div_g = b / g; + if (b_div_g * g != b) { return {0, 0}; } const T mod = m / g; if (mod == 1) { return {0, 1}; } - const T rem = mul_mod(b / g, inv_mod(a / g, mod), mod); + if (b_div_g == 0) { + return {0, mod}; + } + const T rem = mul_mod_normalized(b_div_g, inv_mod(a / g, mod), mod); return {rem, mod}; } } // namespace generalized_garner_internal @@ -165,8 +207,8 @@ std::pair generalized_garner(const std::vector &a, assert(a.size() == m.size()); R r0 = 0; R m0 = 1; - const int n = static_cast(a.size()); - for (int i = 0; i < n; ++i) { + const std::size_t n = a.size(); + for (std::size_t i = 0; i < n; ++i) { const R mi = static_cast(m[i]); assert(mi > 0); const auto [r1, m1] = gg_internal::solve_linear_congruence( diff --git a/other/enumerate-maximum-independent-set-path-sums.hpp b/other/enumerate-maximum-independent-set-path-sums.hpp index 49a17f4..05b5da9 100644 --- a/other/enumerate-maximum-independent-set-path-sums.hpp +++ b/other/enumerate-maximum-independent-set-path-sums.hpp @@ -45,13 +45,13 @@ std::vector marginal_values(const std::vector &a) { while (stack.size() >= 2 && stack.back().first && removable_side(x, stack.back()) && removable_side(stack[stack.size() - 2], stack.back())) { - const P l = stack[stack.size() - 2]; - const P p = stack.back(); - res.push_back(p.second); - stack.pop_back(); - stack.pop_back(); + P l = std::move(stack[stack.size() - 2]); + P p = std::move(stack.back()); x = l.first && x.first ? P(true, l.second - p.second + x.second) : none; + res.push_back(std::move(p.second)); + stack.pop_back(); + stack.pop_back(); } stack.push_back(x); }; @@ -60,17 +60,24 @@ std::vector marginal_values(const std::vector &a) { add(P(true, x)); } add(none); - assert(res.size() == (a.size() + 1) / 2); return res; } -template std::size_t bucket_index(const T x, const T offset) { - assert(-offset <= x && x <= offset); - const std::size_t base = static_cast(offset); +// 呼び出し元で x の絶対値が offset 以下であることを保証する。 +template +std::size_t bucket_index(const T x, const std::size_t offset) { if (T(0) <= x) { - return base + static_cast(x); + return offset + static_cast(x); } - return base - static_cast(T(0) - x); + return offset - static_cast(T(0) - x); +} + +template +T bucket_value(const std::size_t index, const std::size_t offset) { + if (offset <= index) { + return static_cast(index - offset); + } + return T(0) - static_cast(offset - index); } } // namespace enumerate_maximum_independent_set_path_sums_internal @@ -134,37 +141,39 @@ std::vector enumerate_maximum_independent_set_path_sums_bucket_sort( std::vector count(offset * 2 + 1, 0); for (const T &x : xs) { count[enumerate_maximum_independent_set_path_sums_internal:: - bucket_index(x, total)] += 1; + bucket_index(x, offset)] += 1; } std::vector res(1, T()); res.reserve(xs.size() + 1); T sum = T(); if constexpr (maximum) { - for (T x = total;; x -= T(1)) { - const std::size_t c = - count[enumerate_maximum_independent_set_path_sums_internal:: - bucket_index(x, total)]; + for (std::size_t index_plus_one = count.size(); index_plus_one > 0; + index_plus_one -= 1) { + const std::size_t index = index_plus_one - 1; + const std::size_t c = count[index]; + if (c == 0) { + continue; + } + const T x = enumerate_maximum_independent_set_path_sums_internal:: + bucket_value(index, offset); for (std::size_t i = 0; i < c; i += 1) { sum += x; res.push_back(sum); } - if (x == -total) { - break; - } } } else { - for (T x = -total;; x += T(1)) { - const std::size_t c = - count[enumerate_maximum_independent_set_path_sums_internal:: - bucket_index(x, total)]; + for (std::size_t index = 0; index < count.size(); index += 1) { + const std::size_t c = count[index]; + if (c == 0) { + continue; + } + const T x = enumerate_maximum_independent_set_path_sums_internal:: + bucket_value(index, offset); for (std::size_t i = 0; i < c; i += 1) { sum += x; res.push_back(sum); } - if (x == total) { - break; - } } } return res; diff --git a/other/minimum-mod-range-increment-decrement-operations.hpp b/other/minimum-mod-range-increment-decrement-operations.hpp index f8b5f73..e27885e 100644 --- a/other/minimum-mod-range-increment-decrement-operations.hpp +++ b/other/minimum-mod-range-increment-decrement-operations.hpp @@ -5,7 +5,7 @@ // 区間に整数 x を加える操作の最小 sum |x| と同値である。 // a, b は同じ長さで、各要素は [0, m) を仮定する。 // T は 64 bit 整数型、m <= 10^9 を想定する。 -// 計算量 O(N log N)。 +// 平均計算量 O(N)。 #include #include @@ -14,8 +14,9 @@ #include template -T minimum_mod_range_increment_decrement_operations(std::vector a, - std::vector b, T m) { +T minimum_mod_range_increment_decrement_operations(const std::vector &a, + const std::vector &b, + T m) { static_assert(std::is_integral_v, "T must be integer."); static_assert(sizeof(T) >= sizeof(long long), "T must be at least 64-bit integer type."); @@ -65,9 +66,13 @@ T minimum_mod_range_increment_decrement_operations(std::vector a, assert(remainder == 0); assert(negative_count <= n + 1); - std::sort(differences.begin(), differences.end()); - const std::size_t keep_count = n + 1 - negative_count; + if (0 < keep_count && keep_count < differences.size()) { + std::nth_element(differences.begin(), + differences.begin() + + static_cast(keep_count), + differences.end()); + } T answer = 0; for (std::size_t i = 0; i < keep_count; ++i) { answer += differences[i]; diff --git a/other/rectangle-add-max-get.hpp b/other/rectangle-add-max-get.hpp index 6310ffa..2f9761d 100644 --- a/other/rectangle-add-max-get.hpp +++ b/other/rectangle-add-max-get.hpp @@ -10,6 +10,7 @@ #include #include +#include #include #include #include @@ -24,7 +25,7 @@ template struct Rectangle { C w; }; -template void assert_valid_weight(C w) { +template void assert_valid_weight([[maybe_unused]] const C &w) { if constexpr (std::numeric_limits::is_specialized && std::numeric_limits::is_integer && std::numeric_limits::is_signed) { @@ -32,8 +33,8 @@ template void assert_valid_weight(C w) { } } +// 呼び出し元で l <= r が保証されていることを仮定する。 template CoordinateLength coordinate_difference(T l, T r) { - assert(l <= r); return static_cast>(r) - static_cast>(l); } @@ -63,16 +64,16 @@ template struct SegmentTree { SegmentTree() = default; - explicit SegmentTree(const std::vector &leaves) { + // leaves が空でないことを仮定する。 + explicit SegmentTree(std::vector leaves) { const int n = static_cast(leaves.size()); - assert(n > 0); while (size < n) { size <<= 1; } data.assign(size << 1, Node{C(), Length(), T(), T()}); lazy.assign(size << 1, C()); for (int i = 0; i < n; ++i) { - data[size + i] = leaves[i]; + data[size + i] = std::move(leaves[i]); } for (int i = size - 1; i >= 1; --i) { data[i] = merge(data[i << 1], data[i << 1 | 1]); @@ -95,14 +96,20 @@ template struct SegmentTree { return Node{a.max_value, a.length + b.length, a.minimum_y, b.maximum_y}; } - void apply(int l, int r, C w, bool add) { apply(1, 0, size, l, r, w, add); } + void apply(int l, int r, const C &w, bool add) { + if (add) { + apply(1, 0, size, l, r, w); + } else { + apply(1, 0, size, l, r, w); + } + } const Node &all_prod() const { return data[1]; } Node prod(int l, int r) { return prod(1, 0, size, l, r); } - void add_node(int v, C w, bool add) { - if (add) { + template void add_node(int v, const C &w) { + if constexpr (Add) { data[v].max_value += w; lazy[v] += w; } else { @@ -115,23 +122,24 @@ template struct SegmentTree { if (lazy[v] == C()) { return; } - add_node(v << 1, lazy[v], true); - add_node(v << 1 | 1, lazy[v], true); + add_node(v << 1, lazy[v]); + add_node(v << 1 | 1, lazy[v]); lazy[v] = C(); } - void apply(int v, int l, int r, int ql, int qr, C w, bool add) { + template + void apply(int v, int l, int r, int ql, int qr, const C &w) { if (qr <= l || r <= ql) { return; } if (ql <= l && r <= qr) { - add_node(v, w, add); + add_node(v, w); return; } push(v); const int m = (l + r) >> 1; - apply(v << 1, l, m, ql, qr, w, add); - apply(v << 1 | 1, m, r, ql, qr, w, add); + apply(v << 1, l, m, ql, qr, w); + apply(v << 1 | 1, m, r, ql, qr, w); data[v] = merge(data[v << 1], data[v << 1 | 1]); } @@ -172,7 +180,7 @@ template struct CompressedRectangleAddMaxGet { return; } rectangle_add_max_get_internal::assert_valid_weight(w); - rectangles.emplace_back(Rectangle{l, d, r, u, w}); + rectangles.emplace_back(Rectangle{l, d, r, u, std::move(w)}); } std::tuple calc_max_lexicographically_minimum_point() const { @@ -256,7 +264,8 @@ template struct CompressedRectangleAddMaxGet { T d = rectangles[0].d; T r = rectangles[0].r; T u = rectangles[0].u; - for (const auto &rect : rectangles) { + for (std::size_t i = 1; i < rectangles.size(); ++i) { + const auto &rect = rectangles[i]; l = std::min(l, rect.l); d = std::min(d, rect.d); r = std::max(r, rect.r); @@ -275,6 +284,8 @@ template struct CompressedRectangleAddMaxGet { Result calc_impl(const QueryRange &query) const { std::vector xs{query.l, query.r}; std::vector ys{query.d, query.u}; + xs.reserve(rectangles.size() * 2 + 2); + ys.reserve(rectangles.size() * 2 + 2); std::vector clipped; clipped.reserve(rectangles.size()); for (const auto &rect : rectangles) { @@ -323,7 +334,7 @@ template struct CompressedRectangleAddMaxGet { ys[i], ys[i + 1]), ys[i], static_cast(ys[i + 1] - T(1))}); } - SegmentTree seg(leaves); + SegmentTree seg(std::move(leaves)); Result ret; bool found = false; @@ -356,16 +367,9 @@ template struct CompressedRectangleAddMaxGet { } found = true; } else if (ret.max_value == now.max_value) { - if (std::make_pair(minimum_x, now.minimum_y) < - std::make_pair(ret.minimum_x, ret.minimum_y)) { - ret.minimum_x = minimum_x; - ret.minimum_y = now.minimum_y; - } - if (std::make_pair(ret.maximum_x, ret.maximum_y) < - std::make_pair(maximum_x, now.maximum_y)) { - ret.maximum_x = maximum_x; - ret.maximum_y = now.maximum_y; - } + // x 座標を昇順に走査しているため、最小点は更新せず、最大点は常に更新する。 + ret.maximum_x = maximum_x; + ret.maximum_y = now.maximum_y; if constexpr (NeedArea) { ret.max_area += area; } @@ -397,7 +401,7 @@ template struct RectangleAddMaxGet { return; } rectangle_add_max_get_internal::assert_valid_weight(w); - rectangles.emplace_back(Rectangle{l, d, r, u, w}); + rectangles.emplace_back(Rectangle{l, d, r, u, std::move(w)}); } std::tuple calc_max_lexicographically_minimum_point() const { @@ -517,7 +521,8 @@ template struct RectangleAddMaxGet { T d = rectangles[0].d; T r = rectangles[0].r; T u = rectangles[0].u; - for (const auto &rect : rectangles) { + for (std::size_t i = 1; i < rectangles.size(); ++i) { + const auto &rect = rectangles[i]; l = std::min(l, rect.l); d = std::min(d, rect.d); r = std::max(r, rect.r); @@ -566,7 +571,7 @@ template struct RectangleAddMaxGet { } static void add_event(std::vector &head, std::vector &events, - int x, int d, int u, C w) { + int x, int d, int u, const C &w) { events.emplace_back(Event{head[x], d, u, w}); head[x] = static_cast(events.size()) - 1; } @@ -584,13 +589,13 @@ template struct RectangleAddMaxGet { if (cl >= cr || cd >= cu) { continue; } - const int xl = rectangle_add_max_get_internal::checked_size( + const int xl = static_cast( rectangle_add_max_get_internal::coordinate_difference(l, cl)); - const int xr = rectangle_add_max_get_internal::checked_size( + const int xr = static_cast( rectangle_add_max_get_internal::coordinate_difference(l, cr)); - const int yd = rectangle_add_max_get_internal::checked_size( + const int yd = static_cast( rectangle_add_max_get_internal::coordinate_difference(d, cd)); - const int yu = rectangle_add_max_get_internal::checked_size( + const int yu = static_cast( rectangle_add_max_get_internal::coordinate_difference(d, cu)); add_event(event_list.add_head, event_list.events, xl, yd, yu, rect.w); @@ -627,7 +632,7 @@ template struct RectangleAddMaxGet { } template - static void update_result(Result &ret, bool &found, C value, + static void update_result(Result &ret, bool &found, const C &value, T minimum_x, T minimum_y, T maximum_x, T maximum_y, T2 area) { if (!found || ret.max_value < value) { @@ -641,16 +646,9 @@ template struct RectangleAddMaxGet { } found = true; } else if (ret.max_value == value) { - if (std::make_pair(minimum_x, minimum_y) < - std::make_pair(ret.minimum_x, ret.minimum_y)) { - ret.minimum_x = minimum_x; - ret.minimum_y = minimum_y; - } - if (std::make_pair(ret.maximum_x, ret.maximum_y) < - std::make_pair(maximum_x, maximum_y)) { - ret.maximum_x = maximum_x; - ret.maximum_y = maximum_y; - } + // x 座標を昇順に走査しているため、最小点は更新せず、最大点は常に更新する。 + ret.maximum_x = maximum_x; + ret.maximum_y = maximum_y; if constexpr (NeedArea) { ret.max_area += area; } @@ -682,7 +680,6 @@ template struct RectangleAddMaxGet { update_result(ret, found, now.max_value, x, now.minimum_y, x, now.maximum_y, area); } - assert(found); return ret; } @@ -703,10 +700,10 @@ template struct RectangleAddMaxGet { if (query.lower_y[i] == query.upper_y[i]) { continue; } - const int d = rectangle_add_max_get_internal::checked_size( + const int d = static_cast( rectangle_add_max_get_internal::coordinate_difference( query.d, query.lower_y[i])); - const int u = rectangle_add_max_get_internal::checked_size( + const int u = static_cast( rectangle_add_max_get_internal::coordinate_difference( query.d, query.upper_y[i])); const Node now = seg.prod(d, u); @@ -718,7 +715,6 @@ template struct RectangleAddMaxGet { update_result(ret, found, now.max_value, x, now.minimum_y, x, now.maximum_y, area); } - assert(found); return ret; } }; diff --git a/structure/others/dynamic-median.hpp b/structure/others/dynamic-median.hpp index a257187..a9736f6 100644 --- a/structure/others/dynamic-median.hpp +++ b/structure/others/dynamic-median.hpp @@ -10,6 +10,7 @@ #include #include #include +#include enum class DynamicMedianMode { Lower, Upper, Average }; @@ -19,29 +20,26 @@ template struct DynamicMedian { void add(T x) { if (lower_values.empty() || !(*lower_values.rbegin() < x)) { - lower_values.insert(x); + lower_values.insert(std::move(x)); } else { - upper_values.insert(x); + upper_values.insert(std::move(x)); } balance(); } - bool erase(T x) { - auto lower_itr = lower_values.find(x); - if (lower_itr != lower_values.end()) { - lower_values.erase(lower_itr); - balance(); - return true; + bool erase(const T &x) { + if (lower_values.empty()) { + return false; } - auto upper_itr = upper_values.find(x); - if (upper_itr != upper_values.end()) { - upper_values.erase(upper_itr); - balance(); - return true; + auto &values = *lower_values.rbegin() < x ? upper_values : lower_values; + const auto itr = values.find(x); + if (itr == values.end()) { + return false; } - - return false; + values.erase(itr); + balance(); + return true; } template @@ -83,15 +81,13 @@ template struct DynamicMedian { private: void balance() { - while (lower_values.size() < upper_values.size()) { - lower_values.insert(*upper_values.begin()); - upper_values.erase(upper_values.begin()); - } - while (lower_values.size() > upper_values.size() + 1) { + // add または erase の直前に不変条件が成り立つため、移動は高々 1 回でよい。 + if (lower_values.size() < upper_values.size()) { + lower_values.insert(upper_values.extract(upper_values.begin())); + } else if (lower_values.size() > upper_values.size() + 1) { auto itr = lower_values.end(); --itr; - upper_values.insert(*itr); - lower_values.erase(itr); + upper_values.insert(lower_values.extract(itr)); } } }; diff --git a/verify/aoj-cgl-4-c.test.cpp b/verify/aoj-cgl-4-c.test.cpp index 52d4b36..10e81ae 100644 --- a/verify/aoj-cgl-4-c.test.cpp +++ b/verify/aoj-cgl-4-c.test.cpp @@ -5,6 +5,7 @@ #include #include #include +#include #include #include "../geometry/line-convex-polygon-intersection.hpp" @@ -16,9 +17,6 @@ struct Point { Point() : x(0), y(0) {} Point(long long x_, long long y_) : x(x_), y(y_) {} - Point operator+(const Point &other) const { - return Point(x + other.x, y + other.y); - } Point operator-(const Point &other) const { return Point(x - other.x, y - other.y); } @@ -28,16 +26,7 @@ struct RealPoint { long double x; long double y; - RealPoint() : x(0), y(0) {} RealPoint(long double x_, long double y_) : x(x_), y(y_) {} - - RealPoint operator+(const RealPoint &other) const { - return RealPoint(x + other.x, y + other.y); - } - RealPoint operator-(const RealPoint &other) const { - return RealPoint(x - other.x, y - other.y); - } - RealPoint operator*(long double k) const { return RealPoint(x * k, y * k); } }; __int128_t cross_value(const Point &a, const Point &b) { @@ -74,78 +63,28 @@ long double dot_ld(const RealPoint &a, const RealPoint &b) { return a.x * b.x + a.y * b.y; } -bool same_point(const RealPoint &a, const RealPoint &b) { - return std::abs(a.x - b.x) <= 1e-12L && std::abs(a.y - b.y) <= 1e-12L; -} - -void add_unique(std::vector &points, const RealPoint &point) { - for (const auto &q : points) { - if (same_point(q, point)) { - return; - } - } - points.push_back(point); -} - -RealPoint line_intersection(const Point &line_a, const Point &line_b, - const Point &segment_a, const Point &segment_b) { - const RealPoint u = to_real_point(segment_a); - const RealPoint v = to_real_point(segment_b); - const RealPoint dir = to_real_point(line_b - line_a); - const RealPoint edge = v - u; - const long double denominator = cross_ld(edge, dir); - const long double numerator = - cross_ld(to_real_point(line_a - segment_a), dir); - const long double t = numerator / denominator; - return u + edge * t; -} - -std::vector brute_intersections(const std::vector &polygon, - const Point &line_a, - const Point &line_b) { +std::vector convex_cut(const std::vector &polygon, + const Point &line_a, const Point &line_b) { const int n = static_cast(polygon.size()); const Point direction = line_b - line_a; - std::vector points; - for (int i = 0; i < n; ++i) { - const Point &u = polygon[i]; - const Point &v = polygon[(i + 1) % n]; - const __int128_t hu = cross_value(u - line_a, direction); - const __int128_t hv = cross_value(v - line_a, direction); - const int su = sign_value(hu); - const int sv = sign_value(hv); - if (su == 0 && sv == 0) { - add_unique(points, to_real_point(u)); - add_unique(points, to_real_point(v)); - } else if (su == 0) { - add_unique(points, to_real_point(u)); - } else if (sv == 0) { - add_unique(points, to_real_point(v)); - } else if (su != sv) { - add_unique(points, line_intersection(line_a, line_b, u, v)); - } - } - if (points.size() <= 2) { - return points; + const RealPoint real_direction = to_real_point(direction); + + const auto exact_intersections = + line_polygon_intersection(polygon, line_a, line_b); + std::vector intersections; + intersections.reserve(exact_intersections.size()); + for (const auto &point : exact_intersections) { + intersections.emplace_back(point.template x_as(), + point.template y_as()); } - const RealPoint dir = to_real_point(direction); - int min_index = 0; - int max_index = 0; - for (int i = 1; i < static_cast(points.size()); ++i) { - if (dot_ld(points[i], dir) < dot_ld(points[min_index], dir)) { - min_index = i; - } - if (dot_ld(points[i], dir) > dot_ld(points[max_index], dir)) { - max_index = i; - } + if (intersections.size() == 2 && + dot_ld(intersections[1], real_direction) < + dot_ld(intersections[0], real_direction)) { + std::swap(intersections[0], intersections[1]); } - return {points[min_index], points[max_index]}; -} -std::vector convex_cut(const std::vector &polygon, - const Point &line_a, const Point &line_b) { - const int n = static_cast(polygon.size()); - const Point direction = line_b - line_a; std::vector result; + result.reserve(n + 2); for (int i = 0; i < n; ++i) { const Point &now = polygon[i]; const Point &next = polygon[(i + 1) % n]; @@ -157,9 +96,12 @@ std::vector convex_cut(const std::vector &polygon, if (sign_value(current) >= 0) { result.push_back(to_real_point(now)); } - if ((current < 0 && next_value > 0) || - (current > 0 && next_value < 0)) { - result.push_back(line_intersection(line_a, line_b, now, next)); + if (current < 0 && next_value > 0) { + assert(intersections.size() == 2); + result.push_back(intersections[1]); + } else if (current > 0 && next_value < 0) { + assert(intersections.size() == 2); + result.push_back(intersections[0]); } } return result; @@ -192,31 +134,6 @@ remove_collinear_vertices(const std::vector &polygon) { return result; } -void verify_intersection(const std::vector &polygon, const Point &line_a, - const Point &line_b) { - const auto fast = line_polygon_intersection(polygon, line_a, line_b); - const auto slow = brute_intersections(polygon, line_a, line_b); - assert(fast.size() == slow.size()); - if (fast.empty()) { - return; - } - - auto to_real = [](const auto &point) { - return RealPoint(point.template x_as(), - point.template y_as()); - }; - - if (fast.size() == 1) { - assert(same_point(to_real(fast[0]), slow[0])); - return; - } - const RealPoint p0 = to_real(fast[0]); - const RealPoint p1 = to_real(fast[1]); - const bool ok = (same_point(p0, slow[0]) && same_point(p1, slow[1])) || - (same_point(p0, slow[1]) && same_point(p1, slow[0])); - assert(ok); -} - int main() { std::ios::sync_with_stdio(false); std::cin.tie(nullptr); @@ -235,7 +152,6 @@ int main() { while (q--) { const Point line_a = read_point(); const Point line_b = read_point(); - verify_intersection(polygon, line_a, line_b); std::cout << polygon_area(convex_cut(polygon, line_a, line_b)) << '\n'; } diff --git a/verify/standalone-01-on-tree.test.cpp b/verify/standalone-01-on-tree.test.cpp new file mode 100644 index 0000000..e1a5ab7 --- /dev/null +++ b/verify/standalone-01-on-tree.test.cpp @@ -0,0 +1,125 @@ +// competitive-verifier: STANDALONE + +#include +#include +#include +#include +#include + +#include "../graph/tree/01-on-tree.hpp" + +long long brute_01_on_tree(int n, const std::vector> &edges, + const std::vector &c0, + const std::vector &c1, int root) { + std::vector> graph(n); + for (const auto &[u, v] : edges) { + graph[u].push_back(v); + graph[v].push_back(u); + } + + std::vector parent(n, -2), stack{root}; + parent[root] = -1; + while (!stack.empty()) { + const int v = stack.back(); + stack.pop_back(); + for (int to : graph[v]) { + if (to == parent[v]) { + continue; + } + parent[to] = v; + stack.push_back(to); + } + } + + std::vector order(n), position(n); + for (int v = 0; v < n; ++v) { + order[v] = v; + } + + long long answer = std::numeric_limits::max(); + do { + for (int i = 0; i < n; ++i) { + position[order[i]] = i; + } + bool valid = true; + for (int v = 0; v < n; ++v) { + if (parent[v] != -1 && position[parent[v]] > position[v]) { + valid = false; + } + } + if (!valid) { + continue; + } + + long long current = 0; + for (int i = 0; i < n; ++i) { + for (int j = i + 1; j < n; ++j) { + current += c1[order[i]] * c0[order[j]]; + } + } + answer = std::min(answer, current); + } while (std::next_permutation(order.begin(), order.end())); + return answer; +} + +void test_01_on_tree(int n, const std::vector> &edges, + const std::vector &c0, + const std::vector &c1, int root) { + const long long expected = brute_01_on_tree(n, edges, c0, c1, root); + const long long actual = solve_01_on_tree(n, edges, c0, c1, root); + assert(actual == expected); +} + +void test_all_small_counts(int n, + const std::vector> &edges) { + int pattern_count = 1; + for (int i = 0; i < 2 * n; ++i) { + pattern_count *= 3; + } + for (int code = 0; code < pattern_count; ++code) { + int current_code = code; + std::vector c0(n), c1(n); + for (int v = 0; v < n; ++v) { + c0[v] = current_code % 3; + current_code /= 3; + c1[v] = current_code % 3; + current_code /= 3; + } + test_01_on_tree(n, edges, c0, c1, 0); + } +} + +void enumerate_parent_trees(int n, int v, + std::vector> &edges) { + if (v == n) { + test_all_small_counts(n, edges); + return; + } + for (int parent = 0; parent < v; ++parent) { + edges.push_back({parent, v}); + enumerate_parent_trees(n, v + 1, edges); + edges.pop_back(); + } +} + +int main() { + for (int n = 1; n <= 4; ++n) { + std::vector> edges; + enumerate_parent_trees(n, 1, edges); + } + + const std::vector> edges{ + {0, 1}, {0, 2}, {2, 3}, {2, 4}}; + const std::vector c0{0, 3, 0, 2, 1}; + const std::vector c1{0, 0, 4, 1, 0}; + for (int root = 0; root < 5; ++root) { + test_01_on_tree(5, edges, c0, c1, root); + } + + const std::vector large_c0{0, 4'000'000'000ULL}; + const std::vector large_c1{4'000'000'000ULL, 0}; + const auto large_answer = solve_01_on_tree(2, {{0, 1}}, large_c0, large_c1); + assert(large_answer == 16'000'000'000'000'000'000ULL); + + return 0; +} diff --git a/verify/standalone-determinant-of-linear-matrix-polynomial.test.cpp b/verify/standalone-determinant-of-linear-matrix-polynomial.test.cpp new file mode 100644 index 0000000..af8a029 --- /dev/null +++ b/verify/standalone-determinant-of-linear-matrix-polynomial.test.cpp @@ -0,0 +1,216 @@ +// competitive-verifier: STANDALONE + +#include +#include + +#include "../math/matrix/determinant-of-linear-matrix-polynomial.hpp" + +namespace { +struct ModInt101 { + static constexpr int mod = 101; + int value; + + ModInt101(long long value = 0) { + value %= mod; + if (value < 0) { + value += mod; + } + this->value = static_cast(value); + } + + ModInt101 &operator+=(const ModInt101 &other) { + value += other.value; + if (value >= mod) { + value -= mod; + } + return *this; + } + + ModInt101 &operator-=(const ModInt101 &other) { + value -= other.value; + if (value < 0) { + value += mod; + } + return *this; + } + + ModInt101 &operator*=(const ModInt101 &other) { + value = value * other.value % mod; + return *this; + } + + ModInt101 &operator/=(const ModInt101 &other) { + assert(other != ModInt101()); + return *this *= power(other, mod - 2); + } + + friend ModInt101 operator-(ModInt101 lhs, const ModInt101 &rhs) { + return lhs -= rhs; + } + + friend ModInt101 operator*(ModInt101 lhs, const ModInt101 &rhs) { + return lhs *= rhs; + } + + friend ModInt101 operator/(ModInt101 lhs, const ModInt101 &rhs) { + return lhs /= rhs; + } + + friend bool operator==(const ModInt101 &lhs, const ModInt101 &rhs) { + return lhs.value == rhs.value; + } + + friend bool operator!=(const ModInt101 &lhs, const ModInt101 &rhs) { + return lhs.value != rhs.value; + } + + private: + static ModInt101 power(ModInt101 base, int exponent) { + ModInt101 result(1); + while (exponent > 0) { + if (exponent % 2 == 1) { + result *= base; + } + base *= base; + exponent /= 2; + } + return result; + } +}; + +using Matrix = std::vector>; + +ModInt101 determinant(Matrix matrix) { + const int n = static_cast(matrix.size()); + ModInt101 result(1); + for (int column = 0; column < n; ++column) { + int pivot = -1; + for (int row = column; row < n; ++row) { + if (matrix[row][column] != ModInt101()) { + pivot = row; + break; + } + } + if (pivot < 0) { + return ModInt101(); + } + if (pivot != column) { + matrix[pivot].swap(matrix[column]); + result = ModInt101() - result; + } + + result *= matrix[column][column]; + const ModInt101 pivot_inverse = ModInt101(1) / matrix[column][column]; + for (int row = column + 1; row < n; ++row) { + const ModInt101 factor = matrix[row][column] * pivot_inverse; + for (int j = column + 1; j < n; ++j) { + matrix[row][j] -= matrix[column][j] * factor; + } + } + } + return result; +} + +ModInt101 evaluate(const std::vector &polynomial, ModInt101 x) { + ModInt101 result; + for (int i = static_cast(polynomial.size()) - 1; i >= 0; --i) { + result *= x; + result += polynomial[i]; + } + return result; +} + +Matrix matrix_from_mask(int n, int mask) { + Matrix matrix(n, std::vector(n)); + for (int i = 0; i < n; ++i) { + for (int j = 0; j < n; ++j) { + matrix[i][j] = ModInt101((mask >> (i * n + j)) % 2); + } + } + return matrix; +} + +void check_characteristic_polynomial(const Matrix &matrix) { + const int n = static_cast(matrix.size()); + const std::vector polynomial = characteristic_polynomial(matrix); + assert(static_cast(polynomial.size()) == n + 1); + + Matrix hessenberg = matrix; + hessenberg_reduction(hessenberg); + for (int i = 0; i < n; ++i) { + for (int j = 0; j + 1 < i; ++j) { + assert(hessenberg[i][j] == ModInt101()); + } + } + + for (int x = 0; x <= n; ++x) { + Matrix shifted = matrix; + Matrix shifted_hessenberg = hessenberg; + for (int i = 0; i < n; ++i) { + for (int j = 0; j < n; ++j) { + const ModInt101 diagonal = i == j ? ModInt101(x) : ModInt101(); + shifted[i][j] = diagonal - shifted[i][j]; + shifted_hessenberg[i][j] = diagonal - shifted_hessenberg[i][j]; + } + } + const ModInt101 expected = determinant(shifted); + assert(evaluate(polynomial, ModInt101(x)) == expected); + assert(determinant(shifted_hessenberg) == expected); + } +} + +void check_linear_matrix_polynomial(const Matrix &matrix_0, + const Matrix &matrix_1) { + const int n = static_cast(matrix_0.size()); + const std::vector polynomial = + determinant_of_linear_matrix_polynomial(matrix_0, matrix_1); + assert(static_cast(polynomial.size()) == n + 1); + + for (int x = 0; x <= n; ++x) { + Matrix evaluated = matrix_0; + for (int i = 0; i < n; ++i) { + for (int j = 0; j < n; ++j) { + evaluated[i][j] += ModInt101(x) * matrix_1[i][j]; + } + } + assert(evaluate(polynomial, ModInt101(x)) == determinant(evaluated)); + } +} + +void self_test() { + check_characteristic_polynomial({}); + check_linear_matrix_polynomial({}, {}); + + for (int n = 1; n <= 3; ++n) { + const int state_count = 1 << (n * n); + for (int mask = 0; mask < state_count; ++mask) { + check_characteristic_polynomial(matrix_from_mask(n, mask)); + } + } + + for (int n = 1; n <= 2; ++n) { + const int state_count = 1 << (n * n); + for (int mask_0 = 0; mask_0 < state_count; ++mask_0) { + for (int mask_1 = 0; mask_1 < state_count; ++mask_1) { + check_linear_matrix_polynomial(matrix_from_mask(n, mask_0), + matrix_from_mask(n, mask_1)); + } + } + } + + constexpr int state_count = 1 << 9; + for (int mask = 0; mask < state_count; ++mask) { + check_linear_matrix_polynomial(matrix_from_mask(3, mask), + matrix_from_mask(3, 0)); + check_linear_matrix_polynomial( + matrix_from_mask(3, mask), + matrix_from_mask(3, (mask * 137 + 91) & (state_count - 1))); + } +} +} // namespace + +int main() { + self_test(); + + return 0; +} diff --git a/verify/standalone-dynamic-matrix-rank.test.cpp b/verify/standalone-dynamic-matrix-rank.test.cpp index 9a4db09..ffdf1cd 100644 --- a/verify/standalone-dynamic-matrix-rank.test.cpp +++ b/verify/standalone-dynamic-matrix-rank.test.cpp @@ -215,7 +215,7 @@ void self_test() { } for (int h = 1; h <= 3; ++h) { - for (int w = 1; w <= 3; ++w) { + for (int w = 0; w <= 3; ++w) { const int states = 1 << (h * w); for (int mask = 0; mask < states; ++mask) { check_updates(matrix_from_mask(h, w, mask)); diff --git a/verify/standalone-dynamic-median.test.cpp b/verify/standalone-dynamic-median.test.cpp index 4395af7..7da23a0 100644 --- a/verify/standalone-dynamic-median.test.cpp +++ b/verify/standalone-dynamic-median.test.cpp @@ -8,6 +8,32 @@ #include "../structure/others/dynamic-median.hpp" namespace { +struct CopyCountedValue { + static inline int copy_count = 0; + + int value; + + explicit CopyCountedValue(int value_) : value(value_) {} + + CopyCountedValue(const CopyCountedValue &other) : value(other.value) { + copy_count += 1; + } + + CopyCountedValue(CopyCountedValue &&) noexcept = default; + + CopyCountedValue &operator=(const CopyCountedValue &other) { + value = other.value; + copy_count += 1; + return *this; + } + + CopyCountedValue &operator=(CopyCountedValue &&) noexcept = default; + + bool operator<(const CopyCountedValue &other) const { + return value < other.value; + } +}; + template std::vector sorted_values(std::vector values) { std::sort(values.begin(), values.end()); return values; @@ -92,9 +118,32 @@ void check_known_cases() { assert(signed_median.median(DynamicMedianMode::Average) == 0); } +void check_value_transfers() { + DynamicMedian median; + const CopyCountedValue values[] = {CopyCountedValue(4), CopyCountedValue(5), + CopyCountedValue(3), + CopyCountedValue(2)}; + for (const auto &value : values) { + CopyCountedValue::copy_count = 0; + median.add(value); + assert(CopyCountedValue::copy_count == 1); + } + + const CopyCountedValue five(5); + CopyCountedValue::copy_count = 0; + assert(median.erase(five)); + assert(CopyCountedValue::copy_count == 0); + + const CopyCountedValue four(4); + CopyCountedValue::copy_count = 0; + assert(median.erase(four)); + assert(CopyCountedValue::copy_count == 0); +} + void self_test() { check_known_cases(); check_scripted_operations(); + check_value_transfers(); } } // namespace diff --git a/verify/standalone-floor-sum.test.cpp b/verify/standalone-floor-sum.test.cpp new file mode 100644 index 0000000..3e8bf75 --- /dev/null +++ b/verify/standalone-floor-sum.test.cpp @@ -0,0 +1,57 @@ +// competitive-verifier: STANDALONE + +#include +#include + +#include "../internal/int128.hpp" +#include "../math/number-theory/floor-sum.hpp" + +long long floor_div_brute(long long x, long long m) { + assert(m > 0); + + long long q = x / m; + const long long r = x % m; + if (r < 0) { + --q; + } + return q; +} + +long long floor_sum_brute(long long n, long long m, long long a, long long b) { + long long result = 0; + for (long long i = 0; i < n; ++i) { + result += floor_div_brute(a * i + b, m); + } + return result; +} + +int main() { + for (long long n = 0; n <= 20; ++n) { + for (long long m = 1; m <= 20; ++m) { + for (long long a = -20; a <= 20; ++a) { + for (long long b = -20; b <= 20; ++b) { + const long long expected = floor_sum_brute(n, m, a, b); + assert(floor_sum(n, m, a, b) == expected); + + using NicheLibrary::Int128; + const Int128 actual = + floor_sum(Int128(n), Int128(m), Int128(a), Int128(b)); + assert(actual == Int128(expected)); + } + } + } + } + for (std::uint64_t n = 0; n <= 20; ++n) { + for (std::uint64_t m = 1; m <= 20; ++m) { + for (std::uint64_t a = 0; a <= 20; ++a) { + for (std::uint64_t b = 0; b <= 20; ++b) { + const long long expected = floor_sum_brute(n, m, a, b); + assert(floor_sum(n, m, a, b) == + static_cast(expected)); + } + } + } + } + + return 0; +} diff --git a/verify/standalone-generalized-garner.test.cpp b/verify/standalone-generalized-garner.test.cpp index 85f044b..04b20cf 100644 --- a/verify/standalone-generalized-garner.test.cpp +++ b/verify/standalone-generalized-garner.test.cpp @@ -1,8 +1,10 @@ // competitive-verifier: STANDALONE #include +#include #include +#include "../internal/int128.hpp" #include "../math/number-theory/generalized-garner.hpp" namespace { @@ -114,12 +116,12 @@ void self_test_small() { } } -#ifdef __SIZEOF_INT128__ void self_test_int128() { - using i128 = __int128_t; + using i128 = NicheLibrary::Int128; const i128 x0 = 123456789; - const i128 m0 = i128(1) << 70; - const i128 m1 = i128(1) << 69; + const i128 m0 = i128::from_words(std::uint64_t{1} << 6, 0); + const i128 m1 = i128::from_words(std::uint64_t{1} << 5, 0); + const i128 m2 = 1000003; { const std::vector a = {3, 5}; const std::vector b = {3 * x0, 5 * x0}; @@ -128,6 +130,23 @@ void self_test_int128() { assert(res.first == x0); assert(res.second == m0); } + { + const std::vector a = {3, 5}; + const std::vector b = {3 * x0, 5 * x0}; + const std::vector m = {m0, m2}; + const auto res = generalized_garner(a, b, m); + assert(res.first == x0); + assert(res.second == m0 * m2); + } + { + const std::vector a = {1, 1}; + const std::vector b = {5, 3}; + const std::vector m = {m0, m2}; + const auto res = generalized_garner(a, b, m); + assert(res.first % m0 == 5); + assert(res.first % m2 == 3); + assert(res.second == m0 * m2); + } { const std::vector a = {1, 1}; const std::vector b = {0, 1}; @@ -137,14 +156,11 @@ void self_test_int128() { assert(res.second == 0); } } -#endif } // namespace int main() { self_test_small(); -#ifdef __SIZEOF_INT128__ self_test_int128(); -#endif return 0; } diff --git a/verify/standalone-graph-isomorphism.test.cpp b/verify/standalone-graph-isomorphism.test.cpp index d4dab0a..73f4069 100644 --- a/verify/standalone-graph-isomorphism.test.cpp +++ b/verify/standalone-graph-isomorphism.test.cpp @@ -147,6 +147,20 @@ void self_test() { } } + const std::vector> cycle_6{{0, 1}, {1, 2}, {2, 3}, + {3, 4}, {4, 5}, {5, 0}}; + const std::vector> two_triangles{ + {0, 1}, {1, 2}, {2, 0}, {3, 4}, {4, 5}, {5, 3}}; + GraphIsomorphism regular_graphs(6, cycle_6, two_triangles); + assert(!regular_graphs.run()); + assert(!regular_graphs.run()); + + const std::vector> complete_bipartite_3_3{ + {0, 3}, {0, 4}, {0, 5}, {1, 3}, {1, 4}, {1, 5}, {2, 3}, {2, 4}, {2, 5}}; + const std::vector> triangular_prism{ + {0, 1}, {1, 2}, {2, 0}, {3, 4}, {4, 5}, {5, 3}, {0, 3}, {1, 4}, {2, 5}}; + assert(!is_graph_isomorphic(6, complete_bipartite_3_3, triangular_prism)); + const std::vector> edges{{0, 1}, {1, 2}, {2, 3}, {3, 4}, {4, 0}, {2, 2}}; const std::vector permutation{2, 4, 1, 3, 0}; diff --git a/verify/standalone-int128.test.cpp b/verify/standalone-int128.test.cpp index 5564a0e..7f9781b 100644 --- a/verify/standalone-int128.test.cpp +++ b/verify/standalone-int128.test.cpp @@ -134,6 +134,13 @@ constexpr bool test_uint128_div_mod_fast_paths() { return false; } + UInt128::div_mod(UInt128::from_words(1, 0), UInt128(2), quotient, + remainder); + if (quotient != UInt128::from_words(0, std::uint64_t{1} << 63) || + remainder != UInt128{}) { + return false; + } + const UInt128 smaller = UInt128::from_words(1, ~std::uint64_t{}); const UInt128 larger = UInt128::from_words(2, 0); UInt128::div_mod(smaller, larger, quotient, remainder); @@ -161,6 +168,16 @@ int main() { "Int128 must be signed."); static_assert(std::numeric_limits::digits == 127, "Int128 must have 127 value bits."); + static_assert(static_cast(UInt128(123)) == 123.0L, + "UInt128 must convert to long double."); + static_assert(static_cast(Int128(-123)) == -123.0L, + "Int128 must convert to long double."); + static_assert(static_cast(UInt128::from_words(1, 0)) == + 0x1p64L, + "wide UInt128 must convert to long double."); + static_assert(static_cast( + Int128::from_words(~std::uint64_t{}, 0)) == -0x1p64L, + "wide negative Int128 must convert to long double."); const UInt128 unsigned_values[] = { UInt128(0), diff --git a/verify/standalone-line-convex-polygon-intersection.test.cpp b/verify/standalone-line-convex-polygon-intersection.test.cpp index ee03c3b..1bbcbee 100644 --- a/verify/standalone-line-convex-polygon-intersection.test.cpp +++ b/verify/standalone-line-convex-polygon-intersection.test.cpp @@ -1,6 +1,7 @@ // competitive-verifier: STANDALONE #include +#include #include #include @@ -66,5 +67,33 @@ int main() { assert(result.size() == 2); } + { + std::vector hull = {Point(0, 0), Point(4, 0), Point(0, 4)}; + const auto result = + line_convex_hull_intersection(hull, Point(0, 1), Point(3, 2)); + + assert(result.size() == 2); + auto is_point = [&](int index, Int128 x_numerator, Int128 y_numerator, + Int128 denominator) { + return result[index].x_numerator == x_numerator && + result[index].y_numerator == y_numerator && + result[index].denominator == denominator; + }; + const bool expected_order = + is_point(0, Int128(0), Int128(1), Int128(1)) && + is_point(1, Int128(9), Int128(7), Int128(4)); + const bool reverse_order = + is_point(1, Int128(0), Int128(1), Int128(1)) && + is_point(0, Int128(9), Int128(7), Int128(4)); + assert(expected_order || reverse_order); + + const int rational_index = result[0].denominator == Int128(4) ? 0 : 1; + const auto real_point = + result[rational_index] + .template to_point>(); + assert(std::abs(real_point.real() - 2.25L) < 1e-15L); + assert(std::abs(real_point.imag() - 1.75L) < 1e-15L); + } + return 0; } diff --git a/verify/standalone-minimum-mod-range-increment-decrement-operations.test.cpp b/verify/standalone-minimum-mod-range-increment-decrement-operations.test.cpp index 176e198..1a8e946 100644 --- a/verify/standalone-minimum-mod-range-increment-decrement-operations.test.cpp +++ b/verify/standalone-minimum-mod-range-increment-decrement-operations.test.cpp @@ -74,6 +74,10 @@ long long mod_difference_ll(long long from, long long to, long long mod) { } void self_test() { + const std::vector empty; + assert(minimum_mod_range_increment_decrement_operations(empty, empty, + 1LL) == 0); + for (int n = 1; n <= 4; ++n) { for (int mod = 1; mod <= 3; ++mod) { int total = 1; diff --git a/verify/standalone-prime-counting-modulo.test.cpp b/verify/standalone-prime-counting-modulo.test.cpp index 639ae4c..1421356 100644 --- a/verify/standalone-prime-counting-modulo.test.cpp +++ b/verify/standalone-prime-counting-modulo.test.cpp @@ -74,6 +74,17 @@ void self_test() { } } } + + const auto table = prime_counting_modulo_table(500, 7).second; + const auto converted = + prime_counting_modulo_mf_prefix_sum_table(500, 7); + assert(converted.size() == table.size()); + for (int r = 0; r < static_cast(table.size()); ++r) { + assert(converted[r].size() == table[r].size()); + for (int i = 0; i < static_cast(table[r].size()); ++i) { + assert(converted[r][i] == table[r][i]); + } + } } } // namespace diff --git a/verify/yosupo-rooted-tree-topological-order-with-minimum-inversions.test.cpp b/verify/yosupo-rooted-tree-topological-order-with-minimum-inversions.test.cpp index 39b519a..c8b9630 100644 --- a/verify/yosupo-rooted-tree-topological-order-with-minimum-inversions.test.cpp +++ b/verify/yosupo-rooted-tree-topological-order-with-minimum-inversions.test.cpp @@ -19,7 +19,9 @@ build_01_on_tree_order(int n, const std::vector> &edges, std::vector parent(n, -2); parent[root] = -1; - std::vector stack{root}; + std::vector stack; + stack.reserve(n); + stack.push_back(root); while (!stack.empty()) { const int v = stack.back(); stack.pop_back(); @@ -27,9 +29,6 @@ build_01_on_tree_order(int n, const std::vector> &edges, if (to == parent[v]) { continue; } - if (parent[to] != -2) { - continue; - } parent[to] = v; stack.push_back(to); } @@ -44,6 +43,9 @@ build_01_on_tree_order(int n, const std::vector> &edges, struct QueueCompare { bool operator()(const QueueNode &a, const QueueNode &b) const { + if (a.zero == b.zero && a.one == b.one) { + return a.vertex < b.vertex; + } if (a.one == 0 && b.one == 0) { return a.vertex < b.vertex; } @@ -53,6 +55,15 @@ build_01_on_tree_order(int n, const std::vector> &edges, if (b.one == 0) { return true; } + if (a.zero == 0 && b.zero == 0) { + return a.vertex < b.vertex; + } + if (a.zero == 0) { + return true; + } + if (b.zero == 0) { + return false; + } const long long lhs = a.zero * b.one; const long long rhs = b.zero * a.one; if (lhs < rhs) { @@ -65,8 +76,8 @@ build_01_on_tree_order(int n, const std::vector> &edges, } }; - std::vector leader(n), up = parent, version(n, 0), head(n), tail(n), - next(n, -1); + std::vector leader(n), up = std::move(parent), version(n, 0), head(n), + tail(n), next(n, -1); std::vector zero = c0, one = c1; for (int v = 0; v < n; ++v) { leader[v] = v; @@ -82,7 +93,10 @@ build_01_on_tree_order(int n, const std::vector> &edges, return v; }; - std::priority_queue, QueueCompare> que; + std::vector queue_storage; + queue_storage.reserve(n); + std::priority_queue, QueueCompare> que( + QueueCompare{}, std::move(queue_storage)); auto push = [&](int v) { if (up[v] != -1) { que.push(QueueNode{v, version[v], zero[v], one[v]}); @@ -94,17 +108,18 @@ build_01_on_tree_order(int n, const std::vector> &edges, } for (int merge_count = 0; merge_count < n - 1; ++merge_count) { - QueueNode current{}; + int v = -1; while (true) { - current = que.top(); + const QueueNode ¤t = que.top(); + const int current_vertex = current.vertex; + const int current_version = current.version; que.pop(); - const int v = find(current.vertex); - if (v == current.vertex && current.version == version[v]) { + v = find(current_vertex); + if (v == current_vertex && current_version == version[v]) { break; } } - const int v = current.vertex; const int p = find(up[v]); next[tail[p]] = head[v]; tail[p] = tail[v]; diff --git a/verify/yosupo-sum-of-floor-of-linear.test.cpp b/verify/yosupo-sum-of-floor-of-linear.test.cpp index 62c3a0a..d3cd2d3 100644 --- a/verify/yosupo-sum-of-floor-of-linear.test.cpp +++ b/verify/yosupo-sum-of-floor-of-linear.test.cpp @@ -1,79 +1,13 @@ // competitive-verifier: PROBLEM https://judge.yosupo.jp/problem/sum_of_floor_of_linear -#include -#include #include #include "../math/number-theory/floor-sum.hpp" -#ifdef __SIZEOF_INT128__ -static __int128 floor_div_i128(__int128 x, __int128 m) { - assert(m > 0); - __int128 q = x / m; - __int128 r = x % m; - if (r < 0) { - --q; - } - return q; -} - -static __int128 brute_floor_sum_i128(long long n, long long m, long long a, - long long b) { - assert(n >= 0); - assert(m > 0); - __int128 res = 0; - for (long long i = 0; i < n; ++i) { - res += floor_div_i128(static_cast<__int128>(a) * i + b, m); - } - return res; -} - -static std::uint64_t xorshift64(std::uint64_t &x) { - x ^= x << 7; - x ^= x >> 9; - return x; -} - -static void self_check_negative() { - std::uint64_t rng = 88172645463325252ull; - - for (int iter = 0; iter < 5000; ++iter) { - const long long n = static_cast(xorshift64(rng) % 25); - const long long m = static_cast(xorshift64(rng) % 25) + 1; - const long long a = static_cast(xorshift64(rng) % 81) - 40; - const long long b = static_cast(xorshift64(rng) % 81) - 40; - - const __int128 want = brute_floor_sum_i128(n, m, a, b); - const __int128 got1 = floor_sum<__int128>( - static_cast<__int128>(n), static_cast<__int128>(m), - static_cast<__int128>(a), static_cast<__int128>(b)); - assert(want == got1); - - const long long got2 = floor_sum(n, m, a, b); - assert(want == static_cast<__int128>(got2)); - } - - { - const long long n = 0, m = 7, a = -3, b = -5; - assert(floor_sum(n, m, a, b) == 0); - } - { - const long long n = 10, m = 1, a = -7, b = 3; - const __int128 want = brute_floor_sum_i128(n, m, a, b); - const long long got = floor_sum(n, m, a, b); - assert(want == static_cast<__int128>(got)); - } -} -#endif - int main() { std::ios::sync_with_stdio(false); std::cin.tie(nullptr); -#ifdef __SIZEOF_INT128__ - self_check_negative(); -#endif - int T; std::cin >> T; while (T--) { From 1cbc0c7231e51d7dae68a57eda59c9b3bfaf4e9d Mon Sep 17 00:00:00 2001 From: Not_Leonian <75620009+NotLeonian@users.noreply.github.com> Date: Sat, 11 Jul 2026 16:15:43 +0900 Subject: [PATCH 2/4] =?UTF-8?q?=E3=83=89=E3=82=AD=E3=83=A5=E3=83=A1?= =?UTF-8?q?=E3=83=B3=E3=83=88=E3=81=AE=E8=A1=A8=E7=8F=BE=E3=82=92=E4=B8=80?= =?UTF-8?q?=E9=83=A8=E4=BF=AE=E6=AD=A3?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../geometry/line-convex-polygon-intersection.md | 4 ++-- docs/graph/others/graph-isomorphism.md | 2 +- docs/math/matrix/dynamic-matrix-rank.md | 16 ++++++++-------- .../prime-counting-modulo.md | 4 ++-- .../generalized-floor-sum-degree-le-2.md | 4 ++-- docs/math/number-theory/generalized-garner.md | 3 ++- ...m-mod-range-increment-decrement-operations.md | 4 ++-- docs/structure/others/dynamic-median.md | 2 +- 8 files changed, 20 insertions(+), 19 deletions(-) diff --git a/docs/geometry/line-convex-polygon-intersection.md b/docs/geometry/line-convex-polygon-intersection.md index edb38e8..1bf94f9 100644 --- a/docs/geometry/line-convex-polygon-intersection.md +++ b/docs/geometry/line-convex-polygon-intersection.md @@ -68,5 +68,5 @@ documentation_of: geometry/line-convex-polygon-intersection.hpp `hull` のサイズ を $N$ とおく。 -- 時間 $O(\log N)$ -- 空間 $O(1)$ +- 時間計算量: $O(\log N)$ +- 空間計算量: $O(1)$ diff --git a/docs/graph/others/graph-isomorphism.md b/docs/graph/others/graph-isomorphism.md index 7dd19b7..5d93184 100644 --- a/docs/graph/others/graph-isomorphism.md +++ b/docs/graph/others/graph-isomorphism.md @@ -32,7 +32,7 @@ documentation_of: graph/others/graph-isomorphism.hpp また、 $i=1,2$ について、`edges_i` を $E_i$ 、 $M_i=\lvert E_i\rvert$ とする。 - $M_i$ は多重辺を重複込みで数え、自己ループも入力 $1$ 個を $1$ 本として数える。 -- $U$ を、各グラフで同じ頂点対に入った辺を $1$ つと数えたときの、 $2$ つのグラフでの個数の合計とする。自己ループの頂点対も $1$ つと数える。 +- $U$ を、各グラフで同じ頂点対に入った辺を $1$ 本と数えたときの、 $2$ つのグラフでの本数の合計とする。自己ループの頂点対も $1$ つと数える。 - 多重辺は重複を消して同一視するのではなく、各頂点対に入った辺の本数まで比較される。 - $D$ を内部隣接リストの全長の合計とする。 $D\le 2U$ である。 - $L$ を色分割中に作られる signature の最大長とする。 $L\le 2n+2$ である。 diff --git a/docs/math/matrix/dynamic-matrix-rank.md b/docs/math/matrix/dynamic-matrix-rank.md index f46f137..a423746 100644 --- a/docs/math/matrix/dynamic-matrix-rank.md +++ b/docs/math/matrix/dynamic-matrix-rank.md @@ -21,10 +21,10 @@ documentation_of: math/matrix/dynamic-matrix-rank.hpp - 空に構築する。後で `build(matrix)` を呼ぶ。 - `DynamicMatrixRank(const std::vector>& matrix)` - `matrix` で構築する。 - - 前提: `matrix` は長方形であり、`T` は体を成す。 + - 前提: `matrix` は長方形であり、`T` は体をなす。 - `void build(const std::vector>& matrix)` - `matrix` を現在の行列として前処理し直す。 - - 前提: `matrix` は長方形であり、`T` は体を成す。 + - 前提: `matrix` は長方形であり、`T` は体をなす。 - `void build()` - 現在保持している行列から前処理し直す。 - `int rank() const` @@ -39,25 +39,25 @@ documentation_of: math/matrix/dynamic-matrix-rank.hpp - 現在の行列を密行列として返す。 - `int rank_after_rank_one_update(const std::vector& column_vector, const std::vector& row_vector) const` - $A+uv^{\top}$ の階数を返す。 - - 前提: `column_vector` の長さは行数、`row_vector` の長さは列数。 + - 前提: `column_vector` の長さは行数に等しく、`row_vector` の長さは列数に等しい。 - 備考: 内部状態は変更しない。 - `int rank_after_row_replacement(int row_index, const std::vector& new_row) const` - `row_index` 行目を `new_row` に差し替えた行列の階数を返す。 - - 前提: $0\le i& new_column) const` - `column_index` 列目を `new_column` に差し替えた行列の階数を返す。 - - 前提: $0\le j& column_vector, const std::vector& row_vector)` - $A+uv^{\top}$ に内部状態を更新し、その階数を返す。 - - 前提: `column_vector` の長さは行数、`row_vector` の長さは列数。 + - 前提: `column_vector` の長さは行数に等しく、`row_vector` の長さは列数に等しい。 - `int apply_row_replacement(int row_index, const std::vector& new_row)` - `row_index` 行目を `new_row` に差し替え、変更後の階数を返す。 - - 前提: $0\le i& new_column)` - `column_index` 列目を `new_column` に差し替え、変更後の階数を返す。 - - 前提: $0\le j Date: Sat, 11 Jul 2026 16:30:25 +0900 Subject: [PATCH 3/4] =?UTF-8?q?=E3=83=98=E3=83=83=E3=83=80=E3=83=95?= =?UTF-8?q?=E3=82=A1=E3=82=A4=E3=83=AB=E3=81=AE=E3=82=B3=E3=83=A1=E3=83=B3?= =?UTF-8?q?=E3=83=88=E3=82=92=E4=B8=80=E9=83=A8=E4=BF=AE=E6=AD=A3?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- docs/math/number-theory/generalized-garner.md | 2 +- internal/int128.hpp | 7 ++++--- math/number-theory/generalized-garner.hpp | 8 ++++---- other/enumerate-maximum-independent-set-path-sums.hpp | 2 +- other/rectangle-add-max-get.hpp | 4 ++-- 5 files changed, 12 insertions(+), 11 deletions(-) diff --git a/docs/math/number-theory/generalized-garner.md b/docs/math/number-theory/generalized-garner.md index 3213a21..44b0652 100644 --- a/docs/math/number-theory/generalized-garner.md +++ b/docs/math/number-theory/generalized-garner.md @@ -32,5 +32,5 @@ documentation_of: math/number-theory/generalized-garner.hpp 入力する式の個数を $N$ 、全ての法と $2$ の最大値を $V$ とする。 - 標準の 64 bit 以下の整数型では時間 $O(N\log V)$ 。 -- それ以外の整数型で通常乗算を使えない場合は時間 $O(N\log^2 V)$ 。 +- それ以外の整数型で通常の乗算を使えない場合は時間 $O(N\log^2 V)$ 。 - 空間 $O(1)$ 。 diff --git a/internal/int128.hpp b/internal/int128.hpp index c90e898..3a9644c 100644 --- a/internal/int128.hpp +++ b/internal/int128.hpp @@ -182,7 +182,7 @@ class UInt128 { return length; } - // divisor != 0 が呼び出し元で保証されていることを仮定する。 + // divisor ≠ 0 が呼び出し元で保証されていることを仮定する。 static constexpr UInt128 div_mod_uint32(UInt128 value, std::uint32_t divisor, std::uint32_t &remainder) { @@ -365,7 +365,7 @@ class UInt128 { remainder = from_32bit_words(remainder_words); } - // rhs != 0 が呼び出し元で保証されていることを仮定する。 + // rhs ≠ 0 が呼び出し元で保証されていることを仮定する。 static constexpr void div_mod_unchecked(UInt128 lhs, UInt128 rhs, UInt128 "ient, UInt128 &remainder) { @@ -589,7 +589,8 @@ class Int128 { return from_unsigned_unchecked(value, negative); } - // 符号を付けた value が Int128 の範囲に収まることを仮定する。 + // 引数 negative で指定した通りに符号を付けた value が + // Int128 の範囲に収まることを仮定する。 static constexpr Int128 from_unsigned_unchecked(UInt128 value, bool negative) { return from_twos_complement(negative ? -value : value); diff --git a/math/number-theory/generalized-garner.hpp b/math/number-theory/generalized-garner.hpp index 40bd233..813321a 100644 --- a/math/number-theory/generalized-garner.hpp +++ b/math/number-theory/generalized-garner.hpp @@ -7,8 +7,8 @@ // 係数や法は互いに素でなくてよい。解がなければ (0, 0) を返す。 // 標準の 64 bit 以下の整数型では剰余乗算に 128 bit 整数型を用いる。 // それ以外の整数型では加算と 2 倍による剰余乗算を用いる。 -// 式の個数を N、2 と全ての法の最大値を V とすると、計算量 O(N log V)。 -// ただし、通常乗算を使えない型では O(N log^2 V)。 +// 式の個数を N、全ての法と 2 の最大値を V とすると、計算量 O(N log V)。 +// ただし、通常の乗算を使えない型では O(N log^2 V)。 #include #include @@ -117,7 +117,7 @@ template T inv_mod(T a, T m) { const T c = a - q * b; a = b; b = c; - // 初期値と呼び出し元が保証する互いに素という条件から 0 < q < m。 + // 初期値と呼び出し元が互いに素を保証することから 0 < q < m。 const T qy = mul_mod_normalized(q, y, m); const T z = x >= qy ? x - qy : x + (m - qy); x = y; @@ -168,7 +168,7 @@ template std::pair solve_linear_congruence(T a, T b, T m) { if (a == 1) { return {b, m}; } - // 0 < a < m なので、引数順をこの向きにすると最初の剰余計算を省ける。 + // 0 < a < m なので、引数をこの順序にすると最初の剰余の計算を省ける。 const T g = gcd(m, a); const T b_div_g = b / g; if (b_div_g * g != b) { diff --git a/other/enumerate-maximum-independent-set-path-sums.hpp b/other/enumerate-maximum-independent-set-path-sums.hpp index 05b5da9..e986735 100644 --- a/other/enumerate-maximum-independent-set-path-sums.hpp +++ b/other/enumerate-maximum-independent-set-path-sums.hpp @@ -63,7 +63,7 @@ std::vector marginal_values(const std::vector &a) { return res; } -// 呼び出し元で x の絶対値が offset 以下であることを保証する。 +// x の絶対値が offset 以下であることを呼び出し元で保証する。 template std::size_t bucket_index(const T x, const std::size_t offset) { if (T(0) <= x) { diff --git a/other/rectangle-add-max-get.hpp b/other/rectangle-add-max-get.hpp index 2f9761d..7bec6a2 100644 --- a/other/rectangle-add-max-get.hpp +++ b/other/rectangle-add-max-get.hpp @@ -367,7 +367,7 @@ template struct CompressedRectangleAddMaxGet { } found = true; } else if (ret.max_value == now.max_value) { - // x 座標を昇順に走査しているため、最小点は更新せず、最大点は常に更新する。 + // x 座標を昇順に走査しているため、最小の点は更新せず、最大の点は常に更新する。 ret.maximum_x = maximum_x; ret.maximum_y = now.maximum_y; if constexpr (NeedArea) { @@ -646,7 +646,7 @@ template struct RectangleAddMaxGet { } found = true; } else if (ret.max_value == value) { - // x 座標を昇順に走査しているため、最小点は更新せず、最大点は常に更新する。 + // x 座標を昇順に走査しているため、最小の点は更新せず、最大の点は常に更新する。 ret.maximum_x = maximum_x; ret.maximum_y = maximum_y; if constexpr (NeedArea) { From a8bbd39c74f008739bf8225cd0673974c2900ed0 Mon Sep 17 00:00:00 2001 From: Not_Leonian <75620009+NotLeonian@users.noreply.github.com> Date: Sat, 11 Jul 2026 16:36:58 +0900 Subject: [PATCH 4/4] =?UTF-8?q?=E3=83=98=E3=83=83=E3=83=80=E3=83=95?= =?UTF-8?q?=E3=82=A1=E3=82=A4=E3=83=AB=E3=81=AE=E3=82=B3=E3=83=A1=E3=83=B3?= =?UTF-8?q?=E3=83=88=E3=82=92=E3=81=95=E3=82=89=E3=81=AB=E4=B8=80=E9=83=A8?= =?UTF-8?q?=E4=BF=AE=E6=AD=A3?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- math/matrix/dynamic-matrix-rank.hpp | 2 +- structure/others/dynamic-median.hpp | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/math/matrix/dynamic-matrix-rank.hpp b/math/matrix/dynamic-matrix-rank.hpp index 7e84d22..9820cef 100644 --- a/math/matrix/dynamic-matrix-rank.hpp +++ b/math/matrix/dynamic-matrix-rank.hpp @@ -7,7 +7,7 @@ // 1 行差し替え、1 列差し替えはそれぞれ e_i (b-a_i)^T, (b-a_j) e_j^T に帰着する。 // 前処理は O(r + rc min(r, c) + k^2(r + c))、更新や判定は // O((k + 1)(r + c)) である。 -// T は体を成し、零判定と四則演算ができることを仮定する。 +// T は体をなし、零判定と四則演算ができることを仮定する。 #include #include diff --git a/structure/others/dynamic-median.hpp b/structure/others/dynamic-median.hpp index a9736f6..b613155 100644 --- a/structure/others/dynamic-median.hpp +++ b/structure/others/dynamic-median.hpp @@ -81,7 +81,8 @@ template struct DynamicMedian { private: void balance() { - // add または erase の直前に不変条件が成り立つため、移動は高々 1 回でよい。 + // add または erase でこの関数が呼び出される位置と + // この関数を呼ばない関数についての不変条件より、移動は高々 1 回でよい。 if (lower_values.size() < upper_values.size()) { lower_values.insert(upper_values.extract(upper_values.begin())); } else if (lower_values.size() > upper_values.size() + 1) {