Sync dev

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

template/base_template/immintrin.hpp

@@ -0,0 +1,3 @@
+#pragma once
+//  PRIORITY 1
+#include <immintrin.h>

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

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,62 +10,76 @@
 // Verified by AOJ GRL_6_A http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=GRL_6_A
 template <typename Capacity = long long>
 struct dinic {
-    using edge = tuple<size_t, Capacity, size_t>; // {to, cap, rev idx at g[to]}
-    size_t n;
-    vector<vector<edge>> g;
-    vector<size_t> 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<vector<edge>> g = {};
+    lazy_fill_vector<size_t> level = {};
+    lazy_fill_vector<size_t> 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<size_t>::max();
-        fill(ALL(level), ZINF);
-        deque<size_t> q = {s};
+
+    void bfs(const size_t s) {
+        constexpr auto ZINF = numeric_limits<size_t>::max();
+        level.fill(ZINF);
+        auto q = deque<size_t>{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);
                 }
             }
         }
     }
-    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<size_t>::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<Capacity>::max())) > 0) {
+            iter.fill(0);
+            for (
+                auto d = Capacity{};
+                (d = dfs(s, t, numeric_limits<Capacity>::max())) > 0;
+            ) {
                 f += d;
             }
         }

template/flow/ford_fulkerson.hpp

@@ -1,49 +1,60 @@
 #pragma once
 #include "base_template.hpp"
+#include "array/lazy_fill_vector.hpp"
 // Ford-Fulkerson法による最大流
 // Capacityには必ず整数型を指定すること
 // Capacityが整数ならば，停止しないときDFS1回O(E)で必ず1以上フローが新たに流れるためO(FE)
 // Verified by AtC ABC010-D https://atcoder.jp/contests/abc010/tasks/abc010_4
 // Verified by AOJ GRL_6_A http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=GRL_6_A
 template <typename Capacity = long long>
 struct ford_fulkerson {
-    using edge = tuple<size_t, Capacity, size_t>; // {to, cap, rev idx at g[to]}
-    size_t n;
-    vector<vector<edge>> g;
-    vector<uint_fast8_t> used;
+    struct edge {
+        size_t to = {};
+        Capacity cap = {};
+        size_t rev = {};
+    };
+
+    size_t n = {};
+    vector<vector<edge>> g = {};
+    lazy_fill_vector<uint_fast8_t> 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<Capacity>::max());
+            used.fill(false);
+            const auto d = dfs(s, t, numeric_limits<Capacity>::max());
             if (d == 0) return f;
             f += d;
         }

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 <uint64_t n>
+constexpr uint64_t montgomery_reduction(const uint64_t t) {
+    constexpr auto r = uint64_t{1} << 32;
+    constexpr auto ninv = static_cast<uint64_t>((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 <uint64_t n>
+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<n>(
+        montgomery_reduction<n>(a * r2mn) * montgomery_reduction<n>(b * r2mn)
+    );
+}

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<modint> fac_cache; \