Develop feb11

bin/run.py

@@ -44,6 +44,10 @@ def compile(source_path):
         "-march=native",
         "-flto=auto",
         "-fmodules",
+        "-fopenmp",
+        "-pthread",
+        "-lstdc++exp",
+        "-ftrivial-auto-var-init=zero",
         "-fconstexpr-depth=1024",
         "-fconstexpr-loop-limit=524288",
         "-fconstexpr-ops-limit=2097152",

template/array/sahz.hpp

@@ -1,13 +1,37 @@
 #pragma once
 #include "base_template.hpp"
-// 座圧 {写像, 逆写像} が返る O(NlogN)
-// verified by https://atcoder.jp/contests/abc036/tasks/abc036_c
-template <typename T, typename U = size_t, typename M = unordered_map<T, U>>
-pair<M, vector<T>> sahz(vector<T> x) {
-    UNIQUE(x);
-    M r;
-    REP(i, x.size()) {
-        r[x[i]] = i;
+#include "ranges_to.hpp"
+
+// 座圧 {写像, 逆写像} が返る
+// Map が std::unordered_map の場合、平均 O(N) 最悪 O(N^2)
+// Map が std::map の場合、平均 O(NlogN) 最悪 O(NlogN)
+// verified at <<WIP>>
+template <
+    typename R,
+    typename Map = unordered_map<ranges::range_value_t<R>, ptrdiff_t>,
+    typename Comp = less<>,
+    typename Proj = identity
+>
+auto sahz(R&& r, Comp comp = {}, Proj proj = {})
+    -> pair<Map, vector<ranges::range_value_t<R>>>
+{
+    using T = ranges::range_value_t<R>;
+    auto unique_values = r | to_vector<T>{};
+    {
+        ranges::sort(unique_values, comp, proj);
+        auto equiv_comp = [&](const auto& x, const auto& y) -> bool {
+            return !comp(x, y) && !comp(y, x);
+        };
+        const auto erased_range = ranges::unique(unique_values, equiv_comp, proj);
+        unique_values.erase(erased_range.begin(), erased_range.end());
     }
-    return {r, x};
-}
+
+    auto smap = Map{};
+    auto index = ptrdiff_t{0};
+    for (auto&& x : unique_values) {
+        smap[x] = index++;
+    }
+
+    return pair{move(smap), move(unique_values)};
+
+}

template/conv/ntt.hpp

@@ -11,9 +11,11 @@
 // Cooley-Tukey型 高速フーリエ変換 O(NlogN)
 // 原始根 prim_ が指定されない場合は最小の原始根が指定される
 template <typename T, ll prim_ = -1>
-requires mod_integral<typename T::value_type>
-constexpr auto ntt(T&& a, const bool inv = false) -> vector<typename T::value_type> {
-    using V = typename T::value_type;
+requires mod_integral<typename std::remove_reference_t<T>::value_type>
+constexpr auto ntt(T&& a, const bool inv = false)
+    -> vector<typename std::remove_reference_t<T>::value_type>
+{
+    using V = typename std::remove_reference_t<T>::value_type;
     constexpr auto prim = V{prim_ != -1 ? prim_ : primitive_root(V::pdiv)};
     if (inv) {
         return butterfly_inv<V, prim>(forward<T>(a));

template/fps/berlekamp_massey.hpp

@@ -0,0 +1,39 @@
+#pragma once
+#include "base_template.hpp"
+#include "modint/modint.hpp"
+
+// Berlekamp-Massey O(N^2)
+// C-recursive の母関数を FPS 有理式で表したときの分母を返す
+// original: https://nyaannyaan.github.io/library/fps/berlekamp-massey.hpp
+template <mod_integral T>
+vector<T> berlekamp_massey(const vector<T>& s) {
+    const auto n = ssize(s);
+    auto b = vector{T{1}};
+    auto c = vector{T{1}};
+    b.reserve(n + 1);
+    c.reserve(n + 1);
+    auto y = T{1};
+    for (auto ed = 1LL; ed <= n; ed++) {
+        auto x = T{0};
+        assert(ed >= ssize(c));
+        RPE(i, ssize(c)) {
+            const auto offset = ed - ssize(c);
+            x += c[i] * s[offset + i];
+        }
+        b.emplace_back(T{0});
+        if (x == 0) continue;
+        const auto freq = x / y;
+        if (ssize(c) < ssize(b)) {
+            auto cprev = c;
+            c.insert(c.begin(), ssize(b) - ssize(c), T{0});
+            assert(ssize(b) == ssize(c));
+            REP(i, ssize(c)) c.rbegin()[i] -= freq * b.rbegin()[i];
+            b = move(cprev);
+            y = x;
+        } else {
+            REP(i, ssize(b)) c.rbegin()[i] -= freq * b.rbegin()[i];
+        }
+    }
+    ranges::reverse(c);
+    return c;
+}

template/fps/bostan_mori.hpp

@@ -4,18 +4,14 @@
 #include "conv/ntt.hpp"
 
 // Bostan-Mori [x^m] f(x) / g(x) を求める O(NlogNlogM) (N: deg g)
-// verified at https://judge.yosupo.jp/submission/336936
-// verified at https://atcoder.jp/contests/abc436/submissions/71731332
+// verified at https://judge.yosupo.jp/submission/339373
+// verified at https://atcoder.jp/contests/abc436/submissions/71918989
 template <mod_integral T>
 constexpr T bostan_mori(vector<T> f, vector<T> g, ll m) {
     assert(m >= 0);
     assert(!g.empty());
     assert(g[0] != 0);
     assert(ssize(f) < ssize(g));
-    if (m == 0) {
-        if (f.empty()) return 0;
-        return f[0] / g[0];
-    }
     // P(x)Q(-x) = A(x^2) + x B(x^2),
     // Q(x)Q(-x) = R(x^2) ∵ Q(x)Q(-x) は偶
     // 以上のように分解すると
@@ -28,25 +24,113 @@ constexpr T bostan_mori(vector<T> f, vector<T> g, ll m) {
 
     // ステップごとにmが半減するのでダブリングの要領で O(NlogNlogM)
     // P,Qの次数は畳み込んで半減なので±1を除いて変わらず
+
+    // 先にNTT用の長さにしておく
+    const auto fg_size = (ll)bit_ceil<ull>(max(ssize(f), ssize(g)) * 2 - 1);
+    f.resize(fg_size, T{0});
+    g.resize(fg_size, T{0});
 
-    // 分母
-    auto gnx = g;  // g(-x)
-    STEP(i, 1, ssize(gnx) - 1, 2) {
-        gnx[i] *= -1;
-    }
-    auto rx2 = convolve_p(move(g), gnx);
-    auto r = vector<T>{};
-    r.reserve((rx2.size() + 1) / 2);
-    STEP(i, 0, ssize(rx2) - 1, 2) {
-        r.emplace_back(rx2[i]);
+    // 周波数領域ずらし用前計算 (ω_{2n}^BitReverse(i))_i=0^{n-1}
+    auto ntt_coeff_rev = vector(fg_size / 2, T{0});
+    auto ntt_coeff_inv = vector(fg_size / 2, T{0});
+    {
+        constexpr auto prim = T{primitive_root(T::pdiv)};
+        const auto root = prim.pow((T::pdiv - 1) / fg_size);
+        // const auto root_inv = prim.pow((T::pdiv - 1) / fg_size).inv();
+        const auto root_inv = root.inv();
+        auto coeff = T{1};
+        auto coeff_inv = T{1};
+        ntt_coeff_rev[0] = coeff;
+        ntt_coeff_inv[0] = coeff_inv;
+        auto idx = 0LL;
+        REP(i, fg_size / 2 - 1) {
+            // bit-reverse
+            auto bit = ~i & (i + 1);
+            auto rev = (fg_size / 4) / bit;
+            idx ^= (fg_size / 2 - 1) & ~(rev - 1);
+            // non-one coeff
+            ntt_coeff_rev[i + 1] = (coeff *= root);  // ⚠️ NOT bit reversed!
+            ntt_coeff_inv[idx] = (coeff_inv *= root_inv);
+        }
     }
 
-    auto ab = convolve_p(move(f), move(gnx));
-    auto ab2 = vector<T>{};
-    ab2.reserve((ssize(ab) + (m + 1) % 2) / 2);
-    STEP(i, m % 2, ssize(ab) - 1, 2) {
-        ab2.emplace_back(ab[i]);
-    }
+    // 周波数領域に写しておく
+    f = ntt(move(f));
+    g = ntt(move(g));
+
+    auto recurse = [fg_size, &ntt_coeff_rev, &ntt_coeff_inv](
+        const auto recurse, vector<T> f, vector<T> g, ll m
+    ) -> T {
+
+        if (m == 0) {
+            // if (f.empty()) return 0;
+            // return f[0] / g[0];
+            return accumulate(ALL(f), T{0}) / accumulate(ALL(g), T{0});
+        }
+
+        // g(-x)
+        auto gnx = vector(fg_size, T{0});
+        REP(j, fg_size / 2) {
+            gnx[j * 2] = g[j * 2 + 1];
+            gnx[j * 2 + 1] = g[j * 2];
+        }
+
+        // 分母 g(x)g(-x) の偶数部分
+        constexpr auto half = T{1} / 2;
+        REP(i, fg_size / 2) {
+            g[i] = (
+                  g[i * 2    ] * gnx[i * 2    ]
+                + g[i * 2 + 1] * gnx[i * 2 + 1]
+            ) * half;
+        }
+        g.resize(fg_size / 2, T{0});
+
+        // 2倍ゼロ埋めリサイズ　-> 「2 倍の長さの DFT の計算」
+        auto gt = ntt(g, true);
+        g.resize(fg_size, T{0});
+        REP(i, fg_size / 2) {
+            gt[i] *= ntt_coeff_rev[i];
+        }
+        gt = ntt(move(gt));
+        DSRNG(i, fg_size / 2 - 1, 0) {
+            g[i] = g[i];
+            g[i + fg_size / 2] = gt[i];
+        }
+
+        // 分子 f(x)g(-x) の偶数 or 奇数部分
+        if (m % 2) {
+            REP(i, fg_size / 2) {
+                f[i] = (
+                      f[i * 2    ] * gnx[i * 2    ]
+                    - f[i * 2 + 1] * gnx[i * 2 + 1]
+                ) * half * ntt_coeff_inv[i];
+            }
+        } else {
+            REP(i, fg_size / 2) {
+                f[i] = (
+                      f[i * 2    ] * gnx[i * 2    ]
+                    + f[i * 2 + 1] * gnx[i * 2 + 1]
+                ) * half;
+            }
+        }
+        f.resize(fg_size / 2, T{0});
+
+        // 2倍ゼロ埋めリサイズ　-> 「2 倍の長さの DFT の計算」
+        auto ft = ntt(f, true);
+        f.resize(fg_size, T{0});
+        REP(i, fg_size / 2) {
+            ft[i] *= ntt_coeff_rev[i];
+        }
+        ft = ntt(move(ft));
+        DSRNG(i, fg_size / 2 - 1, 0) {
+            f[i] = f[i];
+            f[i + fg_size / 2] = ft[i];
+        }
+
+        return recurse(recurse, move(f), move(g), m / 2);
+
+    };
+
+    return recurse(recurse, move(f), move(g), m);
 
-    return bostan_mori(move(ab2), move(r), m / 2);
 }