diff --git a/bin/run.py b/bin/run.py index 3edcacd..629f76d 100644 --- a/bin/run.py +++ b/bin/run.py @@ -36,16 +36,17 @@ def compile(source_path): debug_print(f"Compiling {source_path} ...") start_time = time.time() compilation = subprocess.run([ - "g++-12", - "-std=gnu++2b", + "g++-15", + "-std=gnu++23", "-O2", "-Wall", "-Wextra", - "-mtune=native", "-march=native", - "-fconstexpr-depth=2147483647", - "-fconstexpr-loop-limit=2147483647", - "-fconstexpr-ops-limit=2147483647", + "-flto=auto", + "-fmodules", + "-fconstexpr-depth=1024", + "-fconstexpr-loop-limit=524288", + "-fconstexpr-ops-limit=2097152", "-DYSN_DEBUG", source_path, "-o", exec_path diff --git a/template/base_template/immintrin.hpp b/template/base_template/immintrin.hpp new file mode 100644 index 0000000..c128b82 --- /dev/null +++ b/template/base_template/immintrin.hpp @@ -0,0 +1,3 @@ +#pragma once +//  PRIORITY 1 +#include \ No newline at end of file diff --git a/template/conv/butterfly_simd.hpp b/template/conv/butterfly_simd.hpp new file mode 100644 index 0000000..22111dc --- /dev/null +++ b/template/conv/butterfly_simd.hpp @@ -0,0 +1,94 @@ +#pragma once +#include "base_template.hpp" +#include "base_template/immintrin.hpp" +#include "modint/modint_petit_p.hpp" + +constexpr uint32_t MOD = 998244353; +constexpr uint64_t R2 = (uint64_t(1) << 32) % MOD * ((uint64_t(1) << 32) % MOD) % MOD; +// INV = -MOD^{-1} mod 2^32 を事前に計算 +constexpr uint32_t INV32 = modinv_extgcd(MOD, uint64_t{1} << 32); + +// ——— モンゴメリ乗算 + 還元 ——— +// 64bit レーン ×4（計 256bit）版 +static inline __m256i montgomery_mul256(__m256i a, __m256i b) { + // acc = a * b の 52bit 積を下位/上位に分けて累積 + __m256i acc = _mm256_setzero_si256(); + acc = _mm256_madd52lo_epu64(acc, a, b); // VPMADD52LUQ + acc = _mm256_madd52hi_epu64(acc, a, b); // VPMADD52HUQ + // acc = z + m*MOD (両命令で既に加算済み) → 64bit 単位で u + // 最後に 2^32 でシフトすれば Montgomery 還元完了 + return _mm256_srli_epi64(acc, 32); +} + +// 8 要素一括バタフライ（32bit 値を 64bit レーンにゼロ拡張して処理） +void butterfly_simd( + uint32_t *a, // モンゴメリ表現済みの配列 + long long n, // 全長 + long long p, // ステップ + uint32_t irot, // ω^p のモンゴメリ表現 + uint32_t iimag // ω^{2p} のモンゴメリ表現 +) { + // 64bit レーンに直接セット + __m256i vIrot = _mm256_set1_epi64x(irot); + __m256i vImag = _mm256_set1_epi64x(iimag); + __m256i vIrot2 = _mm256_set1_epi64x(uint64_t(irot)*irot % MOD); + __m256i vIrot3 = _mm256_set1_epi64x(uint64_t(irot)*irot % MOD * irot % MOD); + __m256i vR2 = _mm256_set1_epi64x(R2); + + for ( long long i = 0; i < n; i += 4*p ) { + for ( long long j = 0; j < p; j += 4 ) { + // 32bit→64bit ゼロ拡張ロード + __m256i a0 = _mm256_cvtepu32_epi64(_mm_loadu_si128((__m128i*)(a + i + j ))); + __m256i a1 = _mm256_cvtepu32_epi64(_mm_loadu_si128((__m128i*)(a + i + j + p ))); + __m256i a2 = _mm256_cvtepu32_epi64(_mm_loadu_si128((__m128i*)(a + i + j + 2*p))); + __m256i a3 = _mm256_cvtepu32_epi64(_mm_loadu_si128((__m128i*)(a + i + j + 3*p))); + + // t0 = (a0 + a1 + a2 + a3) * R^{-1} mod MOD + __m256i sum01 = _mm256_add_epi64(a0, a1); + __m256i sum23 = _mm256_add_epi64(a2, a3); + __m256i t0 = _mm256_add_epi64(sum01, sum23); + t0 = montgomery_mul256(vR2, t0); + + // t1 = (a0 - a1 + (a2 - a3)*iimag) * irot + __m256i sub01 = _mm256_sub_epi64(a0, a1); + __m256i sub23 = _mm256_sub_epi64(a2, a3); + __m256i tmp1 = montgomery_mul256(sub23, vImag); + __m256i t1_pre = _mm256_add_epi64(sub01, tmp1); + __m256i t1 = montgomery_mul256(t1_pre, vIrot); + + // t2 = (a0 + a1 - a2 - a3) * irot2 + __m256i t2_pre = _mm256_sub_epi64(sum01, sum23); + __m256i t2 = montgomery_mul256(t2_pre, vIrot2); + + // t3 = (a0 - a1 - (a2 - a3)*iimag) * irot3 + __m256i t3_pre = _mm256_sub_epi64(sub01, tmp1); + __m256i t3 = montgomery_mul256(t3_pre, vIrot3); + + // 64bit→32bit に切り捨て（値範囲は <2^32 なので安全） + __m128i r0 = _mm256_cvtepi64_epi32(t0); + __m128i r1 = _mm256_cvtepi64_epi32(t1); + __m128i r2 = _mm256_cvtepi64_epi32(t2); + __m128i r3 = _mm256_cvtepi64_epi32(t3); + + // 書き戻し + _mm_storeu_si128((__m128i*)(a + i + j ), r0); + _mm_storeu_si128((__m128i*)(a + i + j + p), r1); + _mm_storeu_si128((__m128i*)(a + i + j +2*p), r2); + _mm_storeu_si128((__m128i*)(a + i + j +3*p), r3); + } + } +} + +void to_montgomery_simd(uint32_t *a, size_t n) { + __m256i vR2 = _mm256_set1_epi64x(R2); + for ( size_t i = 0; i < n; i += 4 ) { + // ① 32bit×4 を 64bit×4 レーンにゼロ拡張ロード + __m128i tmp = _mm_loadu_si128((__m128i*)(a + i)); + __m256i v = _mm256_cvtepu32_epi64(tmp); + // ② 各要素を montgomery_mul256(v, vR2) → x * R mod M + v = montgomery_mul256(v, vR2); + // ③ 64bit→32bit に安全にトランケートして戻す + __m128i out = _mm256_cvtepi64_epi32(v); + _mm_storeu_si128((__m128i*)(a + i), out); + } +} \ No newline at end of file diff --git a/template/flow/dinic.hpp b/template/flow/dinic.hpp index 1041ec0..9352770 100644 --- a/template/flow/dinic.hpp +++ b/template/flow/dinic.hpp @@ -1,5 +1,6 @@ #pragma once #include "base_template.hpp" +#include "array/lazy_fill_vector.hpp" // Dinic法による最大流 // Capacityには必ず整数型を指定すること // BFSとDFSを組み合わせてO(EV^2)だが実用的にはもっと速いらしい @@ -9,30 +10,40 @@ // Verified by AOJ GRL_6_A http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=GRL_6_A template struct dinic { - using edge = tuple; // {to, cap, rev idx at g[to]} - size_t n; - vector> g; - vector level, iter; - dinic() : n(0), g(), level(), iter() {} - dinic(size_t _n) : n(_n), g(n), level(n), iter(n) {} - void add_arc(size_t from, size_t to, Capacity cap) { + struct edge { + size_t to = {}; + Capacity cap = {}; + size_t rev = {}; + }; + + size_t n = {}; + vector> g = {}; + lazy_fill_vector level = {}; + lazy_fill_vector iter = {}; + + dinic() = default; + dinic(const size_t _n) : n(_n), g(n), level(n), iter(n) {} + + void add_arc(const size_t from, const size_t to, const Capacity cap) { assert(from < n && to < n); assert(cap >= 0); - g[from].emplace_back(to , cap, g[to ].size()); - g[to ].emplace_back(from, 0, g[from].size() - 1); + g[from].emplace_back(edge{.to = to, .cap = cap, .rev = g[to].size()}); + g[to].emplace_back(edge{.to = from, .cap = 0, .rev = g[from].size() - 1}); } - void add_edge(size_t u, size_t v, Capacity cap) { + + void add_edge(const size_t u, const size_t v, const Capacity cap) { add_arc(u, v, cap); add_arc(v, u, cap); } - void bfs(size_t s) { - constexpr size_t ZINF = numeric_limits::max(); - fill(ALL(level), ZINF); - deque q = {s}; + + void bfs(const size_t s) { + constexpr auto ZINF = numeric_limits::max(); + level.fill(ZINF); + auto q = deque{s}; level[s] = 0; while (q.size()) { - size_t u = q.front(); q.pop_front(); - for (const auto &[v, cap, _] : g[u]) { + const auto u = q.front(); q.pop_front(); + for (const auto& [v, cap, rev] : g[u]) { if (cap > 0 && level[v] == ZINF) { level[v] = level[u] + 1; q.push_back(v); @@ -40,31 +51,35 @@ struct dinic { } } } - Capacity dfs(size_t v, size_t t, Capacity f) { + + Capacity dfs(const size_t v, const size_t t, const Capacity f) { if (v == t) return f; - for (size_t &i = iter[v]; i < g[v].size(); i++) { + for (auto& i = iter[v]; i < g[v].size(); i++) { auto& [to, cap, rev] = g[v][i]; if (cap <= 0 || level[v] >= level[to]) continue; - Capacity d = dfs(to, t, min(f, cap)); + const auto d = dfs(to, t, min(f, cap)); if (d > 0) { cap -= d; - get<1>(g[to][rev]) += d; + g[to][rev].cap += d; return d; } } return 0; } + Capacity max_flow(size_t s, size_t t) { assert(s < n && t < n); constexpr size_t ZINF = numeric_limits::max(); if (s == t) return 0; - Capacity f = 0; + auto f = Capacity{0}; while (true) { bfs(s); if (level[t] == ZINF) return f; - fill(ALL(iter), 0); - Capacity d; - while ((d = dfs(s, t, numeric_limits::max())) > 0) { + iter.fill(0); + for ( + auto d = Capacity{}; + (d = dfs(s, t, numeric_limits::max())) > 0; + ) { f += d; } } diff --git a/template/flow/ford_fulkerson.hpp b/template/flow/ford_fulkerson.hpp index f2be3c0..6c538f2 100644 --- a/template/flow/ford_fulkerson.hpp +++ b/template/flow/ford_fulkerson.hpp @@ -1,5 +1,6 @@ #pragma once #include "base_template.hpp" +#include "array/lazy_fill_vector.hpp" // Ford-Fulkerson法による最大流 // Capacityには必ず整数型を指定すること // Capacityが整数ならば，停止しないときDFS1回O(E)で必ず1以上フローが新たに流れるためO(FE) @@ -7,43 +8,53 @@ // Verified by AOJ GRL_6_A http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=GRL_6_A template struct ford_fulkerson { - using edge = tuple; // {to, cap, rev idx at g[to]} - size_t n; - vector> g; - vector used; + struct edge { + size_t to = {}; + Capacity cap = {}; + size_t rev = {}; + }; + + size_t n = {}; + vector> g = {}; + lazy_fill_vector used = {}; + ford_fulkerson() : n(0), g(), used() {} - ford_fulkerson(size_t _n) : n(_n), g(n), used(n) {} - void add_arc(size_t from, size_t to, Capacity cap) { - assert(0 <= from && from < n && 0 <= to && to < n); + ford_fulkerson(const size_t _n) : n(_n), g(n), used(n) {} + + void add_arc(const size_t from, const size_t to, const Capacity cap) { + assert(from < n && to < n); assert(cap >= 0); - g[from].emplace_back(to , cap, g[to ].size()); - g[to ].emplace_back(from, 0, g[from].size() - 1); + g[from].emplace_back(edge{.to = to, .cap = cap, .rev = g[to].size()}); + g[to].emplace_back(edge{.to = from, .cap = 0, .rev = g[from].size() - 1}); } - void add_edge(size_t u, size_t v, Capacity cap) { + + void add_edge(const size_t u, const size_t v, const Capacity cap) { add_arc(u, v, cap); add_arc(v, u, cap); } - Capacity dfs(size_t v, size_t t, Capacity f) { + + Capacity dfs(const size_t v, const size_t t, const Capacity f) { if (v == t) return f; used[v] = true; - for (auto &[to, cap, rev] : g[v]) { + for (auto& [to, cap, rev] : g[v]) { if (used[to] || cap <= 0) continue; - Capacity d = dfs(to, t, min(f, cap)); + const auto d = dfs(to, t, min(f, cap)); if (d > 0) { cap -= d; - get<1>(g[to][rev]) += d; + g[to][rev].cap += d; return d; } } return 0; } - Capacity max_flow(size_t s, size_t t) { - assert(0 <= s && s < n && 0 <= t && t < n); + + Capacity max_flow(const size_t s, const size_t t) { + assert(s < n && t < n); if (s == t) return 0; - Capacity f = 0; + auto f = Capacity{0}; while (true) { - fill(ALL(used), false); - Capacity d = dfs(s, t, numeric_limits::max()); + used.fill(false); + const auto d = dfs(s, t, numeric_limits::max()); if (d == 0) return f; f += d; } diff --git a/template/modint/montgomery.hpp b/template/modint/montgomery.hpp new file mode 100644 index 0000000..51767dd --- /dev/null +++ b/template/modint/montgomery.hpp @@ -0,0 +1,24 @@ +#pragma once +#include "base_template.hpp" +#include "modint/modint_petit_p.hpp" + +// t * n^{-1} mod 2^32 +template +constexpr uint64_t montgomery_reduction(const uint64_t t) { + constexpr auto r = uint64_t{1} << 32; + constexpr auto ninv = static_cast((r - modinv_extgcd(n, r)) % (r - 1)); + const auto m = ((t & (r - 1)) * ninv) & (r - 1); + const auto t2 = (t + m * n) >> 32; + return t2 - (t2 >= n ? n : 0); +} + +// a * b * (2^32) mod N を % なしで求める +template +constexpr uint64_t montgomery_prod(const uint64_t a, const uint64_t b) { + constexpr auto r = uint64_t{1} << 32; + constexpr auto rmn = r % n; + constexpr auto r2mn = rmn * rmn % n; + return montgomery_reduction( + montgomery_reduction(a * r2mn) * montgomery_reduction(b * r2mn) + ); +} \ No newline at end of file diff --git a/template/modint/using_modint.hpp b/template/modint/using_modint.hpp index 75a727b..3009705 100644 --- a/template/modint/using_modint.hpp +++ b/template/modint/using_modint.hpp @@ -9,7 +9,7 @@ using modint = moduloint<(PDIV)>; \ \ /* Literal */ \ - constexpr modint operator"" _p(unsigned long long _x) noexcept { return modint(_x); } \ + constexpr modint operator""_p(unsigned long long _x) noexcept { return modint(_x); } \ \ /* Factorial n! */ \ factorial_cache fac_cache; \ diff --git a/template/segment_union.hpp b/template/segment_union.hpp index 30d50a1..bd44a2e 100644 --- a/template/segment_union.hpp +++ b/template/segment_union.hpp @@ -62,7 +62,7 @@ struct segment_union { // 区間 [l, r) を追加; S <- S ∪ [l, r) // 区間数に対してならし対数時間 // 削除した区間を返す - vector insert(const value_type l, const value_type r) { + pair> insert(const value_type l, const value_type r) { auto added_segment = pair{l, r}; auto removed_segments = vector{}; @@ -91,24 +91,25 @@ struct segment_union { s.erase(it); } - s.emplace(move(added_segment)); + s.emplace(added_segment); - return removed_segments; + return {move(added_segment), move(removed_segments)}; } // 区間 [l, r) を意味する seg = {l, r} を追加; S <- S ∪ [l, r) // 区間数に対してならし対数時間 template - void insert(const Seg& seg) { + vector insert(const Seg& seg) { const auto& [l, r] = seg; - insert(l, r); + return insert(l, r); } // 区間 [l, r) を除去; S <- S \ [l, r) // 区間数に対してならし対数時間 - void erase(const value_type l, const value_type r) { + pair, vector> erase(const value_type l, const value_type r) { auto replaced_segments = vector{}; + auto emplaced_segments = vector{}; const auto seg = segment_type{l, r}; // 左 @@ -131,9 +132,17 @@ struct segment_union { // 除去して {l,r} を除いた区間に変換 for (const auto& [x, y] : replaced_segments) { - if (x < l) s.emplace(x, l); - if (r < y) s.emplace(r, y); + if (x < l) { + s.emplace(x, l); + emplaced_segments.emplace_back(x, l); + } + if (r < y) { + s.emplace(r, y); + emplaced_segments.emplace_back(r, y); + } } + + return {move(replaced_segments), move(emplaced_segments)}; } diff --git a/template/segtree/range/linear.hpp b/template/segtree/range/linear.hpp index d7a1d42..4780b13 100644 --- a/template/segtree/range/linear.hpp +++ b/template/segtree/range/linear.hpp @@ -49,7 +49,7 @@ struct segtree_linear { query_type h; REP(i, n * n) h[i] = 0; REP(i, n) REP(k, n) RPE(j, n) { - h[i * n + j] = f[i * n + k] * g[k * n + j]; + h[i * n + j] += f[i * n + k] * g[k * n + j]; } return h; }