diff --git a/bin/run.py b/bin/run.py
index 629f76d..70e5257 100644
--- a/bin/run.py
+++ b/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",
diff --git a/template/array/sahz.hpp b/template/array/sahz.hpp
index e516330..3f4a12c 100644
--- a/template/array/sahz.hpp
+++ b/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)};
+
+}
\ No newline at end of file
diff --git a/template/conv/ntt.hpp b/template/conv/ntt.hpp
index 457e106..87d432b 100644
--- a/template/conv/ntt.hpp
+++ b/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));
diff --git a/template/fps/berlekamp_massey.hpp b/template/fps/berlekamp_massey.hpp
new file mode 100644
index 0000000..1f980aa
--- /dev/null
+++ b/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;
+}
diff --git a/template/fps/bostan_mori.hpp b/template/fps/bostan_mori.hpp
index e60b68d..faf2629 100644
--- a/template/fps/bostan_mori.hpp
+++ b/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);
 }
diff --git a/template/fps/compose.hpp b/template/fps/compose.hpp
new file mode 100644
index 0000000..c83d459
--- /dev/null
+++ b/template/fps/compose.hpp
@@ -0,0 +1,234 @@
+#pragma once
+#include "base_template.hpp"
+#include "modint/modint.hpp"
+#include "conv/ntt.hpp"
+#include "fps/inv.hpp"
+
+// 多項式合成 f(g(x)) Kinoshita-Li algorithm O(N(logN)^2)
+// verified at https://judge.yosupo.jp/submission/339265
+// https://noshi91.hatenablog.com/entry/2024/03/16/224034
+// https://maspypy.com/fps-%e5%90%88%e6%88%90%e3%83%bb%e9%80%86%e9%96%a2%e6%95%b0%e3%81%ae%e8%a7%a3%e8%aa%ac%ef%bc%882%ef%bc%89%e8%bb%a2%e7%bd%ae%e5%8e%9f%e7%90%86%e3%81%ab%e3%82%88%e3%82%8b%e5%90%88%e6%88%90%e3%82%a2
+// https://arxiv.org/abs/2404.05177
+template <mod_integral T>
+constexpr vector<T> fps_compose(
+    vector<T> f_, vector<T> g_
+) {
+    assert(f_.size() == g_.size());
+    assert(f_.size() != 0);
+    const auto n_original = (ll)f_.size();
+    const auto n_ = (ll)bit_ceil(f_.size());
+    f_.resize(n_, T{0});
+    g_.resize(n_, T{0});
+    const auto m_ = n_;
+    assert(ssize(g_) == n_);
+    assert(g_[0] == 0);
+    const auto m = m_ < 0 ? n_ : m_;
+    assert(m <= n_);
+
+    // Power Projection の勾配を求めると答えが求まる
+    // まずpower projection の重み付き和をとる 
+    // 勾配を求めるための不定元 w(x) を入力の片方としておく
+    // \sum_i f_i [x^{n-1}] w(x) g(x)^i
+    // = \sum_i f_i \sum_j w_{n-1-j} [x^j] g(x)^i
+    // そこでw(x)の係数について勾配をとる
+    //   ∇_w \sum_i f_i \sum_j w_{n-1-j} [x^j] g(x)^i
+    // = ∇_w \sum_i f_i \sum_j w_{j} [x^j] g(x)^i
+    // = ∇_w \sum_j w_j \sum_i [x^{n-1-j}] f_i g(x)^i
+    // = \sum_i [x^{n-1-j}] f_i g(x)^i　(j: 自由変数)
+    // = [x^{n-1-j}] \sum_i f_i g(x)^i
+    // = [x^j] rev(f(g(x)))  Wow!
+
+    auto initialize_forward = [&n_](
+        vector<T> f_,
+        const ll n, const ll l
+    ) {
+        // Q(x, y) = 1 - yf(x)
+        auto q = move(f_);
+        REP(i, ssize(q)) q[i] = -q[i];
+        q.resize(n * l, T{0});
+        q[n] = 1;
+        return q;
+    };
+
+    // ∇_g X
+    auto initialize_backward = [&n_](const ll g_ssize, vector<T> grad_p) {
+        auto grad_g_ = move(grad_p);
+        grad_g_.resize(g_ssize, T{0});
+        return grad_g_;
+    };
+
+    // convolve_p の backward; 便宜上第一引数の勾配を求める
+    auto convolve_backward_p = [](
+        vector<T> q,
+        vector<T> grad_pq
+    ) {
+        // middle product
+        const auto p_ssize = ssize(grad_pq) + 1 - ssize(q);
+        assert(p_ssize > 0);
+        // FFTの勾配を考えて高速化、長さ |p| + |q| - 1 に抑える
+        assert(ssize(q) >= 0);
+        const auto h = (ll)bit_ceil(grad_pq.size());
+
+        auto grad_p = move(grad_pq);
+        grad_p.resize(h, T{0});
+
+        // transposed INTT
+        ranges::reverse(grad_p.begin() + 1, grad_p.end());
+        grad_p = ntt(move(grad_p));
+        // 周波数積
+        q.resize(h, T{0});
+        q = ntt(move(q));
+        RPE(i, ssize(grad_p)) grad_p[i] *= q[i];
+        // transposed NTT
+        grad_p = ntt(move(grad_p), true);
+        ranges::reverse(grad_p.begin() + 1, grad_p.end());
+
+        grad_p.resize(p_ssize, T{0});
+        return grad_p;
+    };
+    
+    // Graeffe’s method: Q(x) Q(-x) で n 半々, l 倍々 にしていく
+    auto graeffe_forward =
+        []
+        (
+            vector<T> q,
+            const ll n,
+            const ll l
+        )
+    {
+        assert(ssize(q) / 2 % 2 == 0);
+        assert(q[q.size() / 2] == 1);
+
+        // 周波数領域に写す
+        q = ntt(move(q));
+
+        // Q(-x): FFT後は前半と後半のswapとなるがbit-reverseしていないので隣接swapでOK
+        auto qnx = vector(q.size(), T{0});
+        REP(j, ssize(q) / 2) {
+            qnx[j * 2] = q[j * 2 + 1];
+            qnx[j * 2 + 1] = q[j * 2];
+        }
+        
+        // 分母 Q(x)Q(-x) の x^偶数次 をとる
+        // 偶数次だけ取れればよいので高速化
+        auto qqnx = move(q);
+        constexpr auto half = T{1} / 2;
+        RPE(i, ssize(qqnx) / 2) {
+            qqnx[i] = (
+                  qqnx[i * 2    ] * qnx[i * 2    ]
+                + qqnx[i * 2 + 1] * qnx[i * 2 + 1]
+            ) * half;
+        }
+        qqnx.resize(ssize(qqnx) / 2, T{0});
+        qqnx = ntt(move(qqnx), true);
+        qqnx[0] -= 1;  // 巡回してきてはみ出た分を調整
+        auto q2 = vector(n * l, T{0});
+        REP(i, n / 4) REP(j, l) {
+            q2[i + (n / 2) * j] = qqnx[i + n / 2 * j];
+        }
+        q2[(n / 2) * l] += 1;  // 巡回してきてはみ出た分を調整
+        
+        return tuple{move(q2), move(qnx), n / 2, l * 2};
+    };
+
+    auto graeffe_backward = [&convolve_backward_p](
+        vector<T> qnx,
+        vector<T> grad_p2,
+        const ll n, const ll l
+    ) -> vector<T> {
+        // b[i + (n / 2) * j] = pqnx[i * 2 + 1 + n * j];
+        auto grad_pqnx = vector(n * l, T{0});
+        // DUMP(grad_p2.size(), n * l);
+        REP(i, n / 4) REP(j, l) {
+            grad_pqnx[i * 2 + 1 + n * j] += grad_p2[i + (n / 2) * j];
+        }
+        
+        // grad of NTTs
+        auto grad_p = move(grad_pqnx);
+        assert(grad_p.size() == qnx.size());
+        
+        // pqnx = ntt(move(pqnx), true);
+        ranges::reverse(grad_p.begin() + 1, grad_p.end());
+        grad_p = ntt(move(grad_p));
+        // 周波数積
+        RPE(i, ssize(grad_p)) grad_p[i] *= qnx[i];
+        // auto pqnx = ntt(move(p));
+        grad_p = ntt(move(grad_p), true);
+        ranges::reverse(grad_p.begin() + 1, grad_p.end());
+        // DUMP(grad_p);
+        return grad_p;
+        
+    };
+
+    auto finalize_forward = [&n_](
+        const vector<T> q_,
+        const ll n
+    ) {
+        const auto m = n_;
+        auto q = vector(m, T{0});
+        RPE(i, m) q[i] = q_[(i + n_ - m + 1) * n];
+        ranges::reverse(q);
+        auto invq = fps_inv(q);
+        return move(invq);
+    };
+    auto finalize_backward = [&convolve_backward_p, &n_](
+        vector<T> invq,
+        vector<T> grad_pq,
+        const ll n
+    ) -> vector<T> {
+        const auto m = n_;
+        auto grad_p = convolve_backward_p(move(invq), move(grad_pq));
+        auto grad_p_ = vector(n_ * n * 2, T{0});
+        RPE(i, n_) grad_p_[(i + n_ - m) * n] = grad_p.rbegin()[i];
+        return grad_p_;
+    };
+
+    // [x^i y^j] arr = arr[i + j * n]
+    // 便宜上 y = x^n として実装を進める
+    // Q(x) の末尾は必ず 1 なので定数倍改善の余地あり
+
+    auto nls = vector<pair<ll, ll>>{{n_ * 2, 2LL}};
+    auto qprev = vector<T>{};
+    auto qnxs = vector<vector<T>>{};
+
+    qprev = initialize_forward(g_, nls.back().first, nls.back().second);
+
+    while (nls.back().first > 2) {
+        auto&& [q2, qnx, n2, l2] = 
+            graeffe_forward(
+                move(qprev),
+                nls.back().first, nls.back().second
+            );
+        qprev = move(q2);
+        qnxs.emplace_back(move(qnx));
+        nls.emplace_back(n2, l2);
+    }
+
+    auto invq = finalize_forward(
+        move(qprev), nls.back().first
+    );
+
+    // 勾配算出
+
+    auto grad_pq = move(f_);
+    grad_pq.resize(ssize(invq) * 2 - 1, T{0});
+
+    auto grad_p_ = finalize_backward(
+        move(invq), move(grad_pq), nls.back().first
+    );
+
+    DSRNG(i, ssize(nls) - 2, 0) {
+        auto&& [n, l] = nls[i];
+        grad_p_ = graeffe_backward(
+            move(qnxs[i]), move(grad_p_), n, l
+        );
+    }
+
+    auto grad_w_ = initialize_backward(n_, move(grad_p_));
+
+    ranges::reverse(grad_w_);
+
+    grad_w_.resize(n_original, T{0});
+    return grad_w_;
+    
+}
diff --git a/template/fps/compositional_inv.hpp b/template/fps/compositional_inv.hpp
new file mode 100644
index 0000000..4e21fc1
--- /dev/null
+++ b/template/fps/compositional_inv.hpp
@@ -0,0 +1,55 @@
+#pragma once
+#include "base_template.hpp"
+#include "modint/modint.hpp"
+#include "conv/ntt.hpp"
+#include "fps/log.hpp"
+#include "fps/exp.hpp"
+#include "fps/power_projection.hpp"
+
+// FPS 合成に対する逆関数 O(N(logN)^2)
+// f に対して g(f(x)) = f(g(x)) なる g を返す
+// verified at https://judge.yosupo.jp/submission/339294
+// https://maspypy.com/fps-%E5%90%88%E6%88%90%E3%83%BB%E9%80%86%E9%96%A2%E6%95%B0%E3%81%AE%E8%A7%A3%E8%AA%AC-1-%E9%80%86%E9%96%A2%E6%95%B0%E3%81%A8-power-projection
+template <mod_integral T>
+constexpr vector<T> compositional_inverse(vector<T> f)
+{
+    assert(ssize(f) >= 2);
+    assert(f[0] == 0);
+    assert(f[1] != 0);
+    if (f[1] != 1) {
+        const auto cinv = f[1].inv();
+        REP(i, ssize(f)) f[i] *= cinv;
+        auto result = compositional_inverse(move(f));
+
+        auto base = T{1};
+        REP(i, ssize(result)) {
+            result[i] *= base;
+            if (i + 1 == i_len) break;
+            base *= cinv;
+        }
+        return result;
+    }
+
+    const auto n_ = ssize(f);
+    f.resize(bit_ceil(f.size()), T{0});
+    const auto n = ssize(f);
+    
+    // [x^{n-1-i}] i * (g(x) / x)^{-n} = n * [x^{n-1}] f(x)^i
+    // [x^{n-1-i}] (g(x) / x)^{-n} = n/i * [x^{n-1}] f(x)^i
+    // ... Lagrange Inversion Theorem
+    const auto one = vector<T>{1};
+    auto goxpmn = power_projection(move(f), move(one), n);
+    RANGE(i, 1 + n - n_, n - 1) goxpmn[i] *= modint(n - 1) / i;
+    ranges::reverse(goxpmn);
+    
+    // g(x) / x
+    goxpmn.resize(n_ - 1);
+    auto gox = fps_log(move(goxpmn));
+    const auto m1on = -T{1} / (n - 1);
+    REP(i, ssize(gox)) gox[i] *= m1on;
+    gox = fps_exp(gox);
+    
+    // g(x)
+    gox.emplace(gox.begin(), T{0});
+    return gox;
+}
\ No newline at end of file
diff --git a/template/fps/fps_rational.hpp b/template/fps/fps_rational.hpp
new file mode 100644
index 0000000..99c29ac
--- /dev/null
+++ b/template/fps/fps_rational.hpp
@@ -0,0 +1,23 @@
+#pragma once
+#include "base_template.hpp"
+#include "modint/modint.hpp"
+#include "fps/berlekamp_massey.hpp"
+
+// P(x)/Q(x) = s(x) mod x^n なる P, Q を求める O(N^2)
+// Q(x) の次数を最小とする
+// C-recursive の初項 s から線形漸化式を求めるのに相当する
+template <mod_integral T>
+pair<vector<T>, vector<T>> fps_rational(vector<T> s) {
+
+    auto denominator = berlekamp_massey(s);
+
+    // 初項 s(x) mod x^#s, 分母 f(x) が求まっているときに分子 A(x) を求める
+    // a(x)/f(x) = s(x) mod x^#s ==> a(x) = s(x)f(x) mod x^#s
+    const auto sz = s.size();
+    auto numerator = convolve_p(move(s), denominator);
+    numerator.resize(sz, T{0});
+    while (!numerator.empty() && numerator.back() == 0) {
+        numerator.pop_back();
+    }
+    return pair{move(numerator), move(denominator)};
+}
diff --git a/template/fps/log.hpp b/template/fps/log.hpp
index 5a3e015..96afbdb 100644
--- a/template/fps/log.hpp
+++ b/template/fps/log.hpp
@@ -10,6 +10,7 @@
 template <mod_integral T>
 constexpr vector<T> fps_log(vector<T> f) {
     assert(f[0] == 1);
+    if (ssize(f) == 1) return vector{T{0}};
     const auto n = f.size();
     auto derf = fps_derivative(f);
     auto invf = fps_inv(move(f));
diff --git a/template/fps/power_projection.hpp b/template/fps/power_projection.hpp
new file mode 100644
index 0000000..06ded6e
--- /dev/null
+++ b/template/fps/power_projection.hpp
@@ -0,0 +1,144 @@
+#pragma once
+#include "base_template.hpp"
+#include "modint/modint.hpp"
+#include "conv/ntt.hpp"
+#include "fps/inv.hpp"
+#include "number/primitive_root.hpp"
+
+// Power Projection O(N(logN)^2)
+// verified at https://judge.yosupo.jp/submission/339294
+// ⚠️ n_ は2冪である必要がある
+// ret[i] := [x^{n-1}] g(x) f(x)^i (i = 0, 1, ..., m-1) を求める
+// 上限 m_ defaults to n.
+// ⚠️ 必要 m_ <= n_
+// https://maspypy.com/fps-%E5%90%88%E6%88%90%E3%83%BB%E9%80%86%E9%96%A2%E6%95%B0%E3%81%AE%E8%A7%A3%E8%AA%AC-1-%E9%80%86%E9%96%A2%E6%95%B0%E3%81%A8-power-projection
+template <mod_integral T>
+constexpr vector<T> power_projection(
+    vector<T> f_, vector<T> g_,
+    const ll n_, const ll m_ = -1
+) {
+    assert(popcount((ull)n_) == 1);
+    const auto m = m_ < 0 ? n_ : m_;
+    assert(m <= n_);
+
+    // [x^i y^j] arr = arr[i + j * n]
+    // 便宜上 y = x^n として実装を進める
+    // Q(x) の末尾は必ず 1 なので定数倍改善の余地あり
+
+    // それぞれ畳み込み用に 2 倍しておく
+    auto n = n_ * 2;  // (x の最高次 + 1) * 2
+    auto l = 2LL; // (y の最高次 + 1) * 2
+
+    
+    // P(x, y) = g(x)
+    auto p = move(g_);
+    p.resize(n * l, T{0});
+        
+    // Q(x, y) = 1 - yf(x)
+    // -(rev_y)-> Q'(x, y) = -f(x) + y
+    auto q = move(f_);
+    REP(i, ssize(q)) q[i] = -q[i];
+    q.resize(n * l, T{0});
+    q[n] = 1;
+
+    // 周波数領域ずらし用前計算 (ω_{2n}^BitReverse(i))_i=0^{n-1}
+    auto ntt_coeff = vector(n * l / 2, T{0});
+    auto ntt_coeff_inv = vector(n * l / 2, T{0});
+    {
+        constexpr auto prim = T{primitive_root(T::pdiv)};
+        const auto root_inv = prim.pow((T::pdiv - 1) / (n * l)).inv();
+        auto coeff_inv = T{1};
+        ntt_coeff_inv[0] = coeff_inv;
+        auto idx = 0LL;
+        REP(i, n * l / 2 - 1) {
+            // bit-reverse
+            auto bit = ~i & (i + 1);
+            auto rev = (n * l / 4) / bit;
+            idx ^= (n * l / 2 - 1) & ~(rev - 1);
+            // non-one coeff
+            ntt_coeff_inv[idx] = (coeff_inv *= root_inv);
+        }
+    }
+    
+    // Graeffe’s method: Q(x) Q(-x) で n 半々, l 倍々 にしていく
+    while (n > 2) {
+        assert(ssize(q) / 2 % 2 == 0);
+        assert(q[q.size() / 2] == 1);
+        // 周波数領域に写す
+        p = ntt(move(p));
+        q = ntt(move(q));
+        // Q(-x): FFT後は前半と後半のswapとなるがbit-reverseしていないので隣接swapでOK
+        auto qnx = vector(q.size(), T{0});
+        REP(j, ssize(q) / 2) {
+            qnx[j * 2] = q[j * 2 + 1];
+            qnx[j * 2 + 1] = q[j * 2];
+        }
+        // 分子 P(x)Q(-x) の x^奇数次 をとる
+        auto pqnx = move(p);
+        auto b = vector<T>{};
+        {
+            constexpr auto half = T{1} / 2;
+            REP(i, ssize(pqnx) / 2) {
+                pqnx[i] = (
+                      pqnx[i * 2    ] * qnx[i * 2    ]
+                    - pqnx[i * 2 + 1] * qnx[i * 2 + 1]
+                ) * half * ntt_coeff_inv[i];
+            }
+            pqnx.resize(ssize(pqnx) / 2, T{0});
+            pqnx = ntt(move(pqnx), true);
+            b = move(pqnx);
+            b.resize(n * l, T{0});
+            REP(j, l) {
+                ranges::fill(
+                    b.begin() + (n / 2) * j + n / 4,
+                    b.begin() + (n / 2) * (j + 1),
+                    T{0}
+                );
+            }
+        }
+        // 分母 Q(x)Q(-x) の x^偶数次 をとる
+        auto qqnx = move(q);
+        auto q2 = vector<T>{};
+        {
+            REP(i, ssize(qqnx) / 2) {
+                constexpr auto half = T{1} / 2;
+                qqnx[i] = (
+                    qqnx[i * 2    ] * qnx[i * 2    ]
+                    + qqnx[i * 2 + 1] * qnx[i * 2 + 1]
+                ) * half;
+            }
+            qqnx.resize(ssize(qqnx) / 2, T{0});
+            qqnx = ntt(move(qqnx), true);
+            qqnx[0] -= 1;  // 巡回してきてはみ出た分を調整
+            q2 = move(qqnx);
+            q2.resize(n * l, T{0});
+            REP(j, l) {
+                ranges::fill(
+                    q2.begin() + (n / 2) * j + n / 4,
+                    q2.begin() + (n / 2) * (j + 1),
+                    T{0}
+                );
+            }
+            q2[q2.size() / 2] += 1;  // 巡回してきてはみ出た分を調整
+        }
+        // n 半々, l 倍々
+        p = move(b);
+        q = move(q2);
+        n /= 2;
+        l *= 2;
+    }
+    
+    assert(l == n_ * 2);
+    RPE(i, m) p[i] = p[(i + n_ - m) * n];
+    RPE(i, m) q[i] = q[(i + n_ - m + 1) * n];
+    p.resize(m, T{0});
+    q.resize(m, T{0});
+    ranges::reverse(p);
+    ranges::reverse(q);
+
+    q = fps_inv(q);
+    p = convolve_p(move(p), move(q));
+    p.resize(m);
+    return p;
+    
+}
\ No newline at end of file
diff --git a/template/modint/modint.hpp b/template/modint/modint.hpp
index 19ccc1d..0b419d7 100644
--- a/template/modint/modint.hpp
+++ b/template/modint/modint.hpp
@@ -131,7 +131,7 @@ istream& operator>> (istream& ist, moduloint<pdiv>& a) { auto x = (typename modu
 // DUMP()対応
 template <int64_t pdiv>
 ostream& __dump_single(const moduloint<pdiv> value) {
-    return cerr << "\e[35m" << value.item() << "\e[2m_" << pdiv << "\e[m";
+    return cerr << "\e[35m" << value.item() << "\e[2mp\e[m";
 }
 
 // ranges_to 対応
diff --git a/template/segtree/bit.hpp b/template/segtree/bit.hpp
index 4112883..7b6db4d 100644
--- a/template/segtree/bit.hpp
+++ b/template/segtree/bit.hpp
@@ -25,7 +25,7 @@ class BIT {
     {
         auto i = size_t{0};
         for (auto&& x : r) {
-            addition(i++, forward<decltype(x)>(x));
+            add(i++, forward<decltype(x)>(x));
         }
     }
 
@@ -51,6 +51,6 @@ class BIT {
     T sum(ptrdiff_t l, ptrdiff_t r) {
         chmax(l, 0);
         chmin(r, static_cast<ptrdiff_t>(n));
-        return addition(get(static_cast<T>(r)), negate(get(static_cast<T>(l))));
+        return addition(get(r), negate(get(l)));
     }
 };
diff --git a/template/segtree/lazy.hpp b/template/segtree/lazy.hpp
index 523575b..b01468a 100644
--- a/template/segtree/lazy.hpp
+++ b/template/segtree/lazy.hpp
@@ -6,9 +6,11 @@
 // Query: 作用素
 // Apply: 作用素をモノイド元に適用 (Query, Value) -> Value
 // Composite: 作用素の合成 (Query, Query) -> Query, (f, g) |-> f○g
+// ⚠️ iterative segment tree 型実装をとっていて、葉が [N+1, 2N+1) となっている
 // [References]
 // - https://hcpc-hokudai.github.io/archive/structure_segtree_001.pdf
 // - https://qiita.com/hotman78/items/b52e0b84ed7082761bad
+// - https://codeforces.com/blog/entry/18051
 template<typename Value, typename Op, typename Query, typename Apply, typename Composite>
 class lazy_segtree {
 
@@ -145,6 +147,116 @@ class lazy_segtree {
 
     }
 
+    // iterative segment tree で [l_, r_) を被覆する O(logN) 個の頂点を返す
+    // 格納順は左から右になる; eval(roots[0]) op ... op eval(roots.end()[-1]) で求まる
+    // time も O(logN) のはず
+    vector<size_t> find_roots(ll l_, ll r_) const {
+        auto [l, r] = regularize_segment(l_, r_);
+
+        auto roots = vector<size_t>{};
+        {
+            auto rleft = vector<size_t>{};
+            auto rright = vector<size_t>{};
+            const auto offset = n + 1;
+            l += offset;
+            r += offset;
+            while (l < r) {
+                if (l & 1) rleft.emplace_back(l++);
+                if (r & 1) rright.emplace_back(--r);
+                l >>= 1;
+                r >>= 1;
+            }
+            // roots = rleft + rright[::-1]
+            roots = move(rleft);
+            ranges::copy(rright | views::reverse, back_inserter(roots));
+        }
+        return roots;
+    }
+
+    // begin/l を固定して end/r を二分探索
+    // predは Value を単一の引数として () 演算子が真偽値を返す必要がある
+    // そのとき、pred(op(v[l:r])) が false となる最大の境界 r を返す
+    // そのような境界が見つからなければ、すなわち op(v[l:]) のとき、無効値を返す
+    // ⚠️ ACL と pred の真偽が逆な点に注意（ACLは単調減少で pred(id) == true）
+    // verified at https://atcoder.jp/contests/abc426/submissions/73129083
+    template <typename Pred>
+    requires predicate<Pred, Value>
+    optional<size_t> bisection_end(ll l_, Pred pred) {
+        const auto l = l_ >= 0 ? (size_t)l_ : size_t{0};
+        assert(l <= n);
+        if (pred(identity_element)) return nullopt;
+        const auto roots = find_roots(l, n);
+        auto sum = identity_element;
+        auto advance_and_add_if = [this, &sum, &pred](
+            auto u, bool negates_prod = false
+        ) -> bool {
+            const auto term = op(sum, eval(u));
+            if (pred(term) == negates_prod) {
+                sum = term;
+                return true;
+            }
+            return false;
+        };
+        eval_topdown(l + n + 1);
+        eval_topdown(n * 2 + 1);
+        for (const auto root : roots) {
+            if (advance_and_add_if(root)) continue;
+            auto u = root;
+            const auto num_non_leaves = n + 1;
+            while (u < num_non_leaves) {
+                eval(u);
+                eval(u <<= 1);
+                if (advance_and_add_if(u)) ++u;
+            }
+            // trueの境界が端に来る場合は右端の要素(!=境界)で止まる
+            if (u < n * 2 + 1) {
+                if (advance_and_add_if(u)) ++u;
+            }
+            return u - num_non_leaves;
+        }
+        return n;
+    }
+
+    // verified at https://atcoder.jp/contests/abc426/submissions/73129132
+    template <typename Pred>
+    requires predicate<Pred, Value>
+    optional<size_t> bisection_begin(ll r_, Pred pred) {
+        assert(r_ >= 0);
+        const auto r = r_ >= 0 ? (size_t)r_ : size_t{0};
+        if (pred(identity_element)) return nullopt;
+        const auto roots = find_roots(0, r);
+        auto sum = identity_element;
+        auto chadd_if = [this, &sum, &pred](
+            const size_t u, bool negates_prod = false
+        ) -> bool {
+            const auto term = op(eval(u), sum);
+            if (pred(term) == negates_prod) {
+                sum = term;
+                return true;
+            }
+            return false;
+        };
+        eval_topdown(n + 1);
+        eval_topdown(r + n + 1);
+        for (const auto root : roots | views::reverse) {
+            if (chadd_if(root)) continue;
+            auto u = root;
+            const auto num_non_leaves = n + 1;
+            while (u < num_non_leaves) {
+                eval(u);
+                eval(u <<= 1);
+                if (++u >= n * 2 + 1) --u;
+                if (chadd_if(u)) --u;
+            }
+            // trueの境界が端に来る場合は右端の要素(!=境界)で止まる
+            if (u >= 0) {
+                if (chadd_if(u)) --u;
+            }
+            return u + 1 - num_non_leaves;
+        }
+        return 0;
+    }
+
     // a[i] を返す
     // 範囲外の場合は単位元を返す
     Value fetch_at(size_t idx) {
@@ -196,4 +308,14 @@ class lazy_segtree {
         return apply_range(idx, idx + 1, f);
     }
 
+    // モノイド (値, op) の単位元を返す
+    Value get_identity_value() const noexcept {
+        return identity_element;
+    }
+
+    // モノイド (作用素, composite) の単位元を返す
+    Query get_identity_query() const noexcept {
+        return identity_query;
+    }
+
 };